"""Original bounded two-state uniformization versus independent analytic transition."""
import math,json
from fractions import Fraction as F

def uniform(q,t,nu,n=20):
 if not isinstance(q,list) or len(q)!=2 or any(not isinstance(r,list) or len(r)!=2 for r in q) or any(type(v)!=int or abs(v)>10 for r in q for v in r) or any(sum(r)!=0 for r in q) or q[0][1]<0 or q[1][0]<0:raise ValueError('two state generator')
 if any(type(v) not in [int,float] or not math.isfinite(v) for v in [t,nu]) or not 0<=t<=1 or not 0<=nu<=4 or nu<max(-q[0][0],-q[1][1]) or type(n)!=int or not 1<=n<=32:raise ValueError('bounded valid time/rate/terms')
 if nu==0:
  assert q==[[0,0],[0,0]];return {'matrix':[[1.,0.],[0.,1.]],'trace':[],'omitted_poisson_mass_float':0.,'analytic_poisson_tail_upper':0.,'rowmass':1.}
 r=[[F(int(i==j))+F(q[i][j],1)/F(nu) for j in range(2)] for i in range(2)];power=[[F(int(i==j)) for j in range(2)] for i in range(2)];lam=nu*t;weight=math.exp(-lam);total=[[0.,0.],[0.,0.]];trace=[];mass=0.
 for k in range(n):
  if k:weight*=lam/k
  for i in range(2):
   for j in range(2):total[i][j]+=weight*float(power[i][j])
  mass+=weight;trace.append({'k':k,'poisson_weight':weight,'R_power':[[str(v) for v in row] for row in power]});power=[[sum(power[i][z]*r[z][j] for z in range(2)) for j in range(2)] for i in range(2)]
 # First omitted term and geometric majorant when subsequent ratios <1.
 first=weight*lam/n;ratio=lam/(n+1);bound=first/(1-ratio) if ratio<1 else 1.
 return {'R':[[str(v) for v in row] for row in r],'matrix':total,'trace':trace,'rowmass':mass,'omitted_poisson_mass_float':max(0.,1-mass),'analytic_poisson_tail_upper':bound}
def analytic(t):
 a=2/3*(1-math.exp(-3*t));b=1/3*(1-math.exp(-3*t));return [[1-a,a],[b,1-b]]
def demo():
 q=[[-2,2],[1,-1]];cases=[]
 for t in [0.,.5,1.]:
  for nu in [2,3]:
   result=uniform(q,t,nu);exact=analytic(t);err=max(abs(result['matrix'][i][j]-exact[i][j]) for i in range(2) for j in range(2));assert err<=result['analytic_poisson_tail_upper']+3e-15;assert all(abs(sum(row)-result['rowmass'])<1e-14 for row in result['matrix']);cases.append({'t':t,'nu':nu,'absolute_entry_error':err,'analytic':exact,**result})
 sample=cases[3];assert sample['R']==[['1/3','2/3'],['1/3','2/3']]
 short=uniform(q,.5,3,1);assert short['rowmass']<1 and short['matrix']==[[math.exp(-1.5),0.],[0.,math.exp(-1.5)]]
 zero=uniform([[0,0],[0,0]],1,0);assert zero['matrix']==[[1.,0.],[0.,1.]]
 rejects=[]
 for label,qq,t,nu,n in [('bad_rate',q,.5,1,20),('negative_offdiag',[[1,-1],[1,-1]],.5,3,20),('row_sum',[[0,2],[1,-1]],.5,3,20),('negative_time',q,-1,3,20),('nonfinite_time',q,math.inf,3,20),('zero_nu_nonzero_Q',q,.5,0,20),('term_budget',q,.5,3,33),('bool_time',q,True,3,20)]:
  try:uniform(qq,t,nu,n)
  except ValueError:rejects.append(label)
  else:raise AssertionError(label)
 return {'Q':q,'six_rate_time_checks':cases,'one_term_not_renormalized':short,'allzero_generator_policy':zero,'invalid_cases_rejected':rejects,'scope':'Poisson truncation bound in exact-real theory is distinct from float roundoff; +3e-15 is observed test tolerance, not universal guarantee. 2states, no SciPy expm/queue/customer rate fit. Common tick may be a dummy self-transition, not real jump.'}
if __name__=='__main__':print(json.dumps(demo(),indent=2))
