pmdsubs.c

#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include <nicklib.h>  

#include "nicksam.h"
#include "popsubs.h" 
#include "pmdsubs.h"

extern double pontus_p ;   //   .3
extern double pontus_c ;   //  .01
extern double pontus_pi;   // .001

extern int verbose ;
extern long numpoly ; 

extern int udgtype ;
extern char *regstring ;
extern int  reglen ;
extern int pmdoff ;


extern int pulldownmode ; // might be pullad

static int pmdhlo = -20, pmdhhi = 20, pmdhlen=0, pmdnpops=0 ;
static int **pmdhist  = NULL;
static double **pmdmoments ;

static double *****zdd = NULL, *****zee = NULL, *zdamvals = NULL ;
static int emnumiter = 1000 ;

int setudgtype(char *udgmode) 
{
 int t ;
 if (udgmode==NULL) return -1 ;

 t = strcmp(udgmode, "+") ;   ; if (t==0) return 2 ;
 t = strcmp(udgmode, "plus")  ; if (t==0) return 2 ;
 t = strcmp(udgmode, "0")     ; if (t==0) return 0 ;
 t = strcmp(udgmode, "minus") ; if (t==0) return 0 ;
 t = strcmp(udgmode, "1/2")   ; if (t==0) return 1 ;
 t = strcmp(udgmode, "half")  ; if (t==0) return 1 ;
 t = strcmp(udgmode, "single")  ; if (t==0) return 3 ;
 t = strcmp(udgmode, "nick22")  ; if (t==0) return 10 ;
 t = strcmp(udgmode, "nick25")  ; if (t==0) return 20 ;

 printf("*** warning: udgmode not recognized: pmd not used\n") ;
 return -1 ;

}
int pontus_dis_single(int pp, int k, int len, int strand, int *dx, int *dy) 
// udg processing goes here 
// dis < 0 => base not useful for scoring
{
   int t, kk, x, tmatch ; 
   char rbase, rr ;
   int tt[3], dd[3], kmin, kmax ; 

   if (pp<2) return -3  ;
   if (pp>=(reglen-1)) return -3  ;
   rbase = toupper(regstring[pp-1]) ;

   t = base2num(rbase) ;
   if (t<0) return -2 ;
   if (t==0) return -1 ;
   if (t==3) return -1 ;

   if ((t==1) && (strand==1)) return -1 ;
   if ((t==2) && (strand==0)) return -1 ;

   kk = len-(k+1) ;
   kmax = MAX(k, kk) ;
   kmin = MIN(k, kk) ;

   if (strand==0) { 
    tt[0] = 0 ; 
    tt[1] = len-2 ;
    tt[2] = len-1 ;
    tmatch = 1 ;  // reg: C
   }
   if (strand==1) { 
    tt[0] = len-1 ; 
    tt[1] = 1 ;
    tt[2] = 0 ;
    tmatch = 2 ; // reg : G 
   } 
   dd[0] = 1 ;
   dd[1] = 2 ;
   dd[2] = 1 ;

   if (t==tmatch) {
    x = findfirst(tt, 3, k) ;
    if (x>=0) { 
     *dx = dd[x] ; 
     *dy = kmax+1 ;
     return *dx ; 
    }
   }

/* CpG */

   if (t==1) {   // pp -> C
    rr = toupper(regstring[pp]) ;      
    if (rr != 'G') return -1 ;
     *dx = kmin+1 ;
     *dy = kmax+1 ;
     return *dx ; 
   }

   if (t==2) {   // pp -> G
    rr = toupper(regstring[pp-2]) ;      
    if (rr != 'C') return -1 ;
     *dx = kmin+1 ;
     *dy = kmax+1 ;
     return *dx ; 
   }

   fatalx("(pontus_dis_single)") ;
}
int pontus_match_single(char rbase, char cseq, int strand)   
{

  if (rbase==cseq) return YES ;
  if ((rbase=='C')&&(cseq=='T')&&(strand == 0)) return NO ;
  if ((rbase=='G')&&(cseq=='A')&&(strand == 1)) return NO ;

  return -1 ; 



}
int pontus_dis(int pp, int k, int len, int udgtype) 
// udg processing goes here 
// dis < 0 => base not useful for scoring
{
   int t, kk ; 
   char rbase, rr ;

   if (pp<2) return -3  ;
   if (pp>=(reglen-1)) return -3  ;
   rbase = toupper(regstring[pp-1]) ;

   t = base2num(rbase) ;
   if (t<0) return -2 ;
   if (t==0) return -1 ;
   if (t==3) return -1 ;

   kk = len-(k+1) ;

   if (t==1) {   // pp -> C
    if (udgtype==0) return k+1  ;  // UDG-
    if ((udgtype==1) && (k==0)) return k+1 ;  ///UDG 1/2
    rr = toupper(regstring[pp]) ;      
    if (rr != 'G') return -1 ;
    return k + 1 ;
   }

   if (t==2) {   // pp -> G
    if (udgtype==0) return kk+1  ;  // UDG-
    if ((udgtype==1) && (kk==0)) return kk+1 ;
    rr = toupper(regstring[pp-2]) ;      
    if (rr != 'C') return -1 ;
    return kk + 1 ;
   }

   fatalx("(pontus_dis)") ;
}
int pontus_match(char rbase, char cseq)   
{

  if (rbase==cseq) return YES ;
  if ((rbase=='C')&&(cseq=='T')) return NO ;
  if ((rbase=='G')&&(cseq=='A')) return NO ;

  return -1 ; 



}
double pmdeps(int qual) 
{
  double y ; 

  if (qual<=0) return 1.0/3.0 ;
  y = (double) qual / 10.0 ;
  return pow(0.1, y)/3.0 ;

}
double pmatch(double pi, double dz, double eps) 
{
 double y ; 
 y = (1.0-pi)*(1.0-dz)*(1.0-eps) ;
 y += (1.0-pi)*dz*eps ;
 y += pi*eps*(1.0-dz) ;

 return y ;
}

