[Genabel-commits] r1444 - pkg/ProbABEL/src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Dec 2 19:22:46 CET 2013


Author: lckarssen
Date: 2013-12-02 19:22:46 +0100 (Mon, 02 Dec 2013)
New Revision: 1444

Modified:
   pkg/ProbABEL/src/coxph_data.cpp
   pkg/ProbABEL/src/eigen_mematrix.cpp
   pkg/ProbABEL/src/extract-snp.cpp
   pkg/ProbABEL/src/main.cpp
   pkg/ProbABEL/src/main_functions_dump.cpp
   pkg/ProbABEL/src/main_functions_dump.h
   pkg/ProbABEL/src/reg1.cpp
   pkg/ProbABEL/src/regdata.cpp
   pkg/ProbABEL/src/usage.cpp
Log:
- Fixing a bunch of cpplint errors reported by Jenkins, mostly about white space, so only code layout. Also includes conversion of a few C-style casts to C++ style. 
- Replaced the tabs introduced by mistake in svn r.1433 with spaces


Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/coxph_data.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -55,8 +55,8 @@
     strata      = obj.strata;
     X           = obj.X;
     order       = obj.order;
-    gcount		=0;
-    freq		=0;
+    gcount      = 0;
+    freq        = 0;
     masked_data = new unsigned short int[nids];
 
     for (int i = 0; i < nids; i++)
@@ -96,7 +96,7 @@
     {
         //          X.put(1.,i,0);
         stime[i] = (phed.Y).get(i, 0);
-        sstat[i] = int((phed.Y).get(i, 1));
+        sstat[i] = static_cast<int>((phed.Y).get(i, 1));
         if (sstat[i] != 1 && sstat[i] != 0)
         {
             std::cerr << "coxph_data: status not 0/1 "
@@ -188,65 +188,68 @@
 }
 
 void coxph_data::update_snp(gendata &gend, const int snpnum) {
-	/**
-	 * This is the main part of the fix of bug #1846
-	 * (C) of the fix:
-	 *   UMC St Radboud Nijmegen,
-	 *   Dept of Epidemiology & Biostatistics,
-	 *   led by Prof. B. Kiemeney
-	 *
-	 * Note this sorts by "order"!!!
-	 * Here we deal with transposed X, hence last two arguments are swapped
-	 * compared to the other 'update_snp'
-	 * Also, the starting column-1 is not necessary for cox X therefore
-	 * 'ncov-j' changes to 'ncov-j-1'
-	 **/
-	//reset counter for frequency since it is a new snp
-	gcount=0;
-	freq=0.0;
-	for (int j = 0; j < ngpreds; j++) {
-		double *snpdata = new double[nids];
-		for (int i = 0; i < nids; i++) {
-			masked_data[i] = 0;
-		}
+    /**
+     * This is the main part of the fix of bug #1846
+     * (C) of the fix:
+     *   UMC St Radboud Nijmegen,
+     *   Dept of Epidemiology & Biostatistics,
+     *   led by Prof. B. Kiemeney
+     *
+     * Note this sorts by "order"!!!
+     * Here we deal with transposed X, hence last two arguments are swapped
+     * compared to the other 'update_snp'
+     * Also, the starting column-1 is not necessary for cox X therefore
+     * 'ncov-j' changes to 'ncov-j-1'
+     **/
 
-		gend.get_var(snpnum * ngpreds + j, snpdata);
+    // reset counter for frequency since it is a new snp
+    gcount = 0;
+    freq   = 0.0;
 
-		for (int i = 0; i < nids; i++) {
-			X.put(snpdata[i], (ncov - j - 1), order[i]);
-			if (std::isnan(snpdata[i])) {
-				masked_data[order[i]] = 1;
-				//snp not masked
-			} else {
-				// checck for first predicor
-				if (j == 0) {
-					gcount++;
-					if (ngpreds == 1) {
-						freq += snpdata[i] * 0.5;
-					} else if (ngpreds == 2) {
-						freq += snpdata[i];
-					}
-				} else if (j == 1) {
-					// add second genotype in two predicor data form
-					freq += snpdata[i] * 0.5;
-				}
-			}    //end std::isnan(snpdata[i]) snp
+    for (int j = 0; j < ngpreds; j++) {
+        double *snpdata = new double[nids];
+        for (int i = 0; i < nids; i++) {
+            masked_data[i] = 0;
+        }
 
-		}    //end i for loop
+        gend.get_var(snpnum * ngpreds + j, snpdata);
 
-		delete[] snpdata;
-	}    //end ngpreds loop
-	freq /= static_cast<double>(gcount); // Allele frequency
+        for (int i = 0; i < nids; i++) {
+            X.put(snpdata[i], (ncov - j - 1), order[i]);
+            if (std::isnan(snpdata[i])) {
+                masked_data[order[i]] = 1;
+                // snp not masked
+            } else {
+                // check for first predictor
+                if (j == 0) {
+                    gcount++;
+                    if (ngpreds == 1) {
+                        freq += snpdata[i] * 0.5;
+                    } else if (ngpreds == 2) {
+                        freq += snpdata[i];
+                    }
+                } else if (j == 1) {
+                    // add second genotype in two predicor data form
+                    freq += snpdata[i] * 0.5;
+                }
+            }    // end std::isnan(snpdata[i]) snp
+        }    // end i for loop
 
+        delete[] snpdata;
+    }    // end ngpreds loop
+
+    freq /= static_cast<double>(gcount);  // Allele frequency
 }
 
 
-
+/**
+ * update_snp() adds SNP information to the design matrix. This
+ * function allows you to strip that information from X again. This
+ * is used for example when calculating the null model.
+ *
+ */
 void coxph_data::remove_snp_from_X()
 {
-    // update_snp() adds SNP information to the design matrix. This
-    // function allows you to strip that information from X again.
-    // This is used for example when calculating the null model.
 
     if (ngpreds == 1)
     {
@@ -277,9 +280,10 @@
     //      delete order;
 }
 
+
 coxph_data coxph_data::get_unmasked_data()
 {
-    coxph_data to; // = coxph_data(*this);
+    coxph_data to;  // = coxph_data(*this);
 
     // filter missing data
     int nmeasured = 0;
@@ -290,6 +294,7 @@
             nmeasured++;
         }
     }
+
     to.nids = nmeasured;
     to.ncov = ncov;
     to.ngpreds = ngpreds;
@@ -331,6 +336,7 @@
     return (to);
 }
 
