cleanup, renamed Likelihood to likelihood

PosteriorProb.cpp

@@ -21,7 +21,7 @@
 #include "misc.h"
 #include "recalibration.h"
 #include "ngsBriggs.h"
-#include "Likelihood.h"
+#include "likelihood.h"
 #include "PosteriorProb.h"
 
 typedef unsigned char uchar;
@@ -949,7 +949,7 @@ bam_hdr_t* calc_pp_pmd_prob(char *refName,char *ifname, char* ofname, int mapped
     int refId=-1;
     double num = 0.0;
     size_t max_site;
-    nproc1 = 0;
+    int nproc1 = 0;
 
     uchar * indref = NULL;
     while(((ret=sam_read1(in,hdr,b)))>0) {

Likelihood.cpp → likelihood.cpp

@@ -1,4 +1,4 @@
-#include "Likelihood.h"
+#include "likelihood.h"
 #include <cstring>
 #include <cstdlib>
 #include <cstdio>
@@ -35,46 +35,48 @@ double MIN0 = 1e-8;
 extern tsk_struct *my_tsk_struct;
 
 // Naive likelihood without nick frequencies
-double loglike(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA){
-    double lambda = x[0];
-    double delta = x[1];
-    double delta_s = x[2];
-    double ll = 0;
-    for(int i=0; i<MAXLENGTH; i++){
-        ll += scaleCT[i]*(freqCT[i]*log(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*log(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*log(pow(1-lambda,i+1)/2*delta_s)+(1-freqGA[i])*log(1-pow(1-lambda,i+1)/2*delta_s));
-    }
-    return -ll;
+double loglike(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA,int &counter){
+  counter++;
+  double lambda = x[0];
+  double delta = x[1];
+  double delta_s = x[2];
+  double ll = 0;
+
+  for(int i=0; i<MAXLENGTH; i++){
+    ll += scaleCT[i]*(freqCT[i]*log(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*log(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*log(pow(1-lambda,i+1)/2*delta_s)+(1-freqGA[i])*log(1-pow(1-lambda,i+1)/2*delta_s));
+  }
+  return -ll;
 }
 
 double b_loglike(const double *x, const void *ptr){
-    ncalls++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    return loglike(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA);
+  wrapOne *wo =(wrapOne *) ptr;
+  return loglike(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA,wo->counter[0]);
 }
 
-void loglike_grad(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA){
-    double lambda = x[0];
-    double delta = x[1];
-    double delta_s = x[2];
-    y[0] = 0;
-    y[1] = 0;
-    y[2] = 0;
-    for(int i=0; i<MAXLENGTH; i++){
-        y[0] -= scaleCT[i]*(freqCT[i]*(-(i+1)*pow(1-lambda,i)/2*(delta_s-delta))/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*((i+1)*pow(1-lambda,i)/2*(delta_s-delta))/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*(-(i+1))/(1-lambda)+(1-freqGA[i])*((i+1)*pow(1-lambda,i)/2*delta_s)/(1-pow(1-lambda,i+1)/2*delta_s));
-        y[1] -= scaleCT[i]*(freqCT[i]*(1-pow(1-lambda,i+1)/2)/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*(-1+pow(1-lambda,i+1)/2)/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)));
-        y[2] -= scaleCT[i]*(freqCT[i]*(pow(1-lambda,i+1)/2)/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*(-pow(1-lambda,i+1)/2)/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*1/(delta_s)+(1-freqGA[i])*(-pow(1-lambda,i+1)/2)/(1-pow(1-lambda,i+1)/2*delta_s));
-    }
+void loglike_grad(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA,int &counter){
+  counter++;
+  double lambda = x[0];
+  double delta = x[1];
+  double delta_s = x[2];
+  y[0] = 0;
+  y[1] = 0;
+  y[2] = 0;
+  for(int i=0; i<MAXLENGTH; i++){
+    y[0] -= scaleCT[i]*(freqCT[i]*(-(i+1)*pow(1-lambda,i)/2*(delta_s-delta))/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*((i+1)*pow(1-lambda,i)/2*(delta_s-delta))/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*(-(i+1))/(1-lambda)+(1-freqGA[i])*((i+1)*pow(1-lambda,i)/2*delta_s)/(1-pow(1-lambda,i+1)/2*delta_s));
+    y[1] -= scaleCT[i]*(freqCT[i]*(1-pow(1-lambda,i+1)/2)/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*(-1+pow(1-lambda,i+1)/2)/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)));
+    y[2] -= scaleCT[i]*(freqCT[i]*(pow(1-lambda,i+1)/2)/(delta+pow(1-lambda,i+1)/2*(delta_s-delta))+(1-freqCT[i])*(-pow(1-lambda,i+1)/2)/(1-delta-pow(1-lambda,i+1)/2*(delta_s-delta)))+scaleGA[i]*(freqGA[i]*1/(delta_s)+(1-freqGA[i])*(-pow(1-lambda,i+1)/2)/(1-pow(1-lambda,i+1)/2*delta_s));
+  }
 }
 
 void b_loglike_grad(const double *x,double *y,const void*ptr){
-    ncalls_grad++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    loglike_grad(x,y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA);
+  wrapOne *wo =(wrapOne *) ptr;
+  loglike_grad(x,y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA,wo->counter[1]);
 }
 
 
 //Consider contamination rate eps + sequencing error
-double loglike_complex3_full_b(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps){
+double loglike_complex3_full_b(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps,int &counter){
+  counter++;
     double lambda = x[0];
     double delta = x[1];
     double delta_s = x[2];
@@ -160,7 +162,8 @@ double loglike_complex3_full_b(const double *x, double * freqCT, double * freqGA
 
 
 // Considering Seq Error + Non-biotin model: symmetric pattern
-double loglike_complex3_full_nb(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps){
+double loglike_complex3_full_nb(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps,int &counter){
+  counter++;
     double lambda = x[0];
     double delta = x[1];
     double delta_s = x[2];
@@ -248,18 +251,18 @@ double loglike_complex3_full_nb(const double *x, double * freqCT, double * freqG
 }
 
 double b_loglike_complex3_full(const double *x, const void *ptr){
-    ncalls++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    return loglike_complex3_full_b(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps);
+
+  wrapOne *wo =(wrapOne *) ptr;
+  return loglike_complex3_full_b(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps,wo->counter[0]);
 }
 
 double nb_loglike_complex3_full(const double *x, const void *ptr){
-    ncalls++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    return loglike_complex3_full_nb(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps);
+  wrapOne *wo =( wrapOne *) ptr;
+  return loglike_complex3_full_nb(x, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps,wo->counter[0]);
 }
 
-void loglike_complex3_grad_full_b(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps){
+void loglike_complex3_grad_full_b(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps,int &counter){
+  counter++;
     double lambda = x[0];
     double delta = x[1];
     double delta_s = x[2];
@@ -444,7 +447,8 @@ void loglike_complex3_grad_full_b(const double *x,double *y, double * freqCT, do
     }
 }
 
-void loglike_complex3_grad_full_nb(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps){
+void loglike_complex3_grad_full_nb(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double* LEN, double* freqLEN, double eps,int &counter){
+  counter++;
     double lambda = x[0];
     double delta = x[1];
     double delta_s = x[2];
@@ -635,15 +639,13 @@ void loglike_complex3_grad_full_nb(const double *x,double *y, double * freqCT, d
 
 
 void b_loglike_complex3_grad_full(const double *x,double *y,const void *ptr){
-    ncalls_grad++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    loglike_complex3_grad_full_b(x, y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps);
+  wrapOne *wo =( wrapOne *) ptr;
+    loglike_complex3_grad_full_b(x, y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps,wo->counter[1]);
 }
 
 void nb_loglike_complex3_grad_full(const double *x,double *y,const void *ptr){
-    ncalls_grad++;
-    const wrapOne *wo =(const wrapOne *) ptr;
-    loglike_complex3_grad_full_nb(x, y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps);
+  wrapOne *wo =(wrapOne *) ptr;
+  loglike_complex3_grad_full_nb(x, y, wo->freqCT, wo->freqGA, wo->scaleCT, wo->scaleGA, wo->seqError, wo->BinNum, wo->Bin_Frag_len, wo->Bin_Frag_freq, wo->Contam_eps,wo->counter[1]);
 }
 
 void loglike_complex3_hessian_full_b(const double *x, double ** z, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA, double * seqError, int BinNum, double * LEN, double * freqLEN, double eps){

Likelihood.h → likelihood.h

@@ -11,7 +11,7 @@ extern double PhredErrorAThird[255];
 extern int tsk_nthreads;
 
 
-typedef struct{
+typedef struct{//this struct is used in recalibration
   int len_limit;
   int len_min;
   double eps;
@@ -23,6 +23,7 @@ typedef struct{
   double **llh_result_hess; //Add the hessian matrix feature
   int threadid;
   double **mat;
+  int counter[2];//<- ncall,ncall_gradient
 }tsk_struct;
 
 
@@ -36,11 +37,12 @@ typedef struct{
   double *Bin_Frag_len;
   double *Bin_Frag_freq;
   double Contam_eps;
+  int counter[2];//<- ncall,ncall_gradient
 }wrapOne;
 
 //below are naive versions with out the nick. 
-void loglike_grad(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA);
-double loglike(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA);
+void loglike_grad(const double *x,double *y, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA,int &counter);
+double loglike(const double *x, double * freqCT, double * freqGA, double * scaleCT, double * scaleGA,int &counter);
 double b_loglike(const double *x, const void *);
 
 

misc.cpp

@@ -494,8 +494,6 @@ void parse_sequencingdata1(char *refName,char *fname, int mapped_only,int se_onl
         if(mapq!=-1 && b->core.qual<mapq)
             continue;
 
-	nproc1++;
-
         //then we simply write it to the output
         memset(reconstructedRef,0,512);
         memset(myread,'N',512);

ngsBriggs.cpp

@@ -15,7 +15,7 @@
 #include "read_all_reads.h"
 #include "misc.h"
 #include "recalibration.h"
-#include "Likelihood.h"
+#include "likelihood.h"
 #include "PosteriorProb.h"
 #include "ngsBriggs_cli.h"
 #include "ngsBriggs.h"
@@ -24,10 +24,7 @@
 
 // definining all global variables used across multiple scripts
 
-int nproc1 = 0;//number of reads processed maybe not used
 double l_check = 15;
-int ncalls =0;
-int ncalls_grad =0;
 tsk_struct *my_tsk_struct = NULL;
 
 int tsk_nthreads = -1;
@@ -36,7 +33,7 @@ int tsk_nthreads = -1;
 // defining our main ngsBriggs function
 int main(int argc, char **argv){
 
-  double** mm5p, **mm3p;// **dmm5p, **dmm3p;
+  double** mm5p, **mm3p;
   double Tol = 1.0E-8; // Tolerance
   double **deamRateCT;
   double **deamRateGA;
@@ -288,14 +285,13 @@ int main(int argc, char **argv){
     double ubd[3] = {1,1,1};
     int nbd[3] = {2,2,2};
     double invec1[3] = {0.2,0.1,0.1};
-    ncalls=0;
-    ncalls_grad = 0;
     wrapOne wo;
     wo.freqCT = freqCT;
     wo.freqGA = freqGA;
     wo.scaleCT = scaleCT;
     wo.scaleGA = scaleGA;
     wo.seqError = seqError;
+    wo.counter[0] = wo.counter[1] = 0;
     double withgrad = findmax_bfgs(3,invec1,(void *)&wo,b_loglike,b_loglike_grad,lbd,ubd,nbd,-1);
 
     double lbd1[4] = {1e-8,1e-8,1e-8,1e-8};
@@ -330,14 +326,12 @@ int main(int argc, char **argv){
     FragArrayBin(number, BinNum, Frag_len, Frag_freq, Bin_Frag_len, Bin_Frag_freq);
     //    invec2 = {lambda,deltad,deltas,nu};
     double invec2[4] = {invec1[0],invec1[1],invec1[2],0.01};
-    ncalls=0;
-    ncalls_grad = 0;
     double withgrad2 = 0;
     wo.BinNum = BinNum;
     wo.Bin_Frag_len = Bin_Frag_len;
     wo.Bin_Frag_freq = Bin_Frag_freq;
     wo.Contam_eps = Contam_eps;
-
+    wo.counter[0] =wo.counter[1] = 0;
     if (model==0){
       fprintf(stderr,"%s\n","The chosen model is a biotin model.");
       withgrad2 = findmax_bfgs(4,invec2,(void *) &wo,b_loglike_complex3_full,b_loglike_complex3_grad_full,lbd1,ubd1,nbd1,-1);
@@ -360,7 +354,7 @@ int main(int argc, char **argv){
       fprintf(stderr,"%s","The chosen model is non-biotin model, ");
     }
     fprintf(stderr,"the inferred parameters are given as follows:\n");
-    fprintf(stderr,"lambda: %f (%f,%f), delta: %f (%f,%f), delta_s: %f (%f,%f), nu: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", invec2[0], invec2[0]-1.96*stdvec2[0], invec2[0]+1.96*stdvec2[0] , invec2[1], invec2[1]-1.96*stdvec2[1], invec2[1]+1.96*stdvec2[1], invec2[2], invec2[2]-1.96*stdvec2[2], invec2[2]+1.96*stdvec2[2], invec2[3], invec2[3]-1.96*stdvec2[3], invec2[3]+1.96*stdvec2[3], ncalls,ncalls_grad, withgrad2);
+    fprintf(stderr,"lambda: %f (%f,%f), delta: %f (%f,%f), delta_s: %f (%f,%f), nu: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", invec2[0], invec2[0]-1.96*stdvec2[0], invec2[0]+1.96*stdvec2[0] , invec2[1], invec2[1]-1.96*stdvec2[1], invec2[1]+1.96*stdvec2[1], invec2[2], invec2[2]-1.96*stdvec2[2], invec2[2]+1.96*stdvec2[2], invec2[3], invec2[3]-1.96*stdvec2[3], invec2[3]+1.96*stdvec2[3], wo.counter[0],wo.counter[1], withgrad2);
 
     // Output count table
     char *otabname = mypars->otab;
@@ -430,7 +424,7 @@ int main(int argc, char **argv){
 	  ksprintf(kstr,"%s","The chosen model is non-biotin model, ");
 	}
         ksprintf(kstr,"the inferred parameters are given as follows:\n");
-        ksprintf(kstr,"lambda: %f (%f,%f), delta: %f (%f,%f), delta_s: %f (%f,%f), nu: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", invec2[0], invec2[0]-1.96*stdvec2[0], invec2[0]+1.96*stdvec2[0] , invec2[1], invec2[1]-1.96*stdvec2[1], invec2[1]+1.96*stdvec2[1], invec2[2], invec2[2]-1.96*stdvec2[2], invec2[2]+1.96*stdvec2[2], invec2[3], invec2[3]-1.96*stdvec2[3], invec2[3]+1.96*stdvec2[3], ncalls,ncalls_grad, withgrad2);
+        ksprintf(kstr,"lambda: %f (%f,%f), delta: %f (%f,%f), delta_s: %f (%f,%f), nu: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", invec2[0], invec2[0]-1.96*stdvec2[0], invec2[0]+1.96*stdvec2[0] , invec2[1], invec2[1]-1.96*stdvec2[1], invec2[1]+1.96*stdvec2[1], invec2[2], invec2[2]-1.96*stdvec2[2], invec2[2]+1.96*stdvec2[2], invec2[3], invec2[3]-1.96*stdvec2[3], invec2[3]+1.96*stdvec2[3], wo.counter[0],wo.counter[1], withgrad2);
 	//my_bgzf_write(fp,kstr->s,kstr->l);
         assert(bgzf_write(fp,kstr->s,kstr->l)==kstr->l);
         kstr->l = 0;
@@ -484,9 +478,7 @@ int main(int argc, char **argv){
       double **mat = read_all_reads(in,hdr,seq_ref,len_limit,invec2[0],invec2[1],invec2[2],invec2[3],Tol,ndim,model);
 
         //tsk stop
-        ncalls = 0;
-        ncalls_grad = 0;
-        //double lbd2[4] = {30,1e-8,30,1e-8};
+      //double lbd2[4] = {30,1e-8,30,1e-8};
         double lbd2[4] = {30,15,30,1};
         double ubd2[4] = {(double)len_limit-1,100,(double)len_limit-1,100};
         int nbd2[4] = {2,2,2,2};
@@ -508,6 +500,7 @@ int main(int argc, char **argv){
             my_tsk_struct[ii].threadid = ii;
             my_tsk_struct[ii].from = ii==0?0:my_tsk_struct[ii-1].to;
             my_tsk_struct[ii].to =   my_tsk_struct[ii].from+mynsites;
+	    my_tsk_struct[ii].counter[0] =  my_tsk_struct[ii].counter[1] = 0;
         }
         my_tsk_struct[tsk_nthreads-1].to = ndim;
         double withgrad3;
@@ -526,7 +519,7 @@ int main(int argc, char **argv){
             stdpar[i] = sqrt(covpar[i][i]);
             //cout<<z4[i][i]<<" std "<<stdvec4[i]<<"\n";
         }
-        fprintf(stderr,"mu_anc: %f (%f,%f), sigma_anc: %f (%f,%f), mu_mod: %f (%f,%f), sigma_mod: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", distparam[0], distparam[0]-1.96*stdpar[0], distparam[0]+1.96*stdpar[0], distparam[1], distparam[1]-1.96*stdpar[1], distparam[1]+1.96*stdpar[1], distparam[2], distparam[2]-1.96*stdpar[2], distparam[2]+1.96*stdpar[2], distparam[3], distparam[3]-1.96*stdpar[3], distparam[3]+1.96*stdpar[3], ncalls,ncalls_grad, withgrad3);
+        fprintf(stderr,"mu_anc: %f (%f,%f), sigma_anc: %f (%f,%f), mu_mod: %f (%f,%f), sigma_mod: %f (%f,%f), nfunctioncalls: %d ngradientcalls: %d, llh: %f.\n", distparam[0], distparam[0]-1.96*stdpar[0], distparam[0]+1.96*stdpar[0], distparam[1], distparam[1]-1.96*stdpar[1], distparam[1]+1.96*stdpar[1], distparam[2], distparam[2]-1.96*stdpar[2], distparam[2]+1.96*stdpar[2], distparam[3], distparam[3]-1.96*stdpar[3], distparam[3]+1.96*stdpar[3], my_tsk_struct[0].counter[0], my_tsk_struct[0].counter[1], withgrad3);
         sam_hdr_destroy(hdr);
     }
 

ngsBriggs.h

@@ -3,11 +3,8 @@
 extern int MAXLENGTH;
 extern double l_check;
 
-extern int ncalls;
-extern int ncalls_grad;
 extern char refToChar[256];
 extern char com[5];
-extern int nproc1;//number of reads processed not used
 
 #endif
 

read_all_reads.cpp

@@ -3,7 +3,7 @@
 #include <cassert>
 #include <htslib/bgzf.h>
 
-#include "Likelihood.h"
+#include "likelihood.h"
 #include "read_all_reads.h"
 #include "profile.h"
 #include "ngsBriggs.h"