double pmatch_single(double pi, double dzx, double dzy, double eps) 
// equation S1 of supplementary of Skoglund et al. 
{
 double y ; 
 y = (1.0-pi)*(1.0-dzx)*(1.0-dzy)*(1.0-eps) ;
 y += (1.0-pi)*dzx*(1.0-dzy)*eps ;
 y += (1.0-pi)*dzy*(1.0-dzx)*eps ;
 y += pi*eps*(1.0-dzx)*(1.0-dzy) ;

 return y ;

}

double pmd_single(READ *readpt, int pmdoff) 
// single strand procssing
{

  int qual, k, len, pos, pp, dis, mtch, *rdistx, *rdisty, *match ;
  char rbase, cseq, cref ;
  double yscore = 0.0, yeps, y1, y2, yp, ypi, pmodel, pnull, yc, dzx, dzy ;
  int strand, dx, dy ;

  len = readpt -> readlen ;
  pos = readpt -> pos ;
  strand = readpt -> strand ;

  if (len<=0) return 0.0 ;

  ZALLOC(rdistx, len, int) ;
  ZALLOC(rdisty, len, int) ;
  ZALLOC(match, len, int) ;


  ivclear(rdistx, -1, len) ; // -1 => ignore
  ivclear(rdisty, -1, len) ; // -1 => ignore
  ivclear(match, -1, len) ; // -1 => ignore

  for (k=0; k<len; ++k) { 
   cref = readpt -> refseq[k] ;
   if (pulldownmode && (isupper(cref))) continue ;  // dangerous bend; specific to pulldown
   pp = pos + k ;
   if (pp<0) continue ;
   if (pp>=reglen) continue ;
   dis = pontus_dis_single(pp, k, len, strand, &dx, &dy) ;
   if (dis<0) continue ;  

   rbase = toupper(cref) ;
   cseq = toupper(readpt -> seq[k]) ;  
   mtch = pontus_match_single(rbase, cseq, strand) ;
   if (mtch < 0) continue ; 
   rdistx[k] = dx ;
   rdisty[k] = dy ;
   match[k] = mtch ;
  }

  for (k=0; k<len; ++k) { 
   dx = rdistx[k] ; 
   if (dx<0) continue ;
   dy = rdisty[k] ; 
   yp = pontus_p ; 
   yc = pontus_c ; 
   ypi = pontus_pi ; 
   dzx = pow(1.0-yp, dx)*yp + yc ;
   dzy = pow(1.0-yp, dy)*yp + yc ;
   qual = readpt -> iqual[k] - pmdoff ;
   yeps = pmdeps(qual) ;
   pmodel = pmatch_single(ypi, dzx, dzy, yeps) ;
   pnull =  pmatch_single(ypi, 0.0, 0.0, yeps) ; 
   if (match[k]==NO) { 
    pmodel = 1.0-pmodel ;
    pnull  = 1.0-pnull ;
   }
   yscore += log(pmodel/pnull) ;
  }

  free(rdistx) ;
  free(rdisty) ;
  free(match) ;

  return yscore ;

}
double pmd_half(READ *readpt, char *polys, char *cpgs, int *****htab) 
{

  int qual, k, len, pos, pp, dis, mtch, *rdist, *match ;
  char rbase, cseq, cref ;
  double yscore = 0.0, yeps, y1, y2, yp, ypi, pmodel, pnull, yc, dz ;
  int pmdoff = 0 ; 
  static long ncall = 0 ; 

  char *poly, *cpg ; 


  ++ncall ; 
  len = readpt -> readlen ;
  pos = readpt -> pos ;

  poly = polys + pos - 1 ;
  cpg = cpgs + pos - 1 ;

  if (len<=0) return 0.0 ;

  ZALLOC(rdist, len, int) ;
  ZALLOC(match, len, int) ;
  ivclear(rdist, -1, len) ; // -1 => ignore
  ivclear(match, -1, len) ; // -1 => ignore

  for (k=0; k<len; ++k) { 
   cref = readpt -> refseq[k] ;
// if (pulldownmode && (isupper(cref))) continue ;  // dangerous bend; specific to pulldown
   pp = pos + k ;
   if (pp<0) continue ;
   if (pp>=reglen) continue ;
   dis = pontus_dis(pp, k, len, udgtype) ;
   if (dis<0) continue ;  

   rbase = toupper(cref) ;
   cseq = toupper(readpt -> seq[k]) ;  
   mtch = pontus_match(rbase, cseq) ;
   if (mtch < 0) continue ; 
   rdist[k] = dis ;
   match[k] = mtch ;
  }

   yp = pontus_p ; 
   yc = pontus_c ; 
   ypi = pontus_pi ; 
   dz = pow(1.0-yp, dis)*yp + yc ;
  for (k=0; k<len; ++k) { 
   dis = rdist[k] ; 
   if (dis<0) continue ;
   if (poly[k] == '1') continue ; 
   qual = readpt -> iqual[k] - pmdoff ;
   yeps = pmdeps(qual) ;
   pmodel = pmatch(ypi, dz, yeps) ;
   pnull =  pmatch(ypi, 0.0, yeps) ; 
   if (match[k]==NO) { 
    pmodel = 1.0-pmodel ;
    pnull  = 1.0-pnull ;
   }
   yscore += log(pmodel/pnull) ;
  }

  free(rdist) ;
  free(match) ;

  return yscore ;


}