+
 coxph_reg::coxph_reg(coxph_data &cdatain)
 {
     coxph_data cdata = cdatain.get_unmasked_data();
@@ -342,6 +348,7 @@
     niter = 0;
 }
 
+
 void coxph_reg::estimate(coxph_data &cdatain, const int verbose,
                         int maxiter, double eps,
                         double tol_chol, const int model,
@@ -432,7 +439,7 @@
         MatrixXd imateigen = imat.data;
         VectorXd infs = ueigen.transpose() * imateigen;
         VectorXd betaeigen = beta.data;
-        if( infs.norm() > eps ||
+        if ( infs.norm() > eps ||
             infs.norm() > sqrt(eps) * betaeigen.norm() )
         {
             cerr << "Warning for " << snpinfo.name[cursnp]
@@ -457,9 +464,7 @@
             sebeta[i] = NAN;
             beta[i]   = NAN;
             loglik    = NAN;
-        }
-        else
-        {
+        } else {
             sebeta[i] = sqrt(imat.get(i, i));
             loglik = loglik_int[1];
         }

Modified: pkg/ProbABEL/src/eigen_mematrix.cpp
===================================================================
--- pkg/ProbABEL/src/eigen_mematrix.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/eigen_mematrix.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -66,19 +66,19 @@
         exit(1);
     }
     int column = i % ncol;
-    int row = (int) floor((double) i / ncol);
+    int row = static_cast<int>( floor(static_cast<double>(i / ncol)) );
 
     return data(row, column);
 }
 
