pamp.c

一个神经网络工具箱

   double *T1=space, *P;

   FOR (k, tree.nbranch) {
      P=Ptb+k*n*n;
      FOR (i,n) FOR (j,n) T1[i*n+j]=U[i*n+j]*exp(Root[j]*branch[k]);
      matby (T1, V, P, n, n, n);
      FOR (i,n*n) if (P[i]<0) P[i]=0;
/*
      printf ("\nbranch %d, P(%.5f)", k+1, branch[k]);
      matout (F0, P, n, n);
      testTransP (P, n);
*/
   }
   return (0);
}


int PatternMP (FILE *fout, double Ft[])
{
/* Ft[]: input counts for the F(t) matrix for each branch, output P(t) 
*/
   int n=com.ncode, i,j,k;
   double *Q, *pi, *Root, *U, *V, *branch, *space, *T1, t;

   if((Q=(double*)malloc((n*n*6+tree.nbranch)*sizeof(double)))==NULL)
      error2("PathwayMP: oom");
   pi=Q+n*n; Root=pi+n; U=Root+n; V=U+n*n; branch=V+n*n; 
   space=T1=branch+tree.nbranch;

   for (k=0; k<tree.nbranch; k++) {  /* branch lengths */
      xtoy(Ft+k*n*n, Q, n*n);
      branch[k]=nodes[tree.branches[k][1]].branch=
       DistanceREV (Q, n, 0, Root, U, V, pi, space, &j);
   }
   OutaTreeB (fout);  FPN (fout);
   FOR (i, tree.nbranch) fprintf(fout,"%9.5f", branch[i]);
   fprintf (fout,"\ntree length: %9.5f\n", sum(branch,tree.nbranch));

   /* pattern Q from average F(t) */
   fprintf(fout,"\nF(t)");
   xtoy (Ft+tree.nbranch*n*n, Q, n*n);
   matout2 (fout, Q, n, n, 12, 2);
   DistanceREV (Q, n, 0, Root, U, V, pi, space, &j);
   if (noisy>=3&&j==-1) { puts("F(t) modified in DistanceREV"); }

   OutQ (fout, n, Q, pi, Root, U, V, T1);
   if (com.nhomo==0) 
      PMatBranch (Ft, n, branch, Root, U, V, space);
   else {
      for (k=0; k<tree.nbranch; k++) {
         for (i=0; i<n; i++) {
            t=sum(Ft+k*n*n+i*n, n);
            if (t>1e-5) abyx (1/t, Ft+k*n*n+i*n, n);
            else        Ft[k*n*n+i*n+i]=1;
         }
      }
   }
   free(Q);
   return (0); 
}