/**
HISTS ***bighistlist ;  
HISTS **bighist ;  
int  *bghlen, bighistlen ;
HISTINFO **histinfo ;  
*/
int skat(char ca, char cb) 
{

  if ((ca == 'C') && (cb == 'C')) return 0 ;
  if ((ca == 'C') && (cb == 'T')) return 0 ;
  if ((ca == 'G') && (cb == 'G')) return 1 ;
  if ((ca == 'G') && (cb == 'A')) return 1 ;


 return -1 ; 



}

double deaminf(HISTS *hpt, double ap, double bp, int kcat) 
{
  int r, c; 
  double y1, y2 ; 
  if (kcat < 0) return 0.5 ;

  if (kcat == 0) { 
    r = 1 ; c = 3 ; 
  }

  if (kcat == 1) { 
    r = 2 ; c = 0 ; 
  }
  if (kcat > 1) fatalx("bad bug\n") ;

  y1 = hpt -> hist[4*r+c] + ap ; 
  y2 = hpt -> hist[4*r+r] + bp ;
  if (verbose) printf("zzdeaminf %9.3f %9.3f %d\n", y1, y2, intsum(hpt -> hist, 16)) ;
  return  y1/(y1+y2) ;

}



double pmd_nick25(READ *readpt, HISTINFO **histinfo, HISTS ***bighistlist, char *polys, char *cpgs, int *****htab) 

{

  int t, qual, i, j, k, x, len,  pos, pp, a, abspp ;
  int t0, t1 ;
  char rbase, cseq, cref, cc ;
  double yscore = 0.0, y1, y2 ;
  int rcat, kcat ; 
  HISTS **bighist, *hpt, *hpt0 ;  
  HISTINFO *hinfo ; 
  double **eef, **ddf ; 

  int maxmq, maxbq, mq, bq, strand, readedge ; 
  static long ncall = 0 ; 
  static double aprior, bprior ; 

  char *poly = NULL, *cpg = NULL ; 
  int r, c, u ;
  double dscore[2], damval, yy[2], z0, z1 ; 

  len = readpt -> readlen ;
  pos = readpt -> pos ;

  if (polys != NULL)  poly = polys + pos - 1 ;
  if (cpgs != NULL) cpg = cpgs + pos - 1 ; // not used at present


  if (len<=0) return 0.0 ;
  rcat = readpt -> rgcat ;
  if (rcat<0) return yscore ;

  ++ncall ;

  bighist = bighistlist[rcat] ;
  hinfo = histinfo[rcat] ; 
  

  mq = MIN(hinfo -> maxmq, readpt -> mapq) ;
  if (mq<0) fatalx("bad hinfo  maxmq unset!\n") ;

  strand = readpt -> strand ; 
  readedge = hinfo -> nump ;   
  
  cc = '?' ; 
  damval = zdamvals[rcat] ; 
  dscore[0] = log(1.0-damval) ; 
  dscore[1] = log(damval) ; 
  for (k=0;  k<len; ++k) {
   cref = readpt -> refseq[k] ;
   cseq = readpt -> seq[k] ; 

   if (poly != NULL) cc = poly[k] ;
   if (verbose) printf("zzpnick25 %d %c %c %c %9.3f\n", k, cc, cref, cseq, yscore) ;
   if (cc == '1') { 
    ++numpoly ;  
    continue ;
   }

   cref = toupper(cref) ; 
   cseq = toupper(cseq) ; 
   r = base2num(cref) ; 
   c = base2num(cseq) ; 
   if (r<0) continue ;
   if (c<0) continue ;

   bq = readpt -> iqual[k] ; 
   bq = MIN(hinfo -> maxbq, bq) ; 

   j = readpt -> readlen - (k+1)  ;
   if (k<j) {
    pp = k+1 ;
   }
   else {
    pp = -(j+1) ;
  }
  abspp = abs(pp) ;
  if (abspp>readedge) pp = 99 ; 
  
  a = ptoa(pp, readedge) ; 
  eef = zee[rcat][mq][bq] ; 
  ddf = zdd[rcat][strand][a] ; 
/** 
 from nick25_damage.tex 
\begin{eqnarray}
 z(r,v; m,b,p,s) & = & d z_0(r,v; m,b,p,s) +  (1-d) z_1(r,; m,b,p,s) \ where \\
 z_0(r,v; m,b,p,s) & = & E(r,v, m,b) \\
 z_1(r,v; m,b,p,s) & = & \sum_u E(u, v,  m,b) D(r,  u; p, s) \\
\end{eqnarray}
*/
  z0 = eef[r][c] ; 
  z1 = 0 ; 
  for (u=0; u<=3; ++u)  {  
   z1 += eef[u][c] * ddf[r][u] ; 
  }
  dscore[0] += log(z0) ; 
  dscore[1] += log(z1) ; 

   if (verbose) { 
    printf("zzpmd25 %d %d %d %d ", mq, bq, strand, pp) ; 
    printf(" %c %c ", cref, cseq) ; 
    printf(" %9.3f", z0) ;
    printf(" %9.3f", z1) ;
    printf(" %9.3f", log(z1/z0)) ;
    printnl() ;
   }
 }
 yscore = dscore[1] - dscore[0] ;
 
 return yscore ; 

}


