MQM code: use R_alloc/S_alloc rather than Calloc/R_chk_calloc

 - Trying to fix memory leaks reported by valgrind; simplest solution is
   to just switch to using memory that R controls (and frees on our
   behalf).

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: qtl
-Version: 1.34-12
+Version: 1.34-13
 Date: 21 Oct 2014
 Title: Tools for analyzing QTL experiments
 Author: Karl W Broman <kbroman@biostat.wisc.edu> and Hao Wu, with

src/mqmaugment.cpp

@@ -241,9 +241,6 @@ int mqmaugmentfull(MQMMarkerMatrix* markers,int* nind, int* augmentednind, ivect
       (*nind)= (*nind)+(current_leftover_ind);
       debug_trace("nind:%d,naugmented:%d",(*nind)+(current_leftover_ind),(*augmentednind)+(current_leftover_ind));
       Rprintf("INFO: VALGRIND MEMORY DEBUG BARRIERE TRIGGERED\n", "");
-      delMQMMarkerMatrix(newmarkerset, nmark);    // Free the newmarkerset, this can only be done here since: (*markers) = newmarkerset_all;
-      // Free(new_y_all);
-      // Free(new_ind_all);
     }else{
       if(ind_still_left && augment_strategy == 3){
         if(verbose) Rprintf("INFO: Dropping %d augment_strategy individuals from further analysis\n",ind_still_left);
@@ -254,11 +251,6 @@ int mqmaugmentfull(MQMMarkerMatrix* markers,int* nind, int* augmentednind, ivect
       (*markers) = newmarkerset;
     }
     if(verbose) Rprintf("INFO: Done with augmentation\n");
-    // Free(new_y);                                // Free vector indicating the new phenotypes
-    // Free(new_ind);                              // Free vector indicating the new individuals
-    Free(succes_ind);                           // Free vector indicating the result of round 1 - augmentation
-    Free(position);                             // Free the positions of the markers
-    Free(r);                                    // Free the recombination frequencies
     return 1;
 }
 
@@ -631,12 +623,6 @@ int mqmaugment(const MQMMarkerMatrix marker, const vector y,
   if (verbose) fatal("Recall procedure with larger value for augmentation parameters or increase the parameter minprob\n");
   retvalue = 0;
 cleanup:
-  Free(newy);
-  delMQMMarkerMatrix(newmarker, Nmark+1); //Free(newmarker);
-  Free(newind);
-  Free(newprob);
-  Free(newprobmax);
-  Free(imarker);
   return retvalue;
 }
 
@@ -747,9 +733,6 @@ void R_mqmaugment(int *geno, double *dist, double *pheno, int *auggeno,
     }
     fatal("Data augmentation failed", "");
   }
-  delMQMMarkerMatrix(markers,*Nmark); // [Danny:] This looked suspicious, we were leaking memory here because we didn't clean it
-  Free(mapdistance);
-  Free(chr);
 
   PutRNGstate(); // write R's RNG seed
 

src/mqmdatatypes.cpp

@@ -29,9 +29,9 @@
 
 #include "mqm.h"
 
-/* 
+/*
  * Determine the experimental cross type from the R/qtl dataset. Returns the
- * type. 
+ * type.
  */
 
 MQMCrossType determine_MQMCross(const int Nmark, const int Nind, const int **Geno, const RqtlCrossType rqtlcrosstype) {
@@ -101,15 +101,14 @@ void change_coding(int *Nmark, int *Nind, int **Geno, MQMMarkerMatrix markers, c
   }
 }
 
-/* 
+/*
  * Allocate a memory block using the 'safe' R method calloc_init, but with
  * guarantee all data has been zeroed
  */
 
 void *calloc_init(size_t num, size_t size) {
   void *buf;
-  buf = R_chk_calloc(num,size);
-  if (buf) memset(buf,0,num*size);
+  buf = S_alloc(num,size);
   return buf;
 }
 
@@ -189,30 +188,6 @@ MQMMarkerMatrix newMQMMarkerMatrix(int rows, int cols) {
   return m;
 }
 