-//template<class DT>
-//mematrix<DT> mematrix<DT>::operator+(DT toadd)
-//{
-//    mematrix<DT> temp(nrow, ncol);
-//    for (int i = 0; i < nelements; i++)
-//        temp.data[i] = data[i] + toadd;
-//    return temp;
-//}
+// template<class DT>
+// mematrix<DT> mematrix<DT>::operator+(DT toadd)
+// {
+//     mematrix<DT> temp(nrow, ncol);
+//     for (int i = 0; i < nelements; i++)
+//         temp.data[i] = data[i] + toadd;
+//     return temp;
+// }
 
 template<class DT>
 mematrix<DT> mematrix<DT>::operator+(const mematrix<DT> &M)
@@ -318,14 +318,14 @@
     return temp;
 }
 
-//template<class DT>
-//mematrix<double> todouble(mematrix<DT> &M)
-//{
-//    mematrix<double> temp(M.nrow, M.ncol);
-//    for (int i = 0; i < temp.nelements; i++)
-//        temp.data[i] = double(M.data[i]);
-//    return temp;
-//}
+// template<class DT>
+// mematrix<double> todouble(mematrix<DT> &M)
+// {
+//     mematrix<double> temp(M.nrow, M.ncol);
+//     for (int i = 0; i < temp.nelements; i++)
+//         temp.data[i] = double(M.data[i]);
+//     return temp;
+// }
 
 template<class DT>
 mematrix<DT> invert(const mematrix<DT> &M)
@@ -345,14 +345,15 @@
 template<class DT>
 mematrix<DT> productMatrDiag(const mematrix<DT> &M, const mematrix<DT> &D)
 {
-    //multiply all rows of M by value of first row of D
+    // multiply all rows of M by value of first row of D
     if (M.ncol != D.nrow)
     {
         std::cerr << "productMatrDiag: wrong dimensions";
         exit(1);
     }
     mematrix<DT> temp = M;
-    //make a array of the first row of D in the same way orientation as M.data.row(i).array()
+    // make an array of the first row of D in the same way orientation as
+    // M.data.row(i).array()
     Array<DT, Dynamic, Dynamic> row = D.data.block(0, 0, M.ncol, 1).transpose();
 
     for (int i = 0; i < temp.nrow; i++)

Modified: pkg/ProbABEL/src/extract-snp.cpp
===================================================================
--- pkg/ProbABEL/src/extract-snp.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/extract-snp.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -80,7 +80,8 @@
     string snpname = "";
     do
     {
-        next_option = getopt_long(argc, argv, short_options, long_options, NULL);
+        next_option = getopt_long(argc, argv,
+                                  short_options, long_options, NULL);
         switch (next_option)
         {
         case 'h':
@@ -101,8 +102,7 @@
         case '?': break;
         case -1 : break;
         }
-    }
-    while (next_option != -1);
+    } while (next_option != -1);
 
     if (snpname.compare("") == 0)
     {

Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/main.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -144,7 +144,8 @@
 #elif COXPH
     coxph_reg nrd = coxph_reg(nrgd);
     nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
-                 input_var.getInteraction(), input_var.getNgpreds(), true, 1, mli, 0);
+                 input_var.getInteraction(), input_var.getNgpreds(),
+                 true, 1, mli, 0);
 #endif
     double null_loglik = nrd.loglik;
 
@@ -173,7 +174,7 @@
         {
             create_header(outfile, input_var, phd, interaction_cox);
         }
-    } else //Dosage data: Only additive model => only one output file
+    } else  // Dosage data: Only additive model => only one output file
     {
         outfile.push_back(
             new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
@@ -189,7 +190,7 @@
         {
             create_header(outfile, input_var, phd, interaction_cox);
         }
-    } // END else: we have dosage data => only one file
+    }  // END else: we have dosage data => only one file
 
 
     int maxmod = 5;             // Total number of models (in random
@@ -202,9 +203,9 @@
     int start_pos, end_pos;
 
     std::vector<std::ostringstream *> beta_sebeta;
-    //Han Chen
+    // Han Chen
     std::vector<std::ostringstream *> covvalue;
-    //Oct 26, 2009
+    // Oct 26, 2009
     std::vector<std::ostringstream *> chi2;
 
     // Create string streams for betas, SEs, etc. These are used to