double pmd_nick(READ *readpt, HISTINFO **histinfo, HISTS ***bighistlist, char *polys, char *cpgs, int *****htab) 

{

  int t, qual, i, j, k, x, len,  pos, pp ;
  int t0, t1 ;
  char rbase, cseq, cref, cc ;
  double yscore = 0.0, y1, y2 ;
  int rcat, kcat ; 
  HISTS **bighist, *hpt, *hpt0 ;  
  HISTINFO *hinfo ; 

  int maxmq, maxbq, mq, bq, strand, readedge ; 
  static long ncall = 0 ; 
  static double aprior, bprior ; 

  char *poly = NULL, *cpg = NULL ; 

  len = readpt -> readlen ;
  pos = readpt -> pos ;

  if (polys != NULL)  poly = polys + pos - 1 ;
  if (cpgs != NULL) cpg = cpgs + pos - 1 ; // not used at present


  if (len<=0) return 0.0 ;
  rcat = readpt -> rgcat ;
  if (rcat<0) return yscore ;

  ++ncall ;

  if (ncall==1) { 
   bpars(&aprior, &bprior, 0.01, 0.01*0.01) ; 
   printf("beta prior for deamination prob: %9.3f %9.3f\n", aprior, bprior) ;
// printf("poly: %x\n", poly) ;
  }

  bighist = bighistlist[rcat] ;
  hinfo = histinfo[rcat] ; 
  

  mq = MIN(hinfo -> maxmq, readpt -> mapq) ;
  if (mq<0) fatalx("bad hinfo  maxmq unset!\n") ;

  strand = readpt -> strand ; 
  readedge = hinfo -> nump ;   
  

  
  cc = '?' ; 
  for (k=0;  k<len; ++k) {
   cref = readpt -> refseq[k] ;
   cseq = readpt -> seq[k] ; 

   if (poly != NULL) cc = poly[k] ;
   if (verbose) printf("zzpnick %d %c %c %c %9.3f\n", k, cc, cref, cseq, yscore) ;
   if (cc == '1') { 
    ++numpoly ;  
    continue ;
   }

   cref = toupper(cref) ; 
   cseq = toupper(cseq) ; 

   kcat = skat(cref, cseq) ; 
   if (verbose) printf("zzsk %d %c %c\n", kcat, cref, cseq) ;

  if (kcat < 0) continue ; 
   
   bq = readpt -> iqual[k] ; 
   bq = MIN(hinfo -> maxbq, bq) ; 

   j = readpt -> readlen - (k+1)  ;
   if (k<j) {
    pp = k+1 ;
   }
   else {
    pp = -(j+1) ;
  }
  
   x = ptoa(pp, readedge) ; 
   if (x<0) continue ; // interior of read
   t = htab[mq][bq][strand][x][rcat] ;
   if (t<0) continue ; 
   hpt = bighist[t] ;

   x = ptoa(99, readedge) ; 
   t0 = htab[mq][bq][strand][x][rcat] ;
   if (t0<0) continue ; 
   hpt0 = bighist[t0] ;

   if (hpt == hpt0) continue ;


  y1 = deaminf(hpt, aprior, bprior, kcat) ;
  y2 = deaminf(hpt0, aprior, bprior, kcat) ;

   if (cseq == cref) { 
    y1 = 1.0 - y1 ;
    y2 = 1.0 - y2 ;
   }

   if (y1==0.0) fatalx("pmdbug3\n") ; 
   if (y2==0.0) fatalx("pmdbug4\n") ; 

   yscore += log(y1/y2) ; 
   if (verbose) { 
    printf("zzpmd %d %d %d %d ", mq, bq, strand, pp) ; 
    printf(" %c %c ", cref, cseq) ; 
    printf(" %9.3f", y1) ;
    printf(" %9.3f", y2) ;
    printf(" %9.3f", log(y1/y2)) ;
    printnl() ;
   }
 }

 
 return yscore ; 

}

