## Permission is hereby granted to anyone to copy, use and distribute ## the contents of this file provided that it is not altered in any way. # This version of the code is for use with Maple 6 (including Maple 6.01). # It will not work with previous versions of Maple. if not type(Matrix,type) then S:= kernelopts(version); interface(errorbreak=2); ERROR(`You are running`,substring(S,1.. 7+SearchText(`,`,S,9..-1)),`so use the file gen4.txt`); fi: generator:= proc(Po) # Maple 6 version # Po is a Markov transition matrix # If possible, find and return a generator for P. # Otherwise, find and return best approximate generator # with all entries > -.001. # Otherwise return NULL # If environment variable EnvAllGenerators is true, # return sequence of all approximate generators with # all entries > -.001. # if infolevel[generator] > 0, may return reason # for nonexistence of generator options `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local P, N, odigs; if nargs > 1 then Digits:= args[2] fi; odigs:= Digits; if Digits > evalhf(Digits) then UseHardwareFloats:= false else Digits:= round(evalhf(Digits)) fi; userinfo(2,generator,`Searching for generator with Digits=`||Digits); if type(Po,Matrix) then P:= Po else P:= convert(Po,Matrix) fi; P:= evalf(P); N:= LinearAlgebra:-Dimension(P); # check this is a suitable matrix if nops([N]) <> 2 or N[1] <> N[2] then ERROR(`Expected a square matrix but received`,eval(Po)) fi; if not type(P, 'Matrix'(nonnegative)) then ERROR(`Expected a nonnegative matrix but received`, eval(Po)) fi; `generator/generator`(P,Digits,odigs) end; `generator/generator`:= proc(Po,Di,odi) options `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local P, N, R, i,j, d, L, E, odigs, addigs, k, Kmax, Kmin, fvals, B, Bi, posevs, negevs, marks, c, score, Q, best, diffs, bestj, npairs, a, DD, DK, ke, EQ, qmin, bestQ, res; P:= Po; Digits:= Di; odigs:= odi; N:= (LinearAlgebra:-Dimension(P))[1]; # normalize rows if necessary R:= P . ; for i from 1 to N do if abs(R[i]-1) > .001 then ERROR(`Sum of row`,i,`is not 1 but`,R[i]) fi; LinearAlgebra:-RowOperation(P,i,1/R[i],inplace=true); od; # determinant d:= LinearAlgebra:-Determinant(P); # check necessary conditions if d <= 0 then userinfo(1,generator,`Matrix has non-positive determinant`); RETURN(NULL) elif d >= mul(P[i,i],i=1..N) + Float(1,2-Digits) then userinfo(1,generator, `Determinant exceeds product of diagonal elements`); RETURN(NULL) fi; # check eigenvalues of P E,B:= LinearAlgebra:-Eigenvectors(P); E:= map(proc(x) if abs(Im(x)) < Float(1,2-Digits) then Re(x) else x fi end, convert(E,list)); userinfo(3,generator,`Eigenvalues `,E); diffs:= simplify(fnormal( [seq(seq(E[i]-E[j],j=i+1..N),i=1..N-1)],Digits, Float(1,1-odigs)),zero); if member(0.,diffs) then WARNING("Matrix does not have distinct eigenvalues"); fi; L:= abs(ln(d)); addigs:= 2+round(log10(LinearAlgebra:-ConditionNumber(B))) + round(L/ln(10)); if Digits < odigs + addigs then if odigs + addigs > 100 then userinfo(1,generator,`Eigenvector matrix is singular`); RETURN(NULL); fi; userinfo(2,generator,`Going to`,odigs+addigs,`digits`); RETURN(`generator/generator`(Po,odigs+addigs,odigs)) fi; Bi:= B^(-1); # identify eigenvalues with positive real part, pair with # complex conjugates posevs:= NULL; negevs:= NULL; marks:= [true $ nops(E)]; for i from 1 to nops(E) do if marks[i] and not type(E[i], positive) then c:= conjugate(E[i]); best:= infinity; for j from i+1 to nops(E) do if marks[j] then score:= abs(c-E[j]); if score < best then best:= score; bestj:= j fi fi od; if best > Float(1,3-Digits) then userinfo(1,generator,`Unpaired non-positive eigenvalue`,E[i]); RETURN(NULL) fi; marks[bestj]:= false; if Im(E[i]) >= 0 then posevs:= posevs,i; negevs:= negevs,bestj; if Im(E[bestj])=0 then E[bestj]:= Complex(E[bestj],-0.0) fi; else posevs:= posevs,bestj; negevs:= negevs,i; if Im(E[i])=0 then E[i]:= Complex(E[i],-0.0) fi; fi; fi od; posevs:= [posevs]; negevs:= [negevs]; npairs:= nops(posevs); userinfo(4,generator,`Conjugate pairs`, seq([E[posevs[j]],E[negevs[j]]],j=1..npairs)); # get minimum and maximum k for each eigenvalue Kmax:= [0 $ npairs]; Kmin:= Kmax; for i from 1 to npairs do a:= argument(E[posevs[i]]); Kmax[i]:= trunc((L-a)/(2*Pi)); Kmin[i]:= trunc((-L-a)/(2*Pi)); od; userinfo(4,generator,`Kmax = `,seq(Kmax[i],i=1..npairs)); userinfo(4,generator,`Kmin = `,seq(Kmin[i],i=1..npairs)); best:= -.001; DD:= LinearAlgebra:-DiagonalMatrix(map(ln,E)); DK:= Matrix(N,N); res:= NULL; while `generator/nextk`(k,Kmin,Kmax) do userinfo(3,generator,`k=`,k); for i from 1 to npairs do ke:= evalf(2*Pi*I*k[i]); DK[posevs[i]$2]:= ke; DK[negevs[i]$2]:= -ke; od; Q:= B . (DD+DK) . Bi; Q:= map(fnormal,map(Re,Q),Digits,Float(5,-1-odigs)); qmin:= min(seq(seq(Q[i,j],j=i+1..N),i=1..N-1), seq(seq(Q[j,i],j=i+1..N),i=1..N-1)); if EnvAllGenerators = true then if qmin > -.001 then userinfo(3,generator, `Possible generator with qmin =`,qmin); res:= res, copy(evalf(Q,odigs)) else userinfo(4,generator,`qmin =`,qmin) fi else if qmin >= 0 then userinfo(3,generator,`Found a generator`); RETURN(evalf(Q,odigs)) else if qmin > best then best:= qmin; bestQ:= Q; fi; if qmin > -.001 then userinfo(3,generator, `Approximate generator with qmin =`,qmin) else userinfo(4,generator,`qmin =`,qmin) fi; fi; fi; od; if EnvAllGenerators = true then RETURN(res) fi; WARNING("No completely valid generator found"); if assigned(bestQ) then evalf(bestQ,odigs) else NULL fi; end; `generator/nextk`:= proc(k::evaln, Kmin, Kmax) options `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local i, kl; if not assigned(k) then k:= [0 $ nops(Kmin)]; RETURN(true) fi; if Kmin = [] then RETURN(false) fi; i:= 1; kl:= eval(k); kl[i]:= `generator/nextki`(kl[i]); while kl[i] > Kmax[i] or kl[i] < Kmin[i] do kl[i]:= 0; i:=i+1; if i > nops(kl) then k:= kl; RETURN(false) fi; kl[i]:= `generator/nextki`(kl[i]); od; k:= kl; true end; `generator/nextki`:= proc(k) options `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; if k >= 0 then -1-k else -k fi end;