@@ -214,15 +215,15 @@
     {
         beta_sebeta.push_back(new std::ostringstream());
         beta_sebeta[i]->precision(6);
-        //*beta_sebeta[i] << scientific;
-        //Han Chen
+        // *beta_sebeta[i] << scientific;
+        // Han Chen
         covvalue.push_back(new std::ostringstream());
         covvalue[i]->precision(6);
-        //*covvalue[i] << scientific;
-        //Oct 26, 2009
+        // *covvalue[i] << scientific;
+        // Oct 26, 2009
         chi2.push_back(new std::ostringstream());
         chi2[i]->precision(6);
-        //*chi2[i] << scientific;
+        // *chi2[i] << scientific;
     }
 
 
@@ -256,7 +257,8 @@
         {
             // Dosage data: only additive model
             int file = 0;
-            write_mlinfo(outfile, file, mli, csnp, input_var, rgd.gcount, rgd.freq);
+            write_mlinfo(outfile, file, mli, csnp, input_var,
+                         rgd.gcount, rgd.freq);
             maxmod = 1;         // We can only calculate the additive
                                 // model with dosage data
         }
@@ -319,7 +321,7 @@
                                         << rd.beta[pos]
                                         << input_var.getSep()
                                         << rd.sebeta[pos];
-                    //Han Chen
+                    // Han Chen
 #if !COXPH
                     if (input_var.getInverseFilename() == NULL
                             && !input_var.getAllcov()
@@ -338,28 +340,26 @@
                                             << input_var.getSep()
                                             << rd.covariance[pos - 2];
                                     }
-
-                                } // END ngpreds=2
+                                }  // END ngpreds=2
                                 else
                                 {
                                     *covvalue[model] << rd.covariance[pos - 1];
                                 }
-
-                            } //END model == 0
+                            }  // END model == 0
                             else
                             {
                                 *covvalue[model] << rd.covariance[pos - 1];
-                            } // END model != 0
-                        } // END if pos > start_pos
+                            }  // END model != 0
+                        }  // END if pos > start_pos
                     }
 #endif
-                    //Oct 26, 2009
-                } // END for(pos = start_pos; pos < rd.beta.nrow; pos++)
+                    // Oct 26, 2009
+                }  // END for(pos = start_pos; pos < rd.beta.nrow; pos++)
 
 
-                //calculate chi^2
-                //________________________________
-                //cout <<  rd.loglik<<" "<<input_var.getNgpreds() << "\n";
+                // calculate chi^2
+                // ________________________________
+                // cout <<  rd.loglik<<" "<<input_var.getNgpreds() << "\n";
 
                 if (input_var.getInverseFilename() == NULL)
                 { // Only if we don't have an inv.sigma file can we use LRT
@@ -420,7 +420,7 @@
                         // We want score test output
                         *chi2[model] << rd.chi2_score;
                     }
-                } // END if( inv.sigma == NULL )
+                }  // END if( inv.sigma == NULL )
                 else if (input_var.getInverseFilename() != NULL)
                 {
                     // We can't use the LRT here, because mmscore is a
@@ -440,10 +440,9 @@
                         *chi2[model] << Z * Z;
                     }
                 }
-            } // END first part of if(poly); allele not too rare
+            }  // END first part of if(poly); allele not too rare
             else
             {   // SNP is rare: beta, sebeta, chi2 = nan
-
                 int number_of_rows_or_columns = rgd.X.ncol;
                 start_pos = get_start_position(input_var, model,
                         number_of_rows_or_columns);
@@ -477,7 +476,7 @@
 
                 if (input_var.getNgpreds() == 2)
                 {
-                    //Han Chen
+                    // Han Chen
 #if !COXPH
                     if (!input_var.getAllcov()
                             && input_var.getInteraction() != 0)
@@ -493,10 +492,10 @@
                         }
                     }
 #endif
-                    //Oct 26, 2009
+                    // Oct 26, 2009
                     *chi2[model] << "nan";
                 } else