double pmd_pontus(READ *readpt, int udgtype, int pmdoff) 
{

  int qual, k, len, pos, pp, dis, mtch, *rdist, *match ;
  char rbase, cseq, cref ;
  double yscore = 0.0, yeps, y1, y2, yp, ypi, pmodel, pnull, yc, dz ;

  if (verbose) { 
   printf("pmd: %d\n", udgtype) ;
   printread(readpt) ;
  }
  if (udgtype<0) return 0.0 ;
  if (udgtype==3) return pmd_single(readpt, pmdoff) ;

  if (udgtype > 3) fatalx("bad udgtype\n") ;

  len = readpt -> readlen ;
  pos = readpt -> pos ;

  if (len<=0) return 0.0 ;

  ZALLOC(rdist, len, int) ;
  ZALLOC(match, len, int) ;
  ivclear(rdist, -1, len) ; // -1 => ignore
  ivclear(match, -1, len) ; // -1 => ignore

  for (k=0; k<len; ++k) { 
   cref = readpt -> refseq[k] ;
   if (pulldownmode && (isupper(cref))) continue ;  // dangerous bend; specific to pulldown
   pp = pos + k ;
   if (pp<0) continue ;
   if (pp>=reglen) continue ;
   dis = pontus_dis(pp, k, len, udgtype) ;
   if (dis<0) continue ;  

   rbase = toupper(cref) ;
   cseq = toupper(readpt -> seq[k]) ;  
   mtch = pontus_match(rbase, cseq) ;
   if (mtch < 0) continue ; 
   rdist[k] = dis ;
   match[k] = mtch ;
  }

  for (k=0; k<len; ++k) { 
   dis = rdist[k] ; 
   if (dis<0) continue ;
   yp = pontus_p ; 
   yc = pontus_c ; 
   ypi = pontus_pi ; 
   dz = pow(1.0-yp, dis)*yp + yc ;
   qual = readpt -> iqual[k] - pmdoff ;
   yeps = pmdeps(qual) ;
   pmodel = pmatch(ypi, dz, yeps) ;
   pnull =  pmatch(ypi, 0.0, yeps) ; 
   if (match[k]==NO) { 
    pmodel = 1.0-pmodel ;
    pnull  = 1.0-pnull ;
   }
   yscore += log(pmodel/pnull) ;
  }

  free(rdist) ;
  free(match) ;

  return yscore ;

}

void clearpmdhist(int npops) 
{
 pmdhlen  = pmdhhi - pmdhlo + 1 ;
 pmdnpops = npops ; 
 pmdhist = initarray_2Dint (npops, pmdhlen, 0) ;
 pmdmoments = initarray_2Ddouble (npops, 3, 0) ;
}

void destroypmdhist() 
{
 if (pmdhlen==0) return ;
 free2Dint(&pmdhist, pmdnpops) ;
 free2D(&pmdmoments, pmdnpops) ;
 pmdnpops = pmdhlen = 0 ;
}

void addpmdhist(int ipop, double y) 
{

 int k, x ; 

 if (pmdhlen==0) return ;
 k = nnint(y) ; 
 k = MAX(k, pmdhlo) ;
 k = MIN(k, pmdhhi) ;

 ++pmdhist[ipop][k-pmdhlo] ;
 pmdmoments[ipop][0] += 1 ;
 pmdmoments[ipop][1] += y ;
 pmdmoments[ipop][2] += y * y  ;


}
void printpmdhist()  
{
 double *ww, *w1, y, y0, y1, y2  ; 
 int t, ipop, k, *hh ;
 double ymean, ysig, yvar, ystderr ;
 if (pmdhlen==0) return ;

 ZALLOC(ww, pmdhlen, double) ;
 ZALLOC(w1, pmdhlen, double) ;
 for (ipop=0; ipop < pmdnpops; ++ipop) {
  hh = pmdhist[ipop] ;
  t = intsum(hh, pmdhlen) ;
  if (t==0) {  
   printf("pmdhist: %5d no data\n", ipop) ;
   continue ;
  }
  floatit(ww, hh, pmdhlen) ;
  y = asum(ww, pmdhlen)  ;              
  vst(ww, ww, 1.0/y, pmdhlen) ;
  for (k=0; k<pmdhlen; ++k) { 
   printf("pmdhist: %5d %3d ", ipop, k+pmdhlo) ;
   printf("%12d %12.6f\n", hh[k], ww[k]) ; 
   w1[k] = k + pmdhlo ;
  }

  y0 = pmdmoments[ipop][0] ; 
  y1 = pmdmoments[ipop][1] ; 
  y2 = pmdmoments[ipop][2] ; 

  y0 += 1.0e-10 ;

  ymean = y1/y0 ;
  y2 -= y0*ymean*ymean ; 
  yvar = y2/y0 ;

  ysig = sqrt(yvar + 1.0e-10) ;
  ystderr = ysig/sqrt(y0) ;
  printf("pmdstats: %5d ", ipop) ;
  
  printf("mean, sig, stderr, z: %12.6f %12.6f %12.6f %12.6f\n", 
   ymean, ysig, ystderr, ymean/ystderr) ;
 }

 free(ww) ;
 free(w1) ;
 fflush(stdout) ;
}