int PathwayMP1 (FILE *fout, int *maxchange, int NSiteChange[], 
    double Ft[], double space[], int job)
{
/* Hartigan, JA.  1973.  Minimum mutation fits to a given tree. 
   Biometrics, 29:53-65.
   Yang, Z.  1996.  
   job=0: 1st pass: calculates maxchange, NSiteChange[], and Ft[]
   job=1: 2nd pass: reconstructs ancestral character states (->fout)
*/
   char *pch=(com.seqtype==0?BASEs:(com.seqtype==2?AAs:BINs));
   char *zz[NNODE], visit[NS-1],nodeb[NNODE],bestPath[NNODE-NS],Equivoc[NS-1];
   int n=com.ncode, nid=tree.nbranch-com.ns+1, it,i1,i2, i,j,k, h, npath;
   int *Ftt=NULL, nchange, nchange0;
   double sumpr, bestpr, pr, *pnode=NULL, *pnsite;

   if((pnsite=(double*)malloc((com.ns-1)*n*sizeof(double)))==NULL)
      error2("PathwayMP1: oom");

   PATHWay=(char*)malloc(nid*(n+3)*sizeof(char));
   NCharaCur=PATHWay+nid;  ICharaCur=NCharaCur+nid;  CharaCur=ICharaCur+nid;
   if (job==0) {
      zero(Ft,n*n*(tree.nbranch+1));
      if((Ftt=(int*)malloc(n*n*tree.nbranch*sizeof(int)))==NULL) error2("oom");
   }
   else {
      pnode=(double*)malloc((nid*com.npatt+1)*(sizeof(double)+sizeof(char)));
      FOR (j,nid) zz[com.ns+j]=(char*)(pnode+nid*com.npatt)+j*com.npatt;
      FOR (j,com.ns) zz[j]=com.z[j];
      if (pnode==NULL) error2 ("oom");
   }
   for (j=0,visit[i=0]=tree.root-com.ns; j<tree.nbranch; j++) 
      if (tree.branches[j][1]>=com.ns) visit[++i]=tree.branches[j][1]-com.ns;

   for (h=0; h<com.npatt; h++) {
      if (job==1) {
         fprintf (fout, "\n%4d%6.0f  ", h+1, com.fpatt[h]);
         FOR (j, com.ns) fprintf (fout, "%c", pch[com.z[j][h]]);
         fprintf (fout, ":  ");
         FOR (j,nid*n) pnsite[j]=0;
      }
      FOR (j,com.ns) nodeb[j]=com.z[j][h];
      if (job==0) FOR (j,n*n*tree.nbranch) Ftt[j]=0;

      InteriorStatesMP (1, h, &nchange, NCharaCur, CharaCur, space); 
      ICharaCur[j=tree.root-com.ns]=0;  PATHWay[j]=CharaCur[j*n+0];
      FOR (j,nid) Equivoc[j]=(NCharaCur[j]>1);

      if (nchange>*maxchange) *maxchange=nchange;
      if (nchange>NCATCHANGE-1) error2 ("raise NCATCHANGE");
      NSiteChange[nchange]+=(int)com.fpatt[h];

      DownStates (tree.root);
      for (npath=0,sumpr=bestpr=0; ;) {
         for (j=0,k=visit[nid-1]; j<NCharaCur[k]; j++) {
            PATHWay[k]=CharaCur[k*n+j]; npath++;
            FOR (i,nid) nodeb[i+com.ns]=PATHWay[i];
            if (job==1) {
               FOR (i,nid) fprintf(fout,"%c",pch[PATHWay[i]]); fputc(' ',fout);
               pr=com.pi[(int)nodeb[tree.root]];
               for (i=0; i<tree.nbranch; i++) {
                  i1=nodeb[tree.branches[i][0]]; i2=nodeb[tree.branches[i][1]];
                  pr*=Ft[i*n*n+i1*n+i2];
               }
               sumpr+=pr;       FOR (i,nid) pnsite[i*n+nodeb[i+com.ns]]+=pr;
               if (pr>bestpr) { bestpr=pr; FOR(i,nid) bestPath[i]=PATHWay[i];}
            }
            else {
               for (i=0,nchange0=0; i<tree.nbranch; i++) {
                  i1=nodeb[tree.branches[i][0]]; i2=nodeb[tree.branches[i][1]];
                  nchange0+=(i1!=i2);
                  Ftt[i*n*n+i1*n+i2]++;
               }
               if (nchange0!=nchange) 
               { puts("\a\nerr:PathwayMP"); fprintf(fout,".%d. ", nchange0); }
            }
         }
         for (j=nid-2; j>=0; j--) {
            if(Equivoc[k=visit[j]] == 0) continue;
            if (ICharaCur[k]+1<NCharaCur[k]) {
               PATHWay[k] = CharaCur[k*n + (++ICharaCur[k])];
               DownStates (k+com.ns);
               break;
            }
            else { /* if (next equivocal node is not ancestor) update node k */
               for (i=j-1; i>=0; i--) if (Equivoc[(int)visit[i]]) break;
               if (i>=0) { 
                  for (it=k+com.ns,i=visit[i]+com.ns; ; it=nodes[it].father)
                     if (it==tree.root || nodes[it].father==i) break;
                  if (it==tree.root)
                     DownStatesOneNode (k+com.ns, nodes[k+com.ns].father);
               }
            }
         }
         if (j<0) break;
      }      /* for (npath) */
/*
      printf ("\rsite pattern %4d/%4d: %6d%6d", h+1,com.npatt,npath,nchange);
*/      
      if (job==0) 
         FOR (j,n*n*tree.nbranch) Ft[j]+=(double)Ftt[j]/npath*com.fpatt[h];
      else {
         FOR (i,nid) zz[com.ns+i][h]=bestPath[i];
         FOR (i,nid) pnode[i*com.npatt+h]=pnsite[i*n+bestPath[i]]/sumpr;
         fprintf (fout, " |%4d (%d) | ", npath, nchange);
         if (npath>1) {
            FOR (i,nid) fprintf (fout, "%c", pch[bestPath[i]]);
            fprintf (fout, " (%.3f)", bestpr/sumpr);

         }
      }
   }   /* for (h) */
   free(PATHWay); 
   if (job==0) {
      free(Ftt);
      FOR (i,tree.nbranch) FOR (j,n*n) Ft[tree.nbranch*n*n+j]+=Ft[i*n*n+j];
   }
   else {
      fprintf (fout,"\n\nApprox. relative accuracy at each node\n");
      FOR (h, com.npatt) {
         fprintf (fout,"\n%4d%6.0f  ", h+1, com.fpatt[h]);
         FOR (j, com.ns) fprintf (fout, "%c", pch[com.z[j][h]]);
         fprintf (fout, ":  ");
         FOR (i,nid) if (pnode[i*com.npatt+h]<.99999) break;
         if (i<nid)  FOR (j, nid) 
            fprintf(fout,"%c (%5.3f) ", pch[zz[j][h]],pnode[j*com.npatt+h]);
      }
      Site2Pattern (fout);
      fprintf (fout,"\n\nlist of extant and reconstructed sequences\n\n");
      for(j=0;j<tree.nnode;j++,FPN(fout)) {
         if(j<com.ns) fprintf(fout,"%-20s", com.spname[j]);
         else         fprintf(fout,"node #%-14d", j+1);
         print1seq (fout, zz[j], com.ls, 1, com.pose);
      }
      free (pnode);
   }
   free(pnsite);
   return (0);
}