-                { //ngpreds==1 (and SNP is rare)
+                { // ngpreds==1 (and SNP is rare)
                     if (input_var.getInverseFilename() == NULL)
                     {
                         //                     Han Chen
@@ -507,18 +506,18 @@
                             *covvalue[model] << "nan";
                         }
 #endif
-                        //Oct 26, 2009
-                    } // END if getInverseFilename == NULL
+                        // Oct 26, 2009
+                    }  // END if getInverseFilename == NULL
                     *chi2[model] << "nan";
-                } // END ngpreds == 1 (and SNP is rare)
-            } // END else: SNP is rare
-        } // END of model cycle
+                }  // END ngpreds == 1 (and SNP is rare)
+            }  // END else: SNP is rare
+        }  // END of model cycle
 
 
         // Start writing beta's, se_beta's etc. to file
         if (input_var.getNgpreds() == 2)
         {
-            for(int model=0; model < maxmod; model++)
+            for (int model = 0; model < maxmod; model++)
             {
                 *outfile[model] << beta_sebeta[model]->str()
                                 << input_var.getSep();
@@ -531,9 +530,9 @@
 #endif
                 *outfile[model] << chi2[model]->str()
                                 << "\n";
-            } // END for loop over all models
+            }  // END for loop over all models
         }
-        else // Dose data: only additive model. Only one output file
+        else  // Dose data: only additive model. Only one output file
         {
             *outfile[0] << beta_sebeta[0]->str() << input_var.getSep();
 #if !COXPH
@@ -550,14 +549,14 @@
         for (int model = 0; model < maxmod; model++)
         {
             beta_sebeta[model]->str("");
-            //Han Chen
+            // Han Chen
             covvalue[model]->str("");
-            //Oct 26, 2009
+            // Oct 26, 2009
             chi2[model]->str("");
         }
 
         update_progress_to_cmd_line(csnp, nsnps);
-    } // END for loop over all SNPs
+    }  // END for loop over all SNPs
 
 
     // We're almost done. All computations have finished, time to

Modified: pkg/ProbABEL/src/main_functions_dump.cpp
===================================================================
--- pkg/ProbABEL/src/main_functions_dump.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/main_functions_dump.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -281,7 +281,8 @@
                         << "sebeta_SNP_A1A1_"
                         << phd.model_terms[interaction_cox];
 #if !COXPH
-            if (input_var.getInverseFilename() == NULL && !input_var.getAllcov())
+            if (input_var.getInverseFilename() == NULL &&
+                !input_var.getAllcov())
             {
                 *outfile[0] << input_var.getSep()
                             << "cov_SNP_A1A2_int_SNP_"
@@ -313,11 +314,11 @@
                 //Oct 26, 2009
             }
         }
-        *outfile[0] << input_var.getSep() << "chi2_SNP_2df\n";  // "loglik\n";
-        *outfile[1] << input_var.getSep() << "chi2_SNP_A1\n";   // "loglik\n";
-        *outfile[2] << input_var.getSep() << "chi2_SNP_domA1\n";// "loglik\n";
-        *outfile[3] << input_var.getSep() << "chi2_SNP_recA1\n";// "loglik\n";
-        *outfile[4] << input_var.getSep() << "chi2_SNP_odomA1\n"; // "loglik\n";
+        *outfile[0] << input_var.getSep() << "chi2_SNP_2df\n";
+        *outfile[1] << input_var.getSep() << "chi2_SNP_A1\n";
+        *outfile[2] << input_var.getSep() << "chi2_SNP_domA1\n";
+        *outfile[3] << input_var.getSep() << "chi2_SNP_recA1\n";
+        *outfile[4] << input_var.getSep() << "chi2_SNP_odomA1\n";
     } // End: ngpreds == 2
     else
     {
@@ -423,9 +424,3 @@
 
     return start_pos;
 }
-
-
-
-
-
-