-void freevector(void *v) {
-  Free(v);
-}
-
-void freematrix(void **m, size_t rows) {
-  for (size_t i = 0; i < rows; i++) {
-    Free(m[i]);
-  }
-  Free(m);
-}
-
-void delmatrix(matrix m, size_t rows) {
-  freematrix((void**)m,rows);
-}
-
-void delcmatrix(cmatrix m, size_t rows) {
-  freematrix((void **)m,rows);
-}
-
-void delMQMMarkerMatrix(MQMMarkerMatrix m, size_t rows) {
-  freematrix((void **)m,rows);
-}
-
-
 void copyvector(vector vsource, vector vdestination, int dim) {
   for(int i=0; i<dim; i++) {
     vdestination[i]= vsource[i];
@@ -226,4 +201,3 @@ void copyvector(vector vsource, vector vdestination, int dim) {
   void debug_trace(const char*, ...){}
   void info(const char* s, ...){ Rprintf("INFO: "); Rprintf(s); Rprintf("\n"); }
 #endif
-

src/mqmdatatypes.h

@@ -59,10 +59,10 @@ const double POSITIONUNKNOWN = 999.0;
 // Marker genotypes (scored at marker)
 /*
 const unsigned char MAA       = '0';  // Homozygous AA
-const unsigned char MH        = '1';  // Heterozygous AB 
+const unsigned char MH        = '1';  // Heterozygous AB
 const unsigned char MBB       = '2';  // Homozygous BB
 const unsigned char MNOTAA    = '3';  // Not AA
-const unsigned char MNOTBB    = '4';  // Not BB 
+const unsigned char MNOTBB    = '4';  // Not BB
 const unsigned char MMISSING  = '9';  // Unknown (marker genotype missing)
 const unsigned char MUNUSED   = '-';  // Unused parameter
 */
@@ -100,11 +100,6 @@ void printmatrix(matrix m, int rows, int cols);
 void printcmatrix(cmatrix m, int rows, int cols);
 cmatrix newcmatrix(int rows, int cols);
 MQMMarkerMatrix newMQMMarkerMatrix(int rows, int cols);
-void freematrix(void **m, size_t rows);
-void freevector(void *v);
-void delmatrix(matrix m, size_t rows);
-void delcmatrix(cmatrix m, size_t rows);
-void delMQMMarkerMatrix(MQMMarkerMatrix m,size_t rows);
 void copyvector(vector vsource, vector vdestination, int dim);
 void *calloc_init(size_t num, size_t size);
 

src/mqmeliminate.cpp

@@ -53,7 +53,7 @@ double backward(int Nind, int Nmark, cvector cofactor, MQMMarkerMatrix marker,
   double savelogL    = logLfull;
   double maxlogL     = logLfull-10000.0;
   bool warned        = false;
-  
+
   if (verbose) Rprintf("INFO: Backward elimination of cofactors started\n");
   for (int j=0; j<Nmark; j++) {
     (*newcofactor)[j]= cofactor[j];
@@ -96,7 +96,7 @@ double backward(int Nind, int Nmark, cvector cofactor, MQMMarkerMatrix marker,
       savelogL= maxlogL;
       (*newcofactor)[dropj]= MNOCOF;
       Ncof-=1;
-      if(verbose) 
+      if(verbose)
         Rprintf("INFO: Marker %d is dropped, resulting in reduced model logL = %.3f\n",(dropj+1),ftruncate3(savelogL));
     } else if ( ((*newcofactor)[dropj]==MBB) && (F1> 2.0*(savelogL-maxlogL)) ) {
       savelogL= maxlogL;
@@ -121,6 +121,5 @@ double backward(int Nind, int Nmark, cvector cofactor, MQMMarkerMatrix marker,
   //Map using the model
   maxF = mapQTL(Nind, Nmark, cofactor, (*newcofactor), marker, position, (*mapdistance), y, r, ind, Naug, variance, 'n',
                 informationcontent,Frun,run,REMLorML,fitQTL,dominance, em, windowsize, stepsize, stepmin, stepmax,crosstype,verbose);
-  Free(logL);
   return maxF;
 }

src/mqmmapqtl.cpp

@@ -29,21 +29,20 @@
 
 #include "mqm.h"
 
-/* 
+/*
  * mapQTL moves a QTL along the chromosome and calculated at each map position
  * the QTL likelihood. Uses either all cofactors, or selected cofactors only
  */
-double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor, 
+double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
               MQMMarkerMatrix marker, cvector position, vector mapdistance, vector y,
               vector r, ivector ind, int Naug, double variance, char
               printoutput, vector *informationcontent, matrix *Frun, int run,
               char REMLorML, bool fitQTL, bool dominance, int em, double
-              windowsize, double stepsize, double stepmin, double stepmax, 
+              windowsize, double stepsize, double stepmin, double stepmax,
               MQMCrossType crosstype, int verbose) {
   if(verbose) Rprintf("INFO: mapQTL function called\n");
   int j, jj, jjj=0;
   int Nloci = Nmark+1;
-  vector Fy = newvector(Naug);
   cvector QTLcofactor       = newcvector(Nloci);
   cvector saveQTLcofactor   = newcvector(Nloci);
   bool warned  = false;
@@ -56,8 +55,6 @@ double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
   weight[0]= -1.0;
 
   /* fit QTL on top of markers (but: should also be done with routine QTLmixture() for exact ML) */
-  cvector newcofactor= newcvector(Nmark);
-  cvector direction = newcvector(Nmark);
   vector cumdistance = newvector(Nmark+1);
   double QTLlikelihood=0.0;
 
@@ -73,15 +70,14 @@ double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
   // cout << "please wait (mixture calculus may take quite a lot of time)" << endl;
   /* estimate variance in mixture model with all marker cofactors */
   // cout << "estimate variance in mixture model with all cofactors" << endl;
-  
+
   variance= -1.0;
   savelogL= 2.0*QTLmixture(marker, cofactor, r, position, y, ind, Nind, Naug, Nmark, &variance, em, &weight, REMLorML, fitQTL, dominance, crosstype, &warned, verbose);
   if (verbose) Rprintf("INFO: log-likelihood of full model = %f\n", savelogL/2);
 
   // augment data for missing QTL observations (x 3)
   fitQTL=true;
   int newNaug = 3 * Naug;
-  Free(weight);
   weight           = newvector(newNaug);
   weight[0]        = 1.0;
   vector weight0   = newvector(newNaug);
@@ -90,15 +86,15 @@ double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
   vector QTLr              = newvector(Nloci);
   vector QTLmapdistance    = newvector(Nloci);
   cvector QTLposition      = newcvector(Nloci);
-  MQMMarkerMatrix QTLloci  = (MQMMarkerMatrix)Calloc(Nloci, MQMMarkerVector);
+  MQMMarkerMatrix QTLloci  = (MQMMarkerMatrix)R_alloc(Nloci, sizeof(MQMMarkerVector));
 
   double moveQTL = stepmin;
   char nextinterval= 'n', firsttime='y';
   double maxF=0.0, savebaseNoQTLModel=0.0;
   int baseNoQTLModel=0, step=0;
 
   for (j=0; j<Nmark; j++) {
-    /* 	fit a QTL in two steps:
+    /*  fit a QTL in two steps:
     1. move QTL along marker interval j -> j+1 with steps of stepsize=20 cM, starting from -20 cM up to 220 cM
     2. all marker-cofactors in the neighborhood of the QTL are dropped by using cM='windows' as criterium
     */
@@ -261,7 +257,7 @@ double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
         if (run>0) (*Frun)[step][run]+= QTLlikelihood;
         else (*Frun)[step][0]+= QTLlikelihood;
 
-        /* 	Each individual has condition multilocus probabilities for being 0, 1 or 2 at the QTL.
+        /*  Each individual has condition multilocus probabilities for being 0, 1 or 2 at the QTL.
         Calculate the maximum per individu. Calculate the mean of this maximum, averaging over all individuals
         This is the information content plotted.
         */
@@ -287,21 +283,5 @@ double mapQTL(int Nind, int Nmark, cvector cofactor, cvector selcofactor,
   }
   fitQTL=false;
 
-  freevector(direction);
-  Free(info0);
-  Free(info1);
-  Free(info2);
-  Free(weight);
-  Free(weight0);
-  Free(QTLr);
-  Free(QTLposition);
-  Free(Fy);
-  Free(newcofactor);
-  Free(QTLcofactor);
-  Free(cumdistance);
-  Free(QTLmapdistance);
-  Free(QTLloci);
-  Free(saveQTLcofactor);
   return maxF; //QTLlikelihood;
 }
-

src/mqmmixture.cpp

@@ -131,9 +131,7 @@ double rmixture(MQMMarkerMatrix marker, vector weight, vector r,
       //Rprintf("r(%d)= %f -> %f\n", j, r[j], (*mapdistance)[j]);
     }
   }
-  if (verbose==1) Rprintf("INFO: Re-estimation of the genetic map took %d iterations, to reach a rdelta of %f\n", iem, rdelta); 
-  Free(indweight);
-  freevector(distance);
+  if (verbose==1) Rprintf("INFO: Re-estimation of the genetic map took %d iterations, to reach a rdelta of %f\n", iem, rdelta);
   return maximum;
 }
 
@@ -143,16 +141,16 @@ double rmixture(MQMMarkerMatrix marker, vector weight, vector r,
  * multilocus information, but assuming known recombination frequencies
 */
 double QTLmixture(MQMMarkerMatrix loci, cvector cofactor, vector r, cvector position,
-                  vector y, ivector ind, int Nind, int Naug, int Nloci, double *variance, 
+                  vector y, ivector ind, int Nind, int Naug, int Nloci, double *variance,
                   int em, vector *weight, const bool useREML,const bool fitQTL,const bool dominance, MQMCrossType crosstype, bool* warned, int verbose) {
-                  
+
   //debug_trace("QTLmixture called Nloci=%d Nind=%d Naug=%d, REML=%d em=%d fit=%d domi=%d cross=%c\n",Nloci,Nind,Naug,useREML,em,fitQTL,dominance,crosstype);
   //for (int i=0; i<Nloci; i++){ debug_trace("loci %d : recombfreq=%f\n",i,r[i]); }
   int iem= 0, i, j;
   bool warnZeroDist=false;
   bool biasadj=false;
   double oldlogL = -10000, delta=1.0, calc_i, Pscale=1.75;
-  
+
   vector indweight  = newvector(Nind);
   int newNaug       = ((!fitQTL) ? Naug : 3*Naug);
   vector Fy         = newvector(newNaug);
@@ -221,9 +219,9 @@ double QTLmixture(MQMMarkerMatrix loci, cvector cofactor, vector r, cvector posi
             Ploci[i+2*Naug]*= calc_i;
           }
         else if (cofactor[j]<=MCOF) // locus j+1 == QTL
-          for (i=0; i<Naug; i++) { // QTL==MAA || MH || MBB means: What is the prob of finding an MAA at J=1    
+          for (i=0; i<Naug; i++) { // QTL==MAA || MH || MBB means: What is the prob of finding an MAA at J=1
             calc_i =left_prob(r[j],loci[j][i],MAA,crosstype);     //calc_i = prob(loci, r, i, j, MAA, crosstype, 0);
-            Ploci[i]*= calc_i;      
+            Ploci[i]*= calc_i;
             calc_i = left_prob(r[j],loci[j][i],MH,crosstype);     //calc_i = prob(loci, r, i, j, MH, crosstype, 0);
             Ploci[i+Naug]*= calc_i;
             calc_i = left_prob(r[j],loci[j][i],MBB,crosstype);    //calc_i = prob(loci, r, i, j, MBB, crosstype, 0);
@@ -242,7 +240,7 @@ double QTLmixture(MQMMarkerMatrix loci, cvector cofactor, vector r, cvector posi
     }
   }
   if(warnZeroDist && verbose && !(*warned)){
-    Rprintf("WARNING: 0.0 from Probability calculation! Markers at same cMorgan but different genotype?\n"); 
+    Rprintf("WARNING: 0.0 from Probability calculation! Markers at same cMorgan but different genotype?\n");
     (*warned) = true;
   }
 //	Rprintf("INFO: Done fitting QTL's\n");
@@ -302,7 +300,7 @@ double QTLmixture(MQMMarkerMatrix loci, cvector cofactor, vector r, cvector posi
     delta= fabs(logL-oldlogL);
     oldlogL= logL;
   }
-  
+
   if ((useREML)&&(!varknown)) { // Bias adjustment after finished ML estimation via EM
     *variance=-1.0;
     biasadj=true;
@@ -342,10 +340,5 @@ double QTLmixture(MQMMarkerMatrix loci, cvector cofactor, vector r, cvector posi
     }
   }
   //for (i=0; i<Nind; i++){ debug_trace("IND %d Ploci: %f Fy: %f UNLOG:%f LogL:%f LogL-LogP: %f\n", i, Ploci[i], Fy[i], indL[i], log(indL[i]), log(indL[i])-logP); }
-  Free(Fy);
-  Free(Ploci);
-  Free(indweight);
-  Free(indL);
   return logL;
 }
-

src/mqmregression.cpp

@@ -188,7 +188,7 @@ double regression(int Nind, int Nmark, cvector cofactor, MQMMarkerMatrix marker,
   ludcmp(XtWX, dimx, indx, &d);
   lusolve(XtWX, dimx, indx, XtWY);
 
-  long double* indL = (long double *)Calloc(Nind, long double);
+  long double* indL = (long double *)R_alloc(Nind, sizeof(long double));
   int newNaug       = ((!fitQTL) ? Naug : 3*Naug);
   vector fit        = newvector(newNaug);
   vector resi       = newvector(newNaug);
@@ -286,14 +286,6 @@ double regression(int Nind, int Nmark, cvector cofactor, MQMMarkerMatrix marker,
   for (int i=0; i<Nind; i++) { //Sum up log likelihoods for each individual
     logL+= log(indL[i]);
   }
-  Free(indL);
-  Free(indx);
-  Free(xtQTL);
-  delmatrix(XtWX,dimx_alloc);
-  delcmatrix(Xt,dimx_alloc);
-  Free(fit);
-  Free(resi);
-  freevector((void *)XtWY);
   return (double)logL;
 }
 
@@ -339,8 +331,6 @@ void ludcmp(matrix m, int dim, ivector ndx, int *d) {
     temp=1.0/m[c][c];
     for (r=c+1; r<dim; r++) m[r][c]*=temp;
   }
-  Free(scale);
-  //Free(swap);
 }
 
 /* Solve the set of n linear equations AX=B.