double DistanceREV (double Ft[], int n,double alpha,double Root[],double U[],
   double V[], double pi[], double space[], int *cond)
{
/* input:  Ft, n, alpha
   output: Q(in Ft), t, Root, U, V, and cond
   space[n*n*2]
*/
   int i, j;
   double *Q=Ft, *T1=space, *T2=space+n*n, t, small=1e-6;
   
   *cond=0;
   for (i=0,t=0; i<n; i++) FOR (j,n) if (i-j) t+=Q[i*n+j];
   if (t<small)  { *cond=1; zero(Q,n*n); return (0); }

   for (i=0;i<n;i++) for (j=0;j<i;j++) Q[i*n+j]=Q[j*n+i]=(Q[i*n+j]+Q[j*n+i])/2;
   abyx (1./sum(Q,n*n), Q, n*n);

   FOR (i,n) {
      pi[i]=sum(Q+i*n, n);  
/*
      if (Q[i*n+i]<=small || Q[i*n+i]<pi[i]/4)
*/
      if (Q[i*n+i]<=small)
         {  Q[i*n+i]=1-pi[i]+Q[i*n+i]; *cond=-1; }

      else  abyx(1/pi[i], Q+i*n, n); 
   }
   if (eigen (1, Q, n, Root, T1, U, V, T2)) error2 ("eigen jgl");
   xtoy (U, V, n*n);
   matinv (V, n, n, T1);

   FOR (i,n) {
      if (Root[i]<=0)  {
         printf ("  Root %d:%10.4f", i+1, Root[i]); 
      }
      Root[i]=(alpha<=0?log(Root[i]):gammap(Root[i],alpha));
   }
   FOR (i,n) FOR (j,n) T1[i*n+j]=U[i*n+j]*Root[j];
   matby (T1, V, Q, n, n, n);
   for (i=0,t=0; i<n; i++) t-=pi[i]*Q[i*n+i];
   if (t<=0) puts ("err: DistanceREV");

   FOR (i,n) Root[i]/=t;
   FOR (i, n) FOR (j,n)  { Q[i*n+j]/=t; if (i-j) Q[i*n+j]=max2(0,Q[i*n+j]); }

   return (t);
}


int PatternLS (FILE *fout, double Ft[], double alpha,double space[],int *cond)
{
/* space[n*n*2]
*/
   int n=com.ncode, i,j,k,h, it;
   double *Q=Ft,*Qt=Q+n*n,*Qm=Qt+n*n;
   double *pi,*Root,*U, *V, *T1=space, *branch, t;
   FILE *fdist=gfopen("Distance", "w");
   
   if((pi=(double*)malloc((n*n*3+tree.nbranch)*sizeof(double)))==NULL)
      error2("PatternLS: oom");
   Root=pi+n;  U=Root+n; V=U+n*n; branch=V+n*n;

   *cond=0;
   for (i=0,zero(Qt,n*n),zero(Qm,n*n); i<com.ns; i++) {
      for (j=0; j<i; j++) {
         for (h=0,zero(Q,n*n); h<com.npatt; h++) {
	    Q[(com.z[i][h])*n+com.z[j][h]] += com.fpatt[h]/2;
            Q[(com.z[j][h])*n+com.z[i][h]] += com.fpatt[h]/2;
	 }
         FOR (k,n*n) Qt[k]+=Q[k]/(com.ns*(com.ns-1)/2);
         it=i*(i-1)/2+j;
	 SeqDistance[it]=DistanceREV (Q, n, alpha, Root,U,V, pi, space, &k);

         if (k==-1) { 
            *cond=-1; printf("\n%d&%d: F(t) modified in DistanceREV",i+1,j+1);
         }

	 fprintf(fdist,"%9.5f",SeqDistance[it]);
/*
FOR (k,n) 
if (Q[k*n+k]>0) { printf ("%d %d %.5f\n", i+1, j+1, Q[k*n+k]); }
*/
         FOR (k,n*n) Qm[k]+=Q[k]/(com.ns*(com.ns-1)/2); 
      }
      FPN(fdist);
   }
   fclose (fdist);
   DistanceREV (Qt, n, alpha, Root, U, V, pi, space, &k);
   if (k==-1) { puts ("F(t) modified in DistanceREV"); }

   fprintf (fout, "\n\nQ: from average F over pairwise comparisons");
   OutQ(fout, n, Qt, pi, Root, U, V, T1);
   fprintf (fout, "\n\nQ: average of Qs over pairwise comparisons\n");
   fprintf (fout, "(disregard this if very different from the previous Q)");
   OutQ (fout, n, Qm, pi, Root, U, V, T1);

   if (tree.nbranch) {
      fillxc (branch, 0.1, tree.nbranch);
      LSDistance (&t, branch, testx);
      OutaTreeB (fout);  FPN (fout);
      FOR (i,tree.nbranch) fprintf(fout,"%9.5f", branch[i]);
      PMatBranch (Ft, com.ncode, branch, Root, U, V, space);
   }
   free(pi);
   return (0);
}

int testx (double x[], int np)
{
   int i;
   double tb[]={1e-5, 99};
   FOR(i,np) if(x[i]<tb[0] ||x[i]>tb[1]) return(-1);
   return(0);
}