Modified: pkg/ProbABEL/src/main_functions_dump.h
===================================================================
--- pkg/ProbABEL/src/main_functions_dump.h	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/main_functions_dump.h	2013-12-02 18:22:46 UTC (rev 1444)
@@ -14,17 +14,17 @@
 
 void update_progress_to_cmd_line(int, int);
 void loadInvSigma(cmdvars& input_var, phedata& phd,
-		masked_matrix& invvarmatrix);
+                  masked_matrix& invvarmatrix);
 int create_phenotype(phedata& phd, cmdvars& input_var);
 void create_start_of_header(std::vector<std::ofstream*>& outfile,
-		cmdvars& input_var, phedata& phd);
+                            cmdvars& input_var, phedata& phd);
 void write_mlinfo(const std::vector<std::ofstream*>& outfile, unsigned int file,
-		const mlinfo& mli, int csnp, const cmdvars& input_var, int gcount,
-		double freq);
+                  const mlinfo& mli, int csnp, const cmdvars& input_var,
+                  int gcount, double freq);
 void open_files_for_output(std::vector<std::ofstream*>& outfile,
-		std::string& outfilename_str);
+                           std::string& outfilename_str);
 void create_header(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
-		phedata& phd, int& interaction_cox);
+                   phedata& phd, int& interaction_cox);
 int get_start_position(const cmdvars& input_var, int model,
         int number_of_rows_or_columns);
 #endif /* MAIN_FUNCTIONS_DUMP_H_ */

Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/reg1.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -13,7 +13,7 @@
 // model 3 = recessive 1 df
 // model 4 = over-dominant 1 df
 {
-    if(nullmodel)
+    if (nullmodel)
     {
         // No need to apply any genotypic model when calculating the
         // null model
@@ -385,7 +385,6 @@
         double relative_error = (tXXeigen * betaeigen - tXYeigen).norm() /
             tXYeigen.norm(); // norm() is L2 norm
         cout << "The relative error is:\n" << relative_error << endl;
-
     }
 
     // This one is needed later on in this function

Modified: pkg/ProbABEL/src/regdata.cpp
===================================================================
--- pkg/ProbABEL/src/regdata.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/regdata.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -25,8 +25,8 @@
     noutcomes               = 0;
     is_interaction_excluded = false;
     masked_data             = NULL;
-    gcount					=0;
-    freq					=0;
+    gcount                  = 0;
+    freq                    = 0;
 }
 
 
@@ -45,6 +45,7 @@
     }
 }
 
+
 regdata::regdata(phedata &phed, gendata &gend, int snpnum,
                  bool ext_is_interaction_excluded)
 {
@@ -98,11 +99,13 @@
     is_interaction_excluded = ext_is_interaction_excluded;
 }
 
+
 void regdata::update_snp(gendata &gend, int snpnum)
 {
-	//reset counter for frequency since it is a new snp
-	gcount=0;
-	freq=0.0;
+    // Reset counter for frequency since it is a new SNP
+    gcount = 0;
+    freq = 0.0;
+
     // Add genotypic data (dosage or probabilities) to the design
     // matrix X
     for (int j = 0; j < ngpreds; j++)
@@ -115,45 +118,46 @@
 
         gend.get_var(snpnum * ngpreds + j, snpdata);
 
-		for (int i = 0; i < nids; i++) {
-			X.put(snpdata[i], i, (ncov - j));
-			if (std::isnan(snpdata[i])) {
-				masked_data[i] = 1;
-				//snp not masked
-			} else {
-				// checck for first predicor
-				if (j == 0) {
-					gcount++;
-					if (ngpreds == 1) {
-						freq += snpdata[i] * 0.5;
-					} else if (ngpreds == 2) {
-						freq += snpdata[i];
-					}
-				} else if (j == 1) {
-					// add second genotype in two predicor data form
-					freq += snpdata[i] * 0.5;
-				}
-			}//end std::isnan(snpdata[i]) snp
+        for (int i = 0; i < nids; i++) {
+            X.put(snpdata[i], i, (ncov - j));
+            if (std::isnan(snpdata[i])) {
+                masked_data[i] = 1;
+                // SNP not masked
+            } else {
+                // check for first predictor
+                if (j == 0) {
+                    gcount++;
+                    if (ngpreds == 1) {
+                        freq += snpdata[i] * 0.5;
+                    } else if (ngpreds == 2) {
+                        freq += snpdata[i];
+                    }
+                } else if (j == 1) {
+                    // Add second genotype in two predicor data form
+                    freq += snpdata[i] * 0.5;
+                }
+            }  // End std::isnan(snpdata[i]) snp
+        }  // End i for loop
 
-		}//end i for loop
+        delete[] snpdata;
+    }  // End ngpreds loop
 