int emsetup(double *pdamval, double ******pee, double ******pdd, HISTINFO **histinfo, int numrgcats, int emiter) 
{
  HISTINFO *hipt ; 
  double *****ee, *****dd, **aa  ; 
  int k, i, j, x, s ; 
  int r, c ;
  int maxmq, maxbq, nump, maxx ; 

  emnumiter = emiter ;
  maxmq = maxbq = nump = 0 ; 
  for (k=0; k<numrgcats; ++k)  {  
    hipt = histinfo[k] ;  
    maxmq = MAX(maxmq, hipt -> maxmq) ; 
    maxbq = MAX(maxbq, hipt -> maxbq) ; 
    nump = MAX(nump, hipt -> nump) ; 
  } 
    printf("zzemsetup  %d %d %d\n",  maxmq, maxbq, nump) ; fflush(stdout) ; 
    maxx  = 2*nump ;  
    ee = initarray_5Ddouble(numrgcats, maxmq+1, maxbq+1, 4, 4, 1) ; 
    *pee = ee ; 

   for (k=0; k < numrgcats; ++k) {
    for (i=0; i <= maxmq; ++i)  {
     for (j=0; j <= maxbq; ++j) {
      aa = ee[k][i][j]  ;
      for (r=0; r < 4; ++r) {
       aa[r][r] = 99 ; 
       bal1(aa[r], 4) ; 
      }
     } 
    }
   }
   dd = initarray_5Ddouble(numrgcats, 2, maxx+1, 4, 4, 1) ; 
   *pdd = dd ;  
   for (k=0; k < numrgcats; ++k) {
    for (s=0; s <= 1; ++s)  {
     for (j=0; j <= maxx; ++j) {
      aa = dd[k][s][j]  ;
      for (r=0; r < 4; ++r) {
       aa[r][r] = 99 ; 
       bal1(aa[r], 4) ; 
      }
     }
    }
   }
  vclear(pdamval, 0.01, numrgcats) ;

  return 1  ; 

}

double dogam(double *GX, double **GE, double **GD, 
double dam, double **ee, double **dd, double *ydata)

{

 double g0, g1, yy, y, yl ; 
 double ylike ;
 int r, c, f ; 
 double fudge = .0001 ;

 if (GX != NULL) {
  vzero(GX, 2) ; 
  clear2D(&GE, 4, 4, 0) ;
  clear2D(&GD, 4, 4, 0) ;
 }

 ylike = 0 ;


 for (r=0; r<=3; ++r) {  
  for (c=0; c<=3; ++c) {  
   yy = ydata[4*r+c] ;
   g0 = ee[r][c]*(1-dam) ;
   g1 = 0 ; 
   for (f=0; f<=3; ++f) {  
    g1 += dd[r][f]*ee[f][c] ;
   }
   g1 *= dam ; 
   yl = g0+g1 ;
   ylike += yy*log(yl) ;
   if (GX != NULL) { 
    GX[0] += yy*g0/yl ;
    GX[1] += yy*g1/yl ;
    GE[r][c] += yy*g0/yl ;
    for (f=0; f<=3; ++f) {  
     y = dd[r][f]*ee[f][c]*dam/yl ;
      GD[r][f] += yy*y ;
      GE[f][c] += yy*y ;
    }
   }

 }} 

 return ylike; 

}

double emlike(double *pdamval, double *****eex, double *****ddx, 
 HISTINFO **histinfo, HISTS ***bighistlist, int *****htab, int numrgcats, double ybase, int rcat)  
{
  HISTINFO *hipt, *hinfo ; 
  HISTS **bighist, *hpt, *hpt0 ;  

  double damval, **eef, **ddf ; 
  int k, i, j, x, s, t  ; 
  int r, c ;
  int maxmq, maxbq, nump, maxx ; 
  int mq, bq, strand ;
  double ylike, yy ; 
  double  *Y ;
  double *va, *vb, pdam ;


  ZALLOC(Y, 4*4, double) ;
  
  
  bighist = bighistlist[rcat] ;
  hipt = histinfo[rcat] ; 
  
   maxmq = hipt -> maxmq ; 
   maxbq = hipt -> maxbq ; 
   nump  = hipt -> nump ; 
   maxx  = 2*nump ;  

    ylike = 0 ; 
    for (mq=0; mq<=maxmq; ++mq) { 
     for (bq=0; bq<=maxbq; ++bq) { 
      for (strand=0; strand <2 ; ++strand) {
        for (x=0; x<=maxx; ++x) { 
         eef = eex[rcat][mq][bq] ;
         ddf = ddx[rcat][strand][x] ;
         damval = *pdamval ;

       t = htab[mq][bq][strand][x][rcat] ;
       if (t<0) continue ;
       hpt = bighist[t] ;
       floatit(Y, hpt -> hist, 16) ; 
       ylike += dogam(NULL, NULL, NULL, damval, eef, ddf, Y) ;

     }}}} 

    ylike -= ybase ;  
   

 free(Y) ; 

 return ylike ;
}

