src/pamp.c

   1 /* PAMP.c, Copyright, Ziheng Yang, April 1995.\r
   2    Specify the sequence type in the file pamp.ctl.  Results go into mp.\r
   3 \r
   4                     gcc -o pamp pamp.c tools.o\r
   5                     pamp <ControlFileName>\r
   6 */\r
   7 \r
   8 #include "paml.h"\r
   9 \r
  10 #define NS            2000\r
  11 #define NBRANCH       (NS*2-2)\r
  12 #define NNODE         (NS*2-1)\r
  13 #define MAXNSONS      10\r
  14 #define NGENE         2000\r
  15 #define LSPNAME       30\r
  16 #define NCODE         20\r
  17 #define NCATG         16\r
  18 \r
  19 double DistanceREV (double Ft[], int n, double alpha, double Root[], double U[],\r
  20    double V[], double pi[], double space[], int *cond);\r
  21 int PMatBranch (double Ptb[], int n, double branch[], \r
  22     double Root[], double U[], double V[], double space[]);\r
  23 int PatternLS (FILE *fout, double Ft[],double alpha, double space[], int *cond);\r
  24 int testx (double x[], int np);\r
  25 int GetOptions (char *ctlf);\r
  26 int AlphaMP (FILE* fout);\r
  27 int PatternMP (FILE *fout, double Ft[]);\r
  28 int PathwayMP1 (FILE *fout, int *maxchange, int NSiteChange[], \r
  29     double Ft[], double space[], int job);\r
  30 double lfunAlpha_Sullivan (double x);\r
  31 double lfunAlpha_YK96 (double x);\r
  32 \r
  33 struct CommonInfo {\r
  34    unsigned char *z[NS];\r
  35    char *spname[NS], seqf[256],outf[256],treef[256];\r
  36    int seqtype, ns, ls, ngene, posG[NGENE+1],lgene[NGENE],*pose,npatt, readpattern;\r
  37    int np, ntime, ncode,fix_kappa,fix_rgene,fix_alpha, clock, model, ncatG, cleandata;\r
  38    int print, nhomo;\r
  39    double *fpatt, *conP;\r
  40    /* not used */\r
  41    double lmax,pi[NCODE], kappa,alpha,rou, rgene[NGENE],piG[NGENE][NCODE];\r
  42 }  com;\r
  43 struct TREEB {\r
  44    int nbranch, nnode, root, branches[NBRANCH][2];\r
  45    double lnL;\r
  46 }  tree;\r
  47 struct TREEN {\r
  48    int father, nson, sons[MAXNSONS], ibranch;\r
  49    double branch, age, label, *conP;\r
  50    char *nodeStr, fossil;\r
  51 }  *nodes;\r
  52 \r
  53 \r
  54 #define NCATCHANGE 100\r
  55 extern int noisy, *ancestor;\r
  56 extern double *SeqDistance;\r
  57 int maxchange, NSiteChange[NCATCHANGE];\r
  58 double MuChange;\r
  59 int LASTROUND=0; /* no use for this */\r
  60 \r
  61 #define LSDISTANCE\r
  62 #define REALSEQUENCE\r
  63 #define NODESTRUCTURE\r
  64 #define RECONSTRUCTION\r
  65 #define PAMP\r
  66 #include "treesub.c"\r
  67 \r
  68 int main (int argc, char *argv[])\r
  69 {\r
  70    FILE *ftree, *fout, *fseq;\r
  71    char ctlf[32]="pamp.ctl";\r
  72    char *Seqstr[]={"nucleotide", "", "amino-acid", "Binary"};\r
  73    int itree, ntree, i, j, s3;\r
  74    double *space, *Ft;\r
  75 \r
  76    com.nhomo=1;  com.print=1;\r
  77    noisy=2;  com.ncatG=8;   com.clock=0; com.cleandata=1;\r
  78    starttimer();\r
  79    GetOptions(ctlf);\r
  80    if(argc>1) { strcpy(ctlf, argv[1]); printf("\nctlfile set to %s.\n",ctlf);}\r
  81 \r
  82    printf("PAMP in %s\n", pamlVerStr);\r
  83    if ((fseq=fopen(com.seqf, "r"))==NULL) error2 ("seqfile err.");\r
  84    if ((fout=fopen (com.outf, "w"))==NULL) error2("outfile creation err.");\r
  85    if((fseq=fopen (com.seqf,"r"))==NULL)  error2("No sequence file!");\r
  86    ReadSeq (NULL, fseq, com.cleandata, 0);\r
  87    SetMapAmbiguity();\r
  88    i=(com.ns*2-1)*sizeof(struct TREEN);\r
  89    if((nodes=(struct TREEN*)malloc(i))==NULL) error2("oom");\r
  90 \r
  91    fprintf (fout,"PAMP %15s, %s sequences\n", com.seqf, Seqstr[com.seqtype]);\r
  92    if (com.nhomo) fprintf (fout, "nonhomogeneous model\n");\r
  93 \r
  94    space = (double*)malloc(1000000*sizeof(double));  /* not safe */\r
  95    SeqDistance=(double*)malloc(com.ns*(com.ns-1)/2*sizeof(double));\r
  96    ancestor=(int*)malloc(com.ns*(com.ns-1)/2*sizeof(int));\r
  97    if (SeqDistance==NULL||ancestor==NULL) error2("oom");\r
  98 \r
  99    i = com.ns*(com.ns-1)/2;\r
 100    s3 = sizeof(double)*((com.ns*2-2)*(com.ns*2-2 + 4 + i) + i);\r
 101    s3 = max2(s3, com.ncode*com.ncode*(2*com.ns-2+1)*(int)sizeof(double));\r
 102 \r
 103    Ft = (double*) malloc(s3);\r
 104    if (space==NULL || Ft==NULL)  error2 ("oom space");\r
 105 \r
 106    InitializeBaseAA (fout);\r
 107    if (com.ngene>1) error2 ("option G not allowed yet");\r
 108 \r
 109 /*\r
 110    PatternLS (fout, Ft, 0., space, &i);\r
 111    printf ("\nPairwise estimation of rate matrix done..\n");\r
 112    fflush(fout);\r
 113 */\r
 114    ftree=gfopen (com.treef,"r");\r
 115    fscanf (ftree, "%d%d", &i, &ntree);\r
 116    if (i!=com.ns) error2 ("ns in the tree file");\r
 117 \r
 118    for(itree=0; itree<ntree; itree++) {\r
 119 \r
 120       printf ("\nTREE # %2d\n", itree+1);\r
 121       fprintf (fout,"\nTREE # %2d\n", itree+1);\r
 122 \r
 123       if (ReadTreeN (ftree, &i,&j, 0, 1)) error2 ("err tree..");\r
 124       OutTreeN (F0, 0, 0);    FPN (F0); \r
 125       OutTreeN (fout, 0, 0);  FPN (fout);\r
 126 \r
 127       for (i=0,maxchange=0; i<NCATCHANGE; i++) NSiteChange[i]=0;\r
 128 \r
 129       PathwayMP1 (fout, &maxchange, NSiteChange, Ft, space, 0);\r
 130       printf ("\nHartigan reconstruction done..\n");\r
 131 \r
 132       fprintf (fout, "\n\n(1) Branch lengths and substitution pattern\n");\r
 133       PatternMP (fout, Ft);\r
 134       printf ("pattern done..\n");    fflush(fout);\r
 135 \r
 136       fprintf (fout, "\n\n(2) Gamma parameter\n");\r
 137       AlphaMP (fout);\r
 138       printf ("gamma done..\n");      fflush(fout);\r
 139 \r
 140       fprintf (fout, "\n\n(3) Parsimony reconstructions\n");\r
 141       PathwayMP1 (fout, &maxchange, NSiteChange, Ft, space, 1);\r
 142       printf ("Yang reconstruction done..\n");    fflush(fout);\r
 143    }\r
 144    free(nodes);\r
 145    return (0);\r
 146 }\r
 147 \r
 148 int GetOptions (char *ctlf)\r
 149 {\r
 150    int iopt, nopt=6, i, lline=4096, t;\r
 151    char line[4096], *pline, opt[32], *comment="*#";\r
 152    char *optstr[] = {"seqfile","outfile","treefile", "seqtype", "ncatG", "nhomo"};\r
 153    FILE  *fctl=gfopen (ctlf, "r");\r
 154 \r
 155    if (fctl) {\r
 156       for (;;) {\r
 157          if (fgets (line, lline, fctl) == NULL) break;\r
 158          for (i=0,t=0; i<lline&&line[i]; i++)\r
 159             if (isalnum(line[i]))  { t=1; break; }\r
 160             else if (strchr(comment,line[i])) break;\r
 161          if (t==0) continue;\r
 162          sscanf (line, "%s%*s%d", opt, &t);\r
 163          if ((pline=strstr(line, "="))==NULL) error2 ("option file.");\r
 164 \r
 165          for (iopt=0; iopt<nopt; iopt++) {\r
 166             if (strncmp(opt, optstr[iopt], 8)==0)  {\r
 167                if (noisy>2)\r
 168                   printf ("\n%3d %15s | %-20s %6d", iopt+1,optstr[iopt],opt,t);\r
 169                switch (iopt) {\r
 170                   case ( 0): sscanf(pline+2, "%s", com.seqf);    break;\r
 171                   case ( 1): sscanf(pline+2, "%s", com.outf);    break;\r
 172                   case ( 2): sscanf(pline+2, "%s", com.treef);    break;\r
 173                   case  (3): com.seqtype=t;   break;\r
 174                   case  (4): com.ncatG=t;     break;\r
 175                   case  (5): com.nhomo=t;     break;\r
 176                }\r
 177                break;\r
 178             }\r
 179          }\r
 180          if (iopt==nopt)\r
 181             { printf ("\nopt %s in %s\n", opt, ctlf);  exit (-1); }\r
 182       }\r
 183       fclose (fctl);\r
 184    }\r
 185    else\r
 186       if (noisy) printf ("\nno ctl file..");\r
 187 \r
 188    if (com.seqtype==0)       com.ncode=4;\r
 189    else if (com.seqtype==2)  com.ncode=20;\r
 190    else if (com.seqtype==3)  com.ncode=2;\r
 191    else                      error2("seqtype");\r
 192    if (com.ncatG>NCATG) error2 ("raise NCATG?");\r
 193    return (0);\r
 194 }\r
 195 \r
 196 \r
 197 int AlphaMP (FILE* fout)\r
 198 {\r
 199    int k, ntotal;\r
 200    double x, xb[2], lnL, var;\r
 201 \r
 202    xb[0]=1e-3; xb[1]=99; /* alpha */\r
 203  \r
 204    fprintf (fout, "\n# changes .. # sites");\r
 205    for (k=0,ntotal=0,MuChange=var=0; k<maxchange+1; k++) {\r
 206       fprintf (fout, "\n%6d%10d", k, NSiteChange[k]);\r
 207       ntotal+=NSiteChange[k];  MuChange+=k*NSiteChange[k];   \r
 208       var+=k*k*NSiteChange[k];\r
 209    }\r
 210    MuChange/=ntotal;   \r
 211    var=(var-MuChange*MuChange*ntotal)/(ntotal-1.);\r
 212    x=MuChange*MuChange/(var-MuChange);\r
 213    fprintf (fout, "\n\n# sites%6d,  total changes%6d\nmean-var%9.4f%9.4f",\r
 214             ntotal, (int)(ntotal*MuChange+.5), MuChange, var);\r
 215    fprintf (fout, "\nalpha (method of moments)%9.4f", x);\r
 216    if (x<=0) x=9;\r
 217 \r
 218    LineSearch(lfunAlpha_Sullivan, &lnL, &x, xb, 0.02, 1e-8);\r
 219    fprintf (fout, "\nalpha (Sullivan et al. 1995)%9.4f\n", x);\r
 220 \r
 221    MuChange/=tree.nbranch; \r
 222    LineSearch(lfunAlpha_YK96, &lnL, &x, xb, 0.02, 1e-8);\r
 223    fprintf (fout, "alpha (Yang & Kumar 1995, ncatG= %d)%9.4f\n", com.ncatG,x);\r
 224    return (0);\r
 225 }\r
 226 \r
 227 double lfunAlpha_Sullivan (double x)\r
 228 {\r
 229    int k;\r
 230    double lnL=0, a=x, t;\r
 231 \r
 232    FOR (k, maxchange+1) { \r
 233       if (NSiteChange[k]==0) continue;\r
 234       t=-a*log(1+MuChange/a);\r
 235       if (k)  \r
 236          t+=LnGamma(k+a)-LnGamma(k+1.) - LnGamma(a) \r
 237           + k*log(MuChange/a/(1+MuChange/a));\r
 238       lnL += NSiteChange[k]*t;\r
 239    }\r
 240    return(-lnL);\r
 241 }\r
 242 \r
 243 double lfunAlpha_YK96 (double x)\r
 244 {\r
 245    int k, ir, b=tree.nbranch, n=com.ncode;\r
 246    double lnL=0, prob, a=x, t=MuChange, p;\r
 247    double freqK[NCATG], rK[NCATG];\r
 248 \r
 249    DiscreteGamma (freqK, rK, a, a, com.ncatG, 0);\r
 250    FOR (k, maxchange+1) {\r
 251       if (NSiteChange[k]==0) continue;\r
 252       for (ir=0,prob=0; ir<com.ncatG; ir++) {\r
 253          p=1./n+(n-1.)/n*exp(-n/(n-1.)*rK[ir]*t);\r
 254          prob+=freqK[ir]*pow(p,(double)(b-k))*pow((1-p)/(n-1.),(double)k);\r
 255       }\r
 256       lnL += NSiteChange[k]*log(prob);\r
 257    }\r
 258    return (-lnL);\r
 259 }\r
 260 \r
 261 \r
 262 int OutQ (FILE *fout, int n, double Q[], double pi[], double Root[],\r
 263     double U[], double V[], double space[])\r
 264 {\r
 265    char aa3[4]="";\r
 266    int i,j;\r
 267    double *T1=space, t;\r
 268 \r
 269    fprintf(fout,"\nrate matrix Q: Qij*dt = prob(i->j; dt)\n");\r
 270    if (n<=4) {\r
 271 /*      matout (fout, pi, 1, n); */\r
 272       matout (fout, Q, n, n);\r
 273       if (n==4) {\r
 274          fprintf (fout, "Order: T, C, A, G");\r
 275          t=pi[0]*Q[0*4+1]+pi[1]*Q[1*4+0]+pi[2]*Q[2*4+3]+pi[3]*Q[3*4+2];\r
 276          fprintf (fout, "\nAverage Ts/Tv =%9.4f\n", t/(1-t));\r
 277       }\r
 278    }\r
 279    else if (n==20) {\r
 280       for (i=0; i<n; i++,FPN(fout))\r
 281          FOR (j,n) fprintf (fout, "%6.0f", Q[i*n+j]*100);\r
 282 /*\r
 283       FOR (i,n) {\r
 284          fprintf (fout,"\n%-4s", getAAstr(aa3,i));\r
 285          FOR (j,i) fprintf (fout, "%4.0f", Q[i*n+j]/pi[j]*100);\r
 286          fprintf (fout, "%4.0f", -Q[i*n+i]*100);\r
 287       }\r
 288       fputs("\n     ",fout);  FOR(i,naa) fprintf(fout,"%5s",getAAstr(aa3,i));\r
 289 */\r
 290       fprintf (fout, "\n\nPAM matrix, P(0.01)\n"); \r
 291       FOR (i,n) FOR (j,n) T1[i*n+j]=U[i*n+j]*exp(0.01*Root[j]);\r
 292       matby (T1, V, Q, n, n, n);\r
 293       FOR (i,n*n) if (Q[i]<0) Q[i]=0;\r
 294       FOR (i,n) {\r
 295          fprintf (fout,"\n%-4s", getAAstr(aa3,i));\r
 296          FOR(j,n) fprintf(fout, "%6.0f", Q[i*n+j]*10000);\r
 297       }\r
 298       fputs("\n     ",fout);  FOR(i,n) fprintf(fout,"%5s",getAAstr(aa3,i));\r
 299    }\r
 300    return (0);\r
 301 }\r
 302 \r
 303 int PMatBranch (double Ptb[], int n, double branch[], \r
 304     double Root[], double U[], double V[], double space[])\r
 305 {\r
 306 /* homogeneised transition prob matrix, with one Q assumed for the whole tree\r
 307 */\r
 308    int i, j, k;\r
 309    double *T1=space, *P;\r
 310 \r
 311    FOR (k, tree.nbranch) {\r
 312       P=Ptb+k*n*n;\r
 313       FOR (i,n) FOR (j,n) T1[i*n+j]=U[i*n+j]*exp(Root[j]*branch[k]);\r
 314       matby (T1, V, P, n, n, n);\r
 315       FOR (i,n*n) if (P[i]<0) P[i]=0;\r
 316 /*\r
 317       printf ("\nbranch %d, P(%.5f)", k+1, branch[k]);\r
 318       matout (F0, P, n, n);\r
 319       testTransP (P, n);\r
 320 */\r
 321    }\r
 322    return (0);\r
 323 }\r
 324 \r
 325 \r
 326 int PatternMP (FILE *fout, double Ft[])\r
 327 {\r
 328 /* Ft[]: input counts for the F(t) matrix for each branch, output P(t) \r
 329 */\r
 330    int n=com.ncode, i,j,k;\r
 331    double *Q, *pi, *Root, *U, *V, *branch, *space, *T1, t;\r
 332 \r
 333    if((Q=(double*)malloc((n*n*6+tree.nbranch)*sizeof(double)))==NULL)\r
 334       error2("PathwayMP: oom");\r
 335    pi=Q+n*n; Root=pi+n; U=Root+n; V=U+n*n; branch=V+n*n; \r
 336    space=T1=branch+tree.nbranch;\r
 337 \r
 338    for (k=0; k<tree.nbranch; k++) {  /* branch lengths */\r
 339       xtoy(Ft+k*n*n, Q, n*n);\r
 340       branch[k]=nodes[tree.branches[k][1]].branch=\r
 341          DistanceREV(Q, n, 0, Root, U, V, pi, space, &j);\r
 342    }\r
 343    OutTreeB (fout);  FPN (fout);\r
 344    FOR (i, tree.nbranch) fprintf(fout,"%9.5f", branch[i]);\r
 345    fprintf (fout,"\ntree length: %9.5f\n", sum(branch,tree.nbranch));\r
 346 \r
 347    /* pattern Q from average F(t) */\r
 348    fprintf(fout,"\nF(t)");\r
 349    xtoy (Ft+tree.nbranch*n*n, Q, n*n);\r
 350    matout2 (fout, Q, n, n, 12, 2);\r
 351    DistanceREV (Q, n, 0, Root, U, V, pi, space, &j);\r
 352    if (noisy>=3&&j==-1) { puts("F(t) modified in DistanceREV"); }\r
 353 \r
 354    OutQ (fout, n, Q, pi, Root, U, V, T1);\r
 355    if (com.nhomo==0) \r
 356       PMatBranch (Ft, n, branch, Root, U, V, space);\r
 357    else {\r
 358       for (k=0; k<tree.nbranch; k++) {\r
 359          for (i=0; i<n; i++) {\r
 360             t=sum(Ft+k*n*n+i*n, n);\r
 361             if (t>1e-5) abyx (1/t, Ft+k*n*n+i*n, n);\r
 362             else        Ft[k*n*n+i*n+i]=1;\r
 363          }\r
 364       }\r
 365    }\r
 366    free(Q);\r
 367    return (0); \r
 368 }\r
 369 \r
 370 \r
 371 int PathwayMP1 (FILE *fout, int *maxchange, int NSiteChange[], \r
 372     double Ft[], double space[], int job)\r
 373 {\r
 374 /* Hartigan, JA.  1973.  Minimum mutation fits to a given tree. \r
 375    Biometrics, 29:53-65.\r
 376    Yang, Z.  1996.  \r
 377    job=0: 1st pass: calculates maxchange, NSiteChange[], and Ft[]\r
 378    job=1: 2nd pass: reconstructs ancestral character states (->fout)\r
 379 */\r
 380    char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs));\r
 381    char *zz[NNODE],nodeb[NNODE],bestPath[NNODE-NS],Equivoc[NS-1];\r
 382    int n=com.ncode, nid=tree.nbranch-com.ns+1, it,i1,i2, i,j,k, h, hp,npath;\r
 383    int *Ftt=NULL, nchange, nchange0, visit[NS-1]={0};\r
 384    double sumpr, bestpr, pr, *pnode=NULL, *pnsite;\r
 385 \r
 386 \r
 387    fputs("\nList of most parsimonious reconstructions (MPRs) at each site: #MPRs (#changes)\n",fout);\r
 388    fputs("and then the most likely reconstruction out of the MPRs and its probability\n",fout);\r
 389    if((pnsite=(double*)malloc((com.ns-1)*n*sizeof(double)))==NULL)\r
 390       error2("PathwayMP1: oom");\r
 391 \r
 392    PATHWay=(char*)malloc(nid*(n+3)*sizeof(char));\r
 393    NCharaCur=PATHWay+nid;  ICharaCur=NCharaCur+nid;  CharaCur=ICharaCur+nid;\r
 394    if (job==0) {\r
 395       zero(Ft,n*n*(tree.nbranch+1));\r
 396       if((Ftt=(int*)malloc(n*n*tree.nbranch*sizeof(int)))==NULL) error2("oom");\r
 397    }\r
 398    else {\r
 399       pnode=(double*)malloc((nid*com.npatt+1)*(sizeof(double)+sizeof(char)));\r
 400       FOR (j,nid) zz[com.ns+j]=(char*)(pnode+nid*com.npatt)+j*com.npatt;\r
 401       FOR (j,com.ns) zz[j]=com.z[j];\r
 402       if (pnode==NULL) error2 ("oom");\r
 403    }\r
 404    for (j=0,visit[i=0]=tree.root-com.ns; j<tree.nbranch; j++) \r
 405       if (tree.branches[j][1]>=com.ns) \r
 406          visit[++i]=tree.branches[j][1]-com.ns;\r
 407 \r
 408    for (h=0; h<com.ls; h++) {\r
 409       hp=com.pose[h];\r
 410       if (job==1) {\r
 411          fprintf (fout, "\n%4d  ", h+1);\r
 412          FOR (j, com.ns) fprintf (fout, "%c", pch[com.z[j][hp]]);\r
 413          fprintf (fout, ":  ");\r
 414          FOR (j,nid*n) pnsite[j]=0;\r
 415       }\r
 416       FOR (j,com.ns) nodeb[j]=com.z[j][hp];\r
 417       if (job==0) FOR (j,n*n*tree.nbranch) Ftt[j]=0;\r
 418 \r
 419       InteriorStatesMP (1, hp, &nchange, NCharaCur, CharaCur, space); \r
 420       ICharaCur[j=tree.root-com.ns]=0;  PATHWay[j]=CharaCur[j*n+0];\r
 421       FOR (j,nid) Equivoc[j]=(NCharaCur[j]>1);\r
 422 \r
 423       if (nchange>*maxchange) *maxchange=nchange;\r
 424       if (nchange>NCATCHANGE-1) error2 ("raise NCATCHANGE");\r
 425 \r
 426       NSiteChange[nchange]++;\r
 427       /* NSiteChange[nchange]+=(int)com.fpatt[hp]; */\r
 428 \r
 429       DownStates (tree.root);\r
 430       for (npath=0,sumpr=bestpr=0; ;) {\r
 431          for (j=0,k=visit[nid-1]; j<NCharaCur[k]; j++) {\r
 432             PATHWay[k]=CharaCur[k*n+j]; npath++;\r
 433             FOR (i,nid) nodeb[i+com.ns]=PATHWay[i];\r
 434             if (job==1) {\r
 435                FOR (i,nid) fprintf(fout,"%c",pch[PATHWay[i]]); fputc(' ',fout);\r
 436                pr=com.pi[(int)nodeb[tree.root]];\r
 437                for (i=0; i<tree.nbranch; i++) {\r
 438                   i1=nodeb[tree.branches[i][0]];\r
 439                   i2=nodeb[tree.branches[i][1]];\r
 440                   pr*=Ft[i*n*n+i1*n+i2];\r
 441                }\r
 442                sumpr+=pr;\r
 443                FOR (i,nid) pnsite[i*n+nodeb[i+com.ns]]+=pr;\r
 444                if (pr>bestpr) \r
 445                   { bestpr=pr; FOR(i,nid) bestPath[i]=PATHWay[i];}\r
 446             }\r
 447             else {\r
 448                for (i=0,nchange0=0; i<tree.nbranch; i++) {\r
 449                   i1=nodeb[tree.branches[i][0]]; \r
 450                   i2=nodeb[tree.branches[i][1]];\r
 451                   if(i1!=i2) nchange0++;\r
 452                   Ftt[i*n*n+i1*n+i2]++;\r
 453                }\r
 454                if (nchange0!=nchange) {\r
 455                   printf("\a\nerr:PathwayMP %d != %d", nchange, nchange0); \r
 456                   fprintf(fout,".%d. ", nchange0); /* ??? */\r
 457                }\r
 458             }\r
 459          }\r
 460          for (j=nid-2; j>=0; j--) {\r
 461             if(Equivoc[k=visit[j]] == 0) continue;\r
 462             if (ICharaCur[k]+1<NCharaCur[k]) {\r
 463                PATHWay[k] = CharaCur[k*n + (++ICharaCur[k])];\r
 464                DownStates (k+com.ns);\r
 465                break;\r
 466             }\r
 467             else { /* if (next equivocal node is not ancestor) update node k */\r
 468                for (i=j-1; i>=0; i--) if (Equivoc[(int)visit[i]]) break;\r
 469                if (i>=0) { \r
 470                   for (it=k+com.ns,i=visit[i]+com.ns; ; it=nodes[it].father)\r
 471                      if (it==tree.root || nodes[it].father==i) break;\r
 472                   if (it==tree.root)\r
 473                      DownStatesOneNode (k+com.ns, nodes[k+com.ns].father);\r
 474                }\r
 475             }\r
 476          }\r
 477          if (j<0) break;\r
 478       }      /* for (npath) */\r
 479 /*\r
 480       printf ("\rsite pattern %4d/%4d: %6d%6d", hp+1,com.npatt,npath,nchange);\r
 481 */      \r
 482       if (job==0) \r
 483          FOR (j,n*n*tree.nbranch) Ft[j]+=(double)Ftt[j]/npath*com.fpatt[hp];\r
 484       else {\r
 485          FOR (i,nid) zz[com.ns+i][hp]=bestPath[i];\r
 486          FOR (i,nid) pnode[i*com.npatt+hp]=pnsite[i*n+bestPath[i]]/sumpr;\r
 487          fprintf (fout, " |%4d (%d) | ", npath, nchange);\r
 488          if (npath>1) {\r
 489             FOR (i,nid) fprintf (fout, "%c", pch[bestPath[i]]);\r
 490             fprintf (fout, " (%.3f)", bestpr/sumpr);\r
 491 \r
 492          }\r
 493       }\r
 494    }   /* for (h) */\r
 495    free(PATHWay); \r
 496    if (job==0) {\r
 497       free(Ftt);\r
 498       FOR (i,tree.nbranch) FOR (j,n*n) Ft[tree.nbranch*n*n+j]+=Ft[i*n*n+j];\r
 499    }\r
 500    else {\r
 501       fprintf (fout,"\n\nApprox. relative accuracy at each node, by site\n");\r
 502       FOR (h, com.ls) {\r
 503          hp=com.pose[h];\r
 504          fprintf (fout,"\n%4d  ", h+1);\r
 505          FOR (j, com.ns) fprintf (fout, "%c", pch[com.z[j][hp]]);\r
 506          fprintf (fout, ":  ");\r
 507          FOR (i,nid) if (pnode[i*com.npatt+hp]<.99999) break;\r
 508          if (i<nid)  FOR (j, nid) \r
 509             fprintf(fout,"%c (%5.3f) ", pch[zz[j][hp]],pnode[j*com.npatt+hp]);\r
 510       }\r
 511       /* Site2Pattern (fout); */\r
 512       fprintf (fout,"\n\nlist of extant and reconstructed sequences\n\n");\r
 513       for(j=0;j<tree.nnode;j++,FPN(fout)) {\r
 514          if(j<com.ns) fprintf(fout,"%-20s", com.spname[j]);\r
 515          else         fprintf(fout,"node #%-14d", j+1);\r
 516          print1seq (fout, zz[j], (com.readpattern?com.npatt:com.ls), com.pose);\r
 517       }\r
 518       free(pnode);\r
 519    }\r
 520    free(pnsite);\r
 521    return (0);\r
 522 }\r
 523 \r
 524 double DistanceREV (double Ft[], int n,double alpha,double Root[],double U[],\r
 525    double V[], double pi[], double space[], int *cond)\r
 526 {\r
 527 /* input:  Ft, n, alpha\r
 528    output: Q(in Ft), t, Root, U, V, and cond\r
 529    space[n*n*2]\r
 530 */\r
 531    int i,j, InApplicable;\r
 532    double *Q=Ft, *T1=space, *T2=space+n*n, t, pi_sqrt[20], small=0.1/com.ls;\r
 533 \r
 534    for (i=0,t=0; i<n; i++) FOR (j,n) if (i-j) t+=Q[i*n+j];\r
 535    if (t<small)  { *cond=1; zero(Q,n*n); return (0); }\r
 536 \r
 537    for(i=0;i<n;i++) for (j=0;j<i;j++) \r
 538       Q[i*n+j]=Q[j*n+i]=(Q[i*n+j]+Q[j*n+i])/2;\r
 539 \r
 540    abyx(1./sum(Q,n*n), Q, n*n);\r
 541    for(i=0;i<n;i++) {\r
 542       pi[i]=sum(Q+i*n, n);\r
 543       if(pi[i]>small) \r
 544          abyx(1/pi[i], Q+i*n, n); \r
 545    }\r
 546 \r
 547    eigenQREV(Q, pi, n, Root, U, V, pi_sqrt);\r
 548    for(i=0,InApplicable=0; i<n; i++) {\r
 549       if (Root[i]<=0)  {\r
 550          InApplicable=1;\r
 551          Root[i]=-300;  /* adhockery */\r
 552       }\r
 553       else \r
 554          Root[i]=(alpha<=0?log(Root[i]):gammap(Root[i],alpha));\r
 555    }\r
 556    FOR (i,n) FOR (j,n) T1[i*n+j]=U[i*n+j]*Root[j];\r
 557    matby (T1, V, Q, n, n, n);\r
 558    for (i=0,t=0; i<n; i++) t-=pi[i]*Q[i*n+i];\r
 559 \r
 560    if(noisy>=9 && InApplicable) printf("Root(P)<0.  adhockery invoked\n"); \r
 561    if(t<=0) error2("err: DistanceREV");\r
 562 \r
 563    FOR (i,n) Root[i]/=t;\r
 564    FOR (i, n) FOR (j,n)  { Q[i*n+j]/=t; if (i-j) Q[i*n+j]=max2(0,Q[i*n+j]); }\r
 565 \r
 566    return (t);\r
 567 }\r
 568 \r
 569 \r
 570 int PatternLS (FILE *fout, double Ft[], double alpha,double space[],int *cond)\r
 571 {\r
 572 /* space[n*n*2]\r
 573 */\r
 574    int n=com.ncode, i,j,k,h, it;\r
 575    double *Q=Ft,*Qt=Q+n*n,*Qm=Qt+n*n;\r
 576    double *pi,*Root,*U, *V, *T1=space, *branch, t;\r
 577    FILE *fdist=gfopen("Distance", "w");\r
 578    \r
 579    if((pi=(double*)malloc((n*n*3+tree.nbranch)*sizeof(double)))==NULL)\r
 580       error2("PatternLS: oom");\r
 581    Root=pi+n;  U=Root+n; V=U+n*n; branch=V+n*n;\r
 582 \r
 583    *cond=0;\r
 584    for (i=0,zero(Qt,n*n),zero(Qm,n*n); i<com.ns; i++) {\r
 585       for (j=0; j<i; j++) {\r
 586          for (h=0,zero(Q,n*n); h<com.npatt; h++) {\r
 587             Q[(com.z[i][h])*n+com.z[j][h]] += com.fpatt[h]/2;\r
 588             Q[(com.z[j][h])*n+com.z[i][h]] += com.fpatt[h]/2;\r
 589          }\r
 590          FOR (k,n*n) Qt[k]+=Q[k]/(com.ns*(com.ns-1)/2);\r
 591          it=i*(i-1)/2+j;\r
 592          SeqDistance[it]=DistanceREV (Q, n, alpha, Root,U,V, pi, space, &k);\r
 593 \r
 594          if (k==-1) { \r
 595             *cond=-1; printf("\n%d&%d: F(t) modified in DistanceREV",i+1,j+1);\r
 596          }\r
 597 \r
 598          fprintf(fdist,"%9.5f",SeqDistance[it]);\r
 599 /*\r
 600 FOR (k,n) \r
 601 if (Q[k*n+k]>0) { printf ("%d %d %.5f\n", i+1, j+1, Q[k*n+k]); }\r
 602 */\r
 603          FOR (k,n*n) Qm[k]+=Q[k]/(com.ns*(com.ns-1)/2); \r
 604       }\r
 605       FPN(fdist);\r
 606    }\r
 607    fclose (fdist);\r
 608    DistanceREV (Qt, n, alpha, Root, U, V, pi, space, &k);\r
 609    if (k==-1) { puts ("F(t) modified in DistanceREV"); }\r
 610 \r
 611    fprintf (fout, "\n\nQ: from average F over pairwise comparisons");\r
 612    OutQ(fout, n, Qt, pi, Root, U, V, T1);\r
 613    fprintf (fout, "\n\nQ: average of Qs over pairwise comparisons\n");\r
 614    fprintf (fout, "(disregard this if very different from the previous Q)");\r
 615    OutQ (fout, n, Qm, pi, Root, U, V, T1);\r
 616 \r
 617    if (tree.nbranch) {\r
 618       fillxc (branch, 0.1, tree.nbranch);\r
 619       LSDistance (&t, branch, testx);\r
 620       OutTreeB (fout);  FPN (fout);\r
 621       FOR (i,tree.nbranch) fprintf(fout,"%9.5f", branch[i]);\r
 622       PMatBranch (Ft, com.ncode, branch, Root, U, V, space);\r
 623    }\r
 624    free(pi);\r
 625    return (0);\r
 626 }\r
 627 \r
 628 int testx (double x[], int np)\r
 629 {\r
 630    int i;\r
 631    double tb[]={1e-5, 99};\r
 632    FOR(i,np) if(x[i]<tb[0] ||x[i]>tb[1]) return(-1);\r
 633    return(0);\r
 634 }\r
 635 \r
 636 int SetBranch (double x[])\r
 637 {\r
 638    int i, status=0;\r
 639    double small=1e-5;\r
 640 \r
 641    FOR (i,tree.nnode)\r
 642       if (i!=tree.root && (nodes[i].branch=x[nodes[i].ibranch])<-small)\r
 643          status=-1;\r
 644    return (status);\r
 645 }\r