+    freq /= static_cast<double>(gcount); // Allele frequency
+}
 
-		  delete[] snpdata;
-}//end ngpreds loop
-	freq /= static_cast<double>(gcount); // Allele frequency
 
-}
+/**
+ * update_snp() adds SNP information to the design matrix. This
+ * function allows you to strip that information from X again.
+ * This is used for example when calculating the null model.
+ */
 void regdata::remove_snp_from_X()
 {
-    // update_snp() adds SNP information to the design matrix. This
-    // function allows you to strip that information from X again.
-    // This is used for example when calculating the null model.
-
-    if(ngpreds == 1)
+    if (ngpreds == 1)
     {
         X.delete_column(X.ncol -1);
     }
-    else if(ngpreds == 2)
+    else if (ngpreds == 2)
     {
         X.delete_column(X.ncol -1);
         X.delete_column(X.ncol -1);
@@ -165,16 +169,18 @@
     }
 }
 
+
 regdata::~regdata()
 {
     delete[] regdata::masked_data;
-    //      delete X;
-    //      delete Y;
+    // delete X;
+    // delete Y;
 }
 
+
 regdata regdata::get_unmasked_data()
 {
-    regdata to; // = regdata(*this);
+    regdata to;  // = regdata(*this);
     int nmeasured = 0;
     for (int i = 0; i < nids; i++)
     {
@@ -210,7 +216,7 @@
         }
     }
 
-    //delete [] to.masked_data;
+    // delete [] to.masked_data;
     to.masked_data = new unsigned short int[to.nids];
     for (int i = 0; i < to.nids; i++)
     {
@@ -221,6 +227,7 @@
     return (to);
 }
 
+
 mematrix<double> regdata::extract_genotypes(void)
 {
     mematrix<double> out;

Modified: pkg/ProbABEL/src/usage.cpp
===================================================================
--- pkg/ProbABEL/src/usage.cpp	2013-12-02 17:39:13 UTC (rev 1443)
+++ pkg/ProbABEL/src/usage.cpp	2013-12-02 18:22:46 UTC (rev 1444)
@@ -32,8 +32,9 @@
     cout << "\t --out     : [optional] output file name "
          << "(default is regression.out.txt)"
          << endl;
-    cout << "\t --skipd   : [optional] how many columns to skip in the predictor"
-         << "\n\t             (dose/prob) file (default 2)"
+    cout << "\t --skipd   : [optional] how many columns to skip in the "
+         << "predictor\n"
+         << "\t             (dose/prob) file (default 2)"
          << endl;
 #if COXPH
     cout << "\t --ntraits : [optional] how many traits are analysed (default 2)"
@@ -42,15 +43,16 @@
     cout << "\t --ntraits : [optional] how many traits are analysed (default 1)"
          << endl;
 #endif
-    cout << "\t --ngpreds : [optional] how many predictor columns per marker"
-         <<"\n\t              (default 1 = MLDOSE; else use 2 for MLPROB)"
+    cout << "\t --ngpreds : [optional] how many predictor columns per marker\n"
+         <<"\t              (default 1 = MLDOSE; else use 2 for MLPROB)"
          << endl;
     cout << "\t --separat : [optional] character to separate fields "
          << "(default is space)"
          << endl;
     cout << "\t --score   : use score test" << endl;
     cout << "\t --no-head : do not report header line" << endl;
-    cout << "\t --allcov  : report estimates for all covariates (large outputs!)"
+    cout << "\t --allcov  : report estimates for all covariates "
+         << "(large outputs!)"
          << endl;
     cout << "\t --interaction: Which covariate to use for interaction with "
          << "SNP analysis (default is no interaction, 0)"
@@ -72,7 +74,7 @@
 }
 
 
-void print_version(void){
+void print_version(void) {
     cout << PACKAGE
          << " v. " << PACKAGE_VERSION
          << "\n(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, "



More information about the Genabel-commits mailing list