double emdamage(double *pdamval, double *****eex, double *****ddx, 
 HISTINFO **histinfo, HISTS ***bighistlist, int *****htab, int numrgcats)  
{
  HISTINFO *hipt, *hinfo ; 
  HISTS **bighist, *hpt, *hpt0 ;  

  double ****gee, ****gdd, gdamval[2] ;
  double damval, **eef, **ddf ; 
  int k, i, j, x, s, t, rcat, scat  ; 
  int r, c, pp ;
  int maxmq, maxbq, nump, maxx ; 
  int iter, numiter  ; 
  int mq, bq, strand ;
  double ylike, ybase, yy ; 
  double *GX, **GE, **GD, *Y ;
  double *va, *vb, pdam ;
  static int loop = 0 ;

  ZALLOC(GX, 2, double) ;
  ZALLOC(Y, 4*4, double) ;
  GE = initarray_2Ddouble(4, 4, 0) ; 
  GD = initarray_2Ddouble(4, 4, 0) ; 

  numiter = emnumiter ;
  
  for (rcat = 0; rcat< numrgcats; ++rcat) { 
   ++loop ; 
   if (loop > numrgcats) fatalx("yukk\n") ;


  bighist = bighistlist[rcat] ;
  hipt = histinfo[rcat] ; 
  
   maxmq = hipt -> maxmq ; 
   maxbq = hipt -> maxbq ; 
   nump  = hipt -> nump ; 
   maxx  = 2*nump ;  
   if ((maxmq<0) || (maxbq<0))  { 
    printf("*** warning -- no data for RG cat %d\n", rcat) ;
    continue ;
   }
   gee = initarray_4Ddouble(maxmq+1, maxbq+1, 4, 4, 1.0e-10) ; 

   printf("gee alloc %d %d %d\n", rcat, loop, numrgcats) ;   fflush(stdout) ;
   gdd = initarray_4Ddouble(2,  maxx+1, 4, 4, 1.0e-10) ; 
   printf("gdd alloc\n") ;   fflush(stdout) ;

   for (iter = 1; iter <= numiter; ++iter) {
    vclear(gdamval, 1.0e-20, 2) ; 
    
    for (i=0; i<=maxmq; ++i) {  
     for (j=0; j <= maxbq; ++j) {
      clear2D(&gee[i][j], 4, 4, 1.0e-10) ;
     }
    }
   for (s=0; s<=1; ++s) { 
     for (j=0; j <= maxx; ++j) {
      clear2D(&gdd[s][j], 4, 4, 1.0e-10) ;
     }
   }
    ylike = 0 ; 
    for (mq=0; mq<=maxmq; ++mq) { 
     for (bq=0; bq<=maxbq; ++bq) { 
      for (strand=0; strand <2 ; ++strand) {
        for (x=0; x<=maxx; ++x) { 
         eef = eex[rcat][mq][bq] ;
         ddf = ddx[rcat][strand][x] ;
         damval = pdamval[rcat] ;

       t = htab[mq][bq][strand][x][rcat] ;
       if (t<0) continue ;
       hpt = bighist[t] ;
       floatit(Y, hpt -> hist, 16) ; 
       ylike += dogam(GX, GE, GD, damval, eef, ddf, Y) ;
       vvp(gdamval, gdamval, GX, 2) ; 
       plus2D(gee[mq][bq], gee[mq][bq], GE, 4, 4) ;
       plus2D(gdd[strand][x], gdd[strand][x], GD, 4, 4) ;

     }}}} 

    if (iter==1) ybase = ylike ;  
    ylike -= ybase ;  
    t = iter % 10 ;
    if ((iter == numiter) || (t==0)) {
     printf("rcat: %2d ", rcat) ;
     printf("iter: %4d  damval: %15.9f like:  %15.3f", iter, pdamval[rcat], ylike) ;
     printnl() ;
     fflush(stdout) ;
    }
    
// now reestimates
    bal1(gdamval, 2) ;
    pdamval[rcat] = gdamval[1] ; 
    for (mq=0; mq<=maxmq; ++mq) { 
     for (bq=0; bq<=maxbq; ++bq) { 
        eef = eex[rcat][mq][bq] ;
        for (r=0; r<=3; ++r)  {  
         va = eef[r] ; 
         vb = gee[mq][bq][r] ; 
         bal1(vb, 4) ; 
         copyarr(vb, va, 4) ;
         vclip(va, va, .0001, .999, 4) ;
         bal1(va, 4) ; 
        }
       }}
      for (strand=0; strand <2 ; ++strand) {
        for (x=0; x<=maxx; ++x) { 
         ddf = ddx[rcat][strand][x] ;
         for (r=0; r<=3; ++r)  {  
          copyarr(vb, va, 4) ;
          va = ddf[r] ; 
          vb = gdd[strand][x][r] ; 
          bal1(vb, 4) ; 
          copyarr(vb, va, 4) ;
         }
        }}
   }
   mq = maxmq ; bq = maxbq ;  scat = rcat ;  strand = 0        ;  
   printf("em estimates\n") ;
   eef = eex[scat][mq][bq] ;
   printf("eex\n") ; printmatl2D(eef, 4, 4) ;
   for (x=0 ; x<=maxx; ++x)  {  
    pp = atop(x, nump) ;
    ddf = ddx[scat][0][x] ;
    printf("ddx %d %d %d\n", x, 0, pp) ; printmatl2D(ddf, 4, 4) ;
    ddf = ddx[scat][1][x] ;
    printf("ddx %d %d %d\n", x, 1, pp) ; printmatl2D(ddf, 4, 4) ;
   }


// destructor
    free4D(&gee, maxmq+1, maxbq+1, 4) ;
    free4D(&gdd, 2, maxx+1, 4) ;
   
 }  // rcat

 fflush(stdout) ;

/**
 rcat = 0 ; 
 damval = pdamval[rcat] ; 
 ybase = emlike(&damval, eex, ddx, 
  histinfo, bighistlist, htab,  numrgcats,  0,  rcat)  ;

 for (x=1; x<=99; ++x)  {  
  break ; 
  damval = (double) x / 100.0 ; 
  yy = emlike(&damval, eex, ddx, 
   histinfo, bighistlist, htab,  numrgcats,  ybase,  rcat)  ;
  printf("zzs  damval: %15.9f like:  %15.3f", damval, yy) ;
  printnl() ;
 
 }
*/
 

 free(GX) ; 
 free(Y) ; 
 free2D(&GE, 4) ;
 free2D(&GD, 4) ;

 return ylike ;
}
  
