## 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 V Release 4 to 5.1.
# It will not work with other versions of Maple.  For
# Maple 6 please use the file gen6.txt.

if type(Matrix,type) then
  S:= kernelopts(version);
  interface(errorbreak=2);
  ERROR(`You are running`,substring(S,1..
    SearchText(`,`,S)),`so use the file gen6.txt`);     
fi:

generator:= proc(Po)
# Maple V Release 5 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 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;
  Digits:= max(Digits,round(evalhf(Digits)));
  userinfo(2,generator,`Searching for generator with Digits=`.Digits);
  if type(Po,matrix) then P:= Po
  else P:= evalm(Po)
  fi;
  P:= map(evalf,P);
  N:= linalg[rowdim](P);
# check this is a suitable matrix, and normalize rows if necessary
  if N <> linalg[coldim](P) then 
    ERROR(`Expected a square matrix but received`,eval(Po))
  fi;
  if not type(P, 'matrix'(nonneg)) 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:= linalg[rowdim](P);
# normalize rows if necessary

  R:= evalm(P &* [seq(1,i=1..N)]);  
  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;
  od;
  P:= evalm(linalg[diag](seq(1/R[i],i=1..N)) &* P);
# determinant
  d:= linalg[det](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:= evalf(Eigenvals(P,B));
  E:= fnormal(convert(E,list));
userinfo(3,generator,`Eigenvalues `,E);
  diffs:= fnormal(
     [seq(seq(E[i]-E[j],j=i+1..N),i=1..N-1)],Digits,
      Float(1,1-odigs));
  if member(0.,diffs) then
     WARNING(`Matrix does not have distinct eigenvalues`);
  fi;
  i:= 1;
  while i <= N do
    if has(E[i],I) then 
      B:= linalg[addcol](B,i+1,i,I);
      B:= linalg[mulcol](B,i+1,-2*I);
      B:= linalg[addcol](B,i,i+1,1);
      i:= i+2
    else
      i:= i+1
    fi
  od;  
  L:= abs(ln(d));
  addigs:= traperror(2+round(readlib(log10)(linalg[cond](B)))
   + round(L/ln(10)));
  if addigs = lasterror then
    ERROR(`Eigenvector matrix is singular`)
  fi;
  if Digits < odigs + addigs then
    if odigs + addigs > 100 then
      ERROR(`Eigenvector matrix is singular`)
    fi;
    userinfo(2,generator,`Going to`,odigs+addigs,`digits`);
    RETURN(`generator/generator`(Po,odigs+addigs,odigs))
  fi;
  Bi:= evalm(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 
        ERROR(`Unpaired non-positive eigenvalue`,E[i])
      fi;
      marks[bestj]:= false;
      if Im(E[i]) >= 0 then
        posevs:= posevs,i;
        negevs:= negevs,bestj;
        if Im(E[bestj])=0 then 
          E[bestj]:= E[bestj]-Float(1,-2*Digits)*I
        fi;
      else
        posevs:= posevs,bestj;
        negevs:= negevs,i;
      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:= linalg[diag](op(map(ln,E)));
  DK:=matrix(N,N,0);
  res:= NULL;
  while `generator/nextk`(k,Kmin,Kmax) do
    userinfo(4,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:= evalm(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:= eval(res), evalf(eval(Q),odigs)
      else
         userinfo(4,generator,`qmin =`,qmin)
      fi
    else     
      if qmin >= 0 then
        userinfo(3,generator,`Found a generator`);
        RETURN(evalf(eval(Q),odigs))
      else 
        if qmin > best then
          best:= qmin;
          bestQ:= eval(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(eval(res))
  fi;
  WARNING(`No completely valid generator found`);
  if assigned(bestQ) then evalf(eval(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;  

if not type(WARNING,procedure) then
# Release 4 doesn't have WARNING
  WARNING:= proc(s) 
    if interface(warnlevel) <> 0 then lprint(`Warning,`,s)
    fi end
fi: