src/fisher.cpp
// pcond2 is a 3-dimensional array of probabilities p(I and k and kp) where "I and k and kp" refers
// to states with leaf state I, and compatible with k and kp.
pcond2 is a 3-dimensional array of probabilities p(I and k and kp) where "I and k and kp" refers
to states with leaf state I, and compatible with k and kp.
pcond2 is a 3 - dimensional array of probabilities p ( I and k and kp ) where " I and k and kp " refers to states with leaf state I and compatible with k and kp .
void fill_pcond2(Tree &T, Model &Mod, StateList &sl, Parameters &Par, Array3 &pcond2) {
  long i, j, l, u, v;
  long k, e, a, b;
  long kp, ep, ap, bp;
  double p;
  long npars = T.nedges*Mod.np + Mod.rdf + 1;
  pcond2.resize(T.nstleaves);
  for (i=0; i < T.nstleaves; i++) {
    pcond2[i].resize(npars);
    for(k=0; k < npars; k++) {
      pcond2[i][k].resize(npars);
      for(l=0; l < npars; l++) {
        pcond2[i][k][l] = 0;
      }
    }
  }
  for (i=0; i < T.nstleaves; i++) {
    for(j=0; j < T.nsthidden; j++) {
      for(e=0; e < T.nedges; e++) {
        for(ep=0; ep < T.nedges; ep++) {
          if (e==ep) continue;
          p = Par.r[rootstate(sl.h[j])];
          for(l=0; l < T.nedges; l++) {
            if (l == e || l == ep) continue;
            edgestate(T, l, sl.h[j], sl.l[i], u, v);
            p = p*Par.tm[l][u][v];
          }
          edgestate(T, e, sl.h[j], sl.l[i], a, b);
          k = erc2param(Mod, e, a, b);
          edgestate(T, ep, sl.h[j], sl.l[i], ap, bp);
          kp = erc2param(Mod, ep, ap, bp);
          pcond2[i][k][kp] = pcond2[i][k][kp] + p;
        }
      }
      for(e=0; e < T.nedges; e++) {
        p = 1;
        for(l=0; l < T.nedges; l++) {
          if (l == e) continue;
          edgestate(T, l, sl.h[j], sl.l[i], u, v);
          p = p*Par.tm[l][u][v];
        }
        edgestate(T, e, sl.h[j], sl.l[i], a, b);
        k = erc2param(Mod, e, a, b);
        kp = root2param(T, Mod, rootstate(sl.h[j]));
        pcond2[i][k][kp] = pcond2[i][k][kp] + p;
        pcond2[i][kp][k] = pcond2[i][kp][k] + p;
      }
    }
  }
}
void fill_pcond2 ( Tree & T , Model & Mod , StateList & sl , Parameters & Par , Array3 & pcond2 ) { long i , j , l , u , v ; long k , e , a , b ; long kp , ep , ap , bp ; double p ; long npars = T . nedges * Mod . np + Mod . rdf + 1 ; pcond2 . resize ( T . nstleaves ) ; for ( i = 0 ; i < T . nstleaves ; i ++ ) { pcond2 [ i ] . resize ( npars ) ; for ( k = 0 ; k < npars ; k ++ ) { pcond2 [ i ] [ k ] . resize ( npars ) ; for ( l = 0 ; l < npars ; l ++ ) { pcond2 [ i ] [ k ] [ l ] = 0 ; } } } for ( i = 0 ; i < T . nstleaves ; i ++ ) { for ( j = 0 ; j < T . nsthidden ; j ++ ) { for ( e = 0 ; e < T . nedges ; e ++ ) { for ( ep = 0 ; ep < T . nedges ; ep ++ ) { if ( e == ep ) continue ; p = Par . r [ rootstate ( sl . h [ j ] ) ] ; for ( l = 0 ; l < T . nedges ; l ++ ) { if ( l == e || l == ep ) continue ; edgestate ( T , l , sl . h [ j ] , sl . l [ i ] , u , v ) ; p = p * Par . tm [ l ] [ u ] [ v ] ; } edgestate ( T , e , sl . h [ j ] , sl . l [ i ] , a , b ) ; k = erc2param ( Mod , e , a , b ) ; edgestate ( T , ep , sl . h [ j ] , sl . l [ i ] , ap , bp ) ; kp = erc2param ( Mod , ep , ap , bp ) ; pcond2 [ i ] [ k ] [ kp ] = pcond2 [ i ] [ k ] [ kp ] + p ; } } for ( e = 0 ; e < T . nedges ; e ++ ) { p = 1 ; for ( l = 0 ; l < T . nedges ; l ++ ) { if ( l == e ) continue ; edgestate ( T , l , sl . h [ j ] , sl . l [ i ] , u , v ) ; p = p * Par . tm [ l ] [ u ] [ v ] ; } edgestate ( T , e , sl . h [ j ] , sl . l [ i ] , a , b ) ; k = erc2param ( Mod , e , a , b ) ; kp = root2param ( T , Mod , rootstate ( sl . h [ j ] ) ) ; pcond2 [ i ] [ k ] [ kp ] = pcond2 [ i ] [ k ] [ kp ] + p ; pcond2 [ i ] [ kp ] [ k ] = pcond2 [ i ] [ kp ] [ k ] + p ; } } } }
pcond2 is a 3-dimensional array of probabilities p(I and k and kp) where "I and k and kp" refers
to states with leaf state I, and compatible with k and kp.
pcond2 is a 3 - dimensional array of probabilities p ( I and k and kp ) where " I and k and kp " refers to states with leaf state I and compatible with k and kp .
// Initializations // Loop over all states // Loop over pairs of edges for the edge-edge part of pcond2 // when e==ep, the corresponding entry must be zero. // Compute cond2 probabilities // Loop a single edge for the edge-root part of pcond2 // Take care of root parameters in pcond1 // Recall: In pcond2 with two root parameters k and kp, // we always have p(I and k and kp) = 0.