void setmaxx(HISTINFO **histinfo, int numrgcats, int *pmaxmq, int *pmaxbq, int *pnump) 
{

  int rcat, maxmq, maxbq, nump ; 
  HISTINFO *hipt ; 
  HISTS **bighist ;  
  
  maxmq = maxbq = nump = 0 ; 

  for (rcat = 0; rcat< numrgcats; ++rcat) { 

  hipt = histinfo[rcat] ; 
  
   maxmq = MAX(maxmq, hipt -> maxmq) ; 
   maxbq = MAX(maxbq, hipt -> maxbq) ; 
   nump  = MAX(nump, hipt -> nump) ; 

  }
  *pmaxmq = maxmq ; 
  *pmaxbq = maxbq ; 
  *pnump = nump ;

}

void emcreatepars(double *pdamvals, double *****eex, double *****ddx, 
 HISTINFO **histinfo, HISTS ***bighistlist, int *****htab, int numrgcats, int xtype)  
{
  HISTINFO *hipt, *hinfo ; 
  HISTS **bighist, *hpt, *hpt0 ;  

  int k, i, j, x, s, rcat, t, tt  ; 
  int r, c ;
  int maxmq, maxbq, nump, maxx ; 
  
  setmaxx(histinfo, numrgcats, &maxmq, &maxbq, &nump) ; 
  maxx = 2*nump ; 

  t = numrgcats * 2 ;  // numrgcats... stores indel pars
  
   if (zee == NULL) { 
    zee = initarray_5Ddouble(t, maxmq+1, maxbq+1, 4, 4, 0) ; 
    zdd = initarray_5Ddouble(t, 2, maxx+1,  4, 4, 0) ; 
    ZALLOC(zdamvals, t, double) ;
   }
   tt = 0 ;               
   if (xtype == 1) tt = numrgcats; 
   copyarr5D(eex, zee + tt, numrgcats, maxmq +1, maxbq +1, 4, 4) ;
   copyarr5D(ddx, zdd + tt, numrgcats, 2, maxx +1,  4, 4) ;
   copyarr(pdamvals,  zdamvals +tt, numrgcats)  ; // damagescore
   vclip(zdamvals+tt,  zdamvals +tt, .001, .999, numrgcats)  ; 
}

double drfest(double *cprob, double *baseprob, int rcat, double *cscore) 
{
// cprob is base distribution from "prev" or "post" 
// baseprob is overall base distribution 



  double  w2[4], ww[4], adis[4], bdis[4]  ;
  int s, t ; 
  double y, yl, yprob[4], theta ; 
  double damval ;

  damval = zdamvals[rcat] ; 
  copyarr(baseprob, adis, 4) ; 
  bal1(adis, 4) ; 
  copyarr(cprob, ww, 4) ; 
  bal1(ww,  4) ; 
  vst(w2, adis, 1-damval, 4) ; 
  vvm(bdis, ww, w2, 4) ; 
  vclip(bdis, bdis, .00001, 999, 4) ; 
  bal1(bdis, 4) ;
  
  printf("zzabdis\n") ; 
  printmat(adis, 1, 4) ; 
  printmat(cprob, 1, 4) ; 
  printmat(bdis, 1, 4) ;

  vvd(ww, bdis, adis, 4) ; 
  vlog(cscore, ww, 4) ;

  y = vdot(adis, cscore, 4) ; // will be negative
  return y ;

}