[Genabel-commits] r1033 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 28 13:34:05 CET 2012
Author: lckarssen
Date: 2012-11-28 13:34:05 +0100 (Wed, 28 Nov 2012)
New Revision: 1033
Modified:
pkg/ProbABEL/src/cholesky.cpp
pkg/ProbABEL/src/main.cpp
pkg/ProbABEL/src/phedata.cpp
pkg/ProbABEL/src/phedata.h
pkg/ProbABEL/src/regdata.cpp
pkg/ProbABEL/src/usage.cpp
Log:
ProbABEL: some more code beautification. No functional changes. Copyright for these changes lies with the Erasmus Medical Centre.
Modified: pkg/ProbABEL/src/cholesky.cpp
===================================================================
--- pkg/ProbABEL/src/cholesky.cpp 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/cholesky.cpp 2012-11-28 12:34:05 UTC (rev 1033)
@@ -4,6 +4,11 @@
* Created on: Mar 15, 2012
* Author: mkooyman
*/
+#include <string>
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+
#if EIGEN
#include "eigen_mematrix.h"
#include "eigen_mematrix.cpp"
@@ -12,10 +17,6 @@
#include "mematri1.h"
#endif
-#include <string>
-#include <cstdarg>
-#include <cstdio>
-#include <cstdlib>
/* SCCS @(#)cholesky2.c 5.2 10/27/98
** subroutine to do Cholesky decompostion on a matrix: C = FDF'
@@ -73,7 +74,8 @@
matrix[i * n + i] = 0;
if (pivot < -8 * eps)
nonneg = -1;
- } else
+ }
+ else
{
rank++;
for (j = (i + 1); j < n; j++)
@@ -132,7 +134,8 @@
matrix[j * n + i] = 0;
for (j = i; j < n; j++)
matrix[i * n + j] = 0;
- } else
+ }
+ else
{
for (j = (i + 1); j < n; j++)
{
@@ -149,4 +152,3 @@
for (int row = 0; row < col; row++)
matrix[col * n + row] = matrix[row * n + col];
}
-
Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/main.cpp 2012-11-28 12:34:05 UTC (rev 1033)
@@ -59,10 +59,11 @@
if (csnp == 0)
{
std::cout << "Analysis: "
- << setw (5) << 100. * static_cast<double>(csnp)
+ << setw(5) << 100. * static_cast<double>(csnp)
/ static_cast<double>(nsnps)
<< "%...";
- } else
+ }
+ else
{
std::cout << "\b\b\b\b\b\b\b\b\b"
<< setw (5) << 100. * static_cast<double>(csnp)
@@ -116,8 +117,8 @@
input_var.getInteraction() > phd.ncov ||
interaction_cox > phd.ncov)
{
- std::cerr << "error: Interaction parameter is out of range (interaction="
- << input_var.getInteraction() << ") \n";
+ std::cerr << "error: Interaction parameter is out of range "
+ << "(interaction=" << input_var.getInteraction() << ") \n";
exit(1);
}
@@ -232,19 +233,20 @@
phedata phd, int interaction_cox)
{
create_start_of_header(outfile, input_var, phd);
- *outfile[0] << input_var.getSep() << "beta_SNP_add" << input_var.getSep()
- << "sebeta_SNP_add";
+ *outfile[0] << input_var.getSep() << "beta_SNP_add"
+ << input_var.getSep() << "sebeta_SNP_add";
if (input_var.getInteraction() != 0)
{
*outfile[0] << input_var.getSep() << "beta_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_" << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "sebeta_SNP_"
+ << phd.model_terms[interaction_cox];
}
if (input_var.getInverseFilename() == NULL)
- //Han Chen
{
+ //Han Chen
#if !COXPH
if (input_var.getInteraction() != 0 && !input_var.getAllcov())
{
@@ -266,7 +268,8 @@
input_var.printinfo();
// if (allcov && ngpreds>1)
// {
- // std::cout << "\n\nWARNING: --allcov allowed only for 1 predictor (MLDOSE)\n";
+ // cout << "\n\n"
+ // << "WARNING: --allcov allowed only for 1 predictor (MLDOSE)\n";
// allcov = 0;
// }
mlinfo mli(input_var.getMlinfofilename(), input_var.getMapfilename());
@@ -275,21 +278,32 @@
int interaction_cox = create_phenotype(phd, input_var);
//interaction--;
- // if(input_var.getInverseFilename()!= NULL && phd.ncov > 1)
- // {
- // std::cerr<<"Error: In mmscore you can not use any covariates. You phenotype file must conatin id column and trait (residuals) only\n";
- // exit(1);
- // }
- // if(input_var.getInverseFilename()!= NULL && (allcov == 1 || score == 1 || input_var.getInteraction()!= 0 || ngpreds==2))
- // {
- // std::cerr<<"Error: In mmscore you can use additive model without any inetractions only\n";
- // exit(1);
- // }
+ // if (input_var.getInverseFilename()!= NULL && phd.ncov > 1)
+ // {
+ // std::cerr << "Error: In mmscore you can not use any covariates."
+ // << " You phenotype file must conatin id column and "
+ // << "trait (residuals) only\n";
+ // exit(1);
+ // }
+ // if (input_var.getInverseFilename()!= NULL &&
+ // (allcov == 1 || score == 1
+ // || input_var.getInteraction()!= 0
+ // || ngpreds==2))
+ // {
+ // std::cerr << "Error: In mmscore you can use additive model "
+ // << "without any inetractions only\n";
+ // exit(1);
+ // }
masked_matrix invvarmatrix;
+
/*
* now should be possible... delete this part later when everything works
#if LOGISTIC
- if(input_var.getInverseFilename()!= NULL) {std::cerr<<"ERROR: mmscore is forbidden for logistic regression\n";exit(1);}
+ if (input_var.getInverseFilename()!= NULL)
+ {
+ std::cerr << "ERROR: mmscore is forbidden for logistic regression\n";
+ exit(1);
+ }
#endif
*/
@@ -297,7 +311,6 @@
if (input_var.getInverseFilename() != NULL)
{
loadInvSigma(input_var, phd, invvarmatrix);
- // matrix.print();
}
gendata gtd;
@@ -316,10 +329,14 @@
std::cout << " loaded genotypic data ..." << std::flush;
/**
if (input_var.getIsFvf())
- gendata gtd (str_genfilename,nsnps,input_var.getNgpreds(),phd.nids_all,phd.allmeasured,phd.idnames);
- else
- gendata gtd (input_var.getGenfilename(),nsnps,input_var.getNgpreds(),phd.nids_all,phd.nids,phd.allmeasured,skipd,phd.idnames);
+ gendata gtd(str_genfilename, nsnps, input_var.getNgpreds(),
+ phd.nids_all, phd.allmeasured, phd.idnames);
+ else
+ gendata gtd(input_var.getGenfilename(), nsnps,
+ input_var.getNgpreds(), phd.nids_all, phd.nids,
+ phd.allmeasured, skipd, phd.idnames);
**/
+
// estimate null model
#if COXPH
coxph_data nrgd = coxph_data(phd, gtd, -1);
@@ -357,14 +374,13 @@
#endif
std::cout << " formed regression object ...";
-
std::cout << " done\n" << std::flush;
//________________________________________________________________
//Maksim, 9 Jan, 2009
-
std::string outfilename_str(input_var.getOutfilename());
std::vector<std::ofstream*> outfile;
+
//All models output.One file per each model
if (input_var.getNgpreds() == 2)
{
@@ -373,7 +389,8 @@
{
create_header_1(outfile, input_var, phd, interaction_cox);
}
- } else //Only additive model. Only one output file
+ }
+ else //Only additive model. Only one output file
{
outfile.push_back(
new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
@@ -409,19 +426,32 @@
}
if (input_var.getNgpreds()==2)
{
- outfile << input_var.getSep() << "beta_SNP_A1A2" << input_var.getSep() << "beta_SNP_A1A1" << input_var.getSep()
- << "sebeta_SNP_A1A2" << input_var.getSep() << "sebeta_SNP_a1A1" << input_var.getSep() << "chi2_SNP_2df"
- << input_var.getSep() << "beta_SNP_addA1" << input_var.getSep() << "sebeta_SNP_addA1" << input_var.getSep() << "chi2_SNP_addA1"
- << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep() << "sebeta_SNP_domA1" << input_var.getSep() << "chi2_SNP_domA1"
- << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep() << "sebeta_SNP_recA1" << input_var.getSep() << "chi2_SNP_recA1"
- << input_var.getSep() << "beta_SNP_odom" << input_var.getSep() << "sebeta_SNP_odom" << input_var.getSep() << "chi2_SNP_odom\n";
+ outfile << input_var.getSep() << "beta_SNP_A1A2"
+ << input_var.getSep() << "beta_SNP_A1A1"
+ << input_var.getSep() << "sebeta_SNP_A1A2"
+ << input_var.getSep() << "sebeta_SNP_a1A1"
+ << input_var.getSep() << "chi2_SNP_2df"
+ << input_var.getSep() << "beta_SNP_addA1"
+ << input_var.getSep() << "sebeta_SNP_addA1"
+ << input_var.getSep() << "chi2_SNP_addA1"
+ << input_var.getSep() << "beta_SNP_domA1"
+ << input_var.getSep() << "sebeta_SNP_domA1"
+ << input_var.getSep() << "chi2_SNP_domA1"
+ << input_var.getSep() << "beta_SNP_recA1"
+ << input_var.getSep() << "sebeta_SNP_recA1"
+ << input_var.getSep() << "chi2_SNP_recA1"
+ << input_var.getSep() << "beta_SNP_odom"
+ << input_var.getSep() << "sebeta_SNP_odom"
+ << input_var.getSep() << "chi2_SNP_odom\n";
}
else
{
- outfile << input_var.getSep() << "beta_SNP_add" << input_var.getSep() << "sebeta_SNP_add" << input_var.getSep() << "chi2_SNP_add\n";
+ outfile << input_var.getSep() << "beta_SNP_add"
+ << input_var.getSep() << "sebeta_SNP_add"
+ << input_var.getSep() << "chi2_SNP_add\n";
}
- }
- */
+
+ */
// exit(1);
//________________________________________________________________
//Maksim, 9 Jan, 2009
@@ -462,7 +492,8 @@
gcount++;
freq += snpdata1[ii] + snpdata2[ii] * 0.5;
}
- } else
+ }
+ else
{
// freq = (gtd.G).column_mean(csnp)/2.;
gtd.get_var(csnp, snpdata1);
@@ -579,7 +610,8 @@
<< input_var.getSep()
<< rd.covariance[pos - 2];
}
- } else
+ }
+ else
{
*covvalue[model] << rd.covariance[pos - 1];
}
@@ -595,13 +627,15 @@
{
//*chi2[model] << 2.*(rd.loglik-null_loglik);
*chi2[model] << rd.loglik;
- } else
+ }
+ else
{
//*chi2[model] << rd.chi2_score;
*chi2[model] << "nan";
}
//________________________________
- } else //beta, sebeta = nan
+ }
+ else //beta, sebeta = nan
{
if (!input_var.getAllcov() && model == 0
&& input_var.getInteraction() == 0)
@@ -621,7 +655,8 @@
if (model == 0)
{
end_pos = rgd.X.ncol;
- } else
+ }
+ else
{
end_pos = rgd.X.ncol - 1;
}
@@ -643,7 +678,8 @@
{
*covvalue[model] << "nan"
<< input_var.getSep() << "nan";
- } else
+ }
+ else
{
*covvalue[model] << "nan";
}
@@ -700,7 +736,8 @@
#endif
*outfile[4] << chi2[4]->str() << "\n";
//Oct 26, 2009
- } else //Only additive model. Only one output file
+ }
+ else //Only additive model. Only one output file
{
//Write mlinfo to output:
*outfile[0] << mli.name[csnp]
@@ -771,7 +808,8 @@
cout << rd.loglik << " <-logliken\n";
cout << rd.sigma2 << " <-sigma2\n";
#endif
- } else
+ }
+ else
{
// if(input_var.getInverseFilename()== NULL)
// {
@@ -815,8 +853,9 @@
//}
//else
//{
- //rd.mmscore(rgd,0,CHOLTOL,model, input_var.getInteraction(),
- //input_var.getNgpreds(), invvarmatrix);
+ // rd.mmscore(rgd, 0, CHOLTOL, model,
+ // input_var.getInteraction(),
+ // input_var.getNgpreds(), invvarmatrix);
//}
}
#elif COXPH
@@ -829,11 +868,13 @@
if (!input_var.getAllcov() && input_var.getInteraction() == 0)
{
start_pos = rd.beta.nrow - 1;
- } else if (!input_var.getAllcov()
- && input_var.getInteraction() != 0)
+ }
+ else if (!input_var.getAllcov()
+ && input_var.getInteraction() != 0)
{
start_pos = rd.beta.nrow - 2;
- } else
+ }
+ else
{
start_pos = 0;
}
@@ -869,13 +910,15 @@
if (input_var.getScore() == 0)
{
*chi2[0] << rd.loglik; //2.*(rd.loglik-null_loglik);
- } else
+ }
+ else
{
*chi2[0] << "nan"; //rd.chi2_score;
}
}
//________________________________
- } else //beta, sebeta = nan
+ }
+ else //beta, sebeta = nan
{
if (!input_var.getAllcov() && input_var.getInteraction() == 0)
start_pos = rgd.X.ncol - 1;
@@ -927,7 +970,8 @@
#endif
*outfile[0] << chi2[model]->str() << "\n";
//Oct 26, 2009
- } else
+ }
+ else
{
*outfile[0] << beta_sebeta[0]->str() << "\n";
#if DEBUG
Modified: pkg/ProbABEL/src/phedata.cpp
===================================================================
--- pkg/ProbABEL/src/phedata.cpp 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/phedata.cpp 2012-11-28 12:34:05 UTC (rev 1033)
@@ -36,7 +36,6 @@
// std::cout << line << "\n ";
while (line_stream >> tmp)
{
-
nphenocols++;
// std::cout << tmp << " " << nphenocols << " ";
}
@@ -180,7 +179,9 @@
// allocate objects
int ntmpcov = 1;
if (ncov > 0)
+ {
ntmpcov = ncov;
+ }
idnames = new std::string[nids];
X.reinit(nids, ntmpcov);
@@ -219,10 +220,13 @@
m++;
}
else
+ {
for (int j = 0; j < nphenocols; j++)
infile >> tmp;
+ }
infile.close();
}
+
phedata::~phedata()
{
// delete X;
Modified: pkg/ProbABEL/src/phedata.h
===================================================================
--- pkg/ProbABEL/src/phedata.h 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/phedata.h 2012-11-28 12:34:05 UTC (rev 1033)
@@ -17,7 +17,6 @@
#endif
class phedata {
-
public:
phedata()
{
@@ -54,4 +53,5 @@
~phedata();
};
+
#endif /* PHEDATA_H_ */
Modified: pkg/ProbABEL/src/regdata.cpp
===================================================================
--- pkg/ProbABEL/src/regdata.cpp 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/regdata.cpp 2012-11-28 12:34:05 UTC (rev 1033)
@@ -19,16 +19,15 @@
regdata::regdata()
{
- nids = 0;
- ncov = 0;
- ngpreds = 0;
- noutcomes = 0;
+ nids = 0;
+ ncov = 0;
+ ngpreds = 0;
+ noutcomes = 0;
is_interaction_excluded = false;
- masked_data = NULL;
-
+ masked_data = NULL;
}
-;
+
regdata::regdata(const regdata &obj) : X(obj.X), Y(obj.Y)
{
nids = obj.nids;
@@ -37,13 +36,15 @@
noutcomes = obj.noutcomes;
is_interaction_excluded = obj.is_interaction_excluded;
masked_data = new unsigned short int[nids];
+
for (int i = 0; i < nids; i++)
{
masked_data[i] = 0;
}
}
+
regdata::regdata(phedata &phed, gendata &gend, int snpnum,
- bool ext_is_interaction_excluded)
+ bool ext_is_interaction_excluded)
{
nids = gend.nids;
masked_data = new unsigned short int[nids];
@@ -69,9 +70,15 @@
X.put(1., i, 0);
Y.put((phed.Y).get(i, 0), i, 0);
}
+
for (int j = 1; j <= phed.ncov; j++)
+ {
for (int i = 0; i < nids; i++)
+ {
X.put((phed.X).get(i, j - 1), i, j);
+ }
+ }
+
if (snpnum > 0)
for (int j = 0; j < ngpreds; j++)
{
@@ -80,19 +87,21 @@
for (int i = 0; i < nids; i++)
X.put(snpdata[i], i, (ncov - ngpreds + 1 + j));
}
- // for (int i=0;i<nids;i++)
- // for (int j=0;j<ngpreds;j++)
- // X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+1+j));
+ // for (int i=0;i<nids;i++)
+ // for (int j=0;j<ngpreds;j++)
+ // X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+1+j));
is_interaction_excluded = ext_is_interaction_excluded;
+}
-}
void regdata::update_snp(gendata &gend, int snpnum)
{
for (int j = 0; j < ngpreds; j++)
{
float snpdata[nids];
for (int i = 0; i < nids; i++)
+ {
masked_data[i] = 0;
+ }
gend.get_var(snpnum * ngpreds + j, snpdata);
@@ -104,6 +113,7 @@
}
}
}
+
regdata::~regdata()
{
delete[] regdata::masked_data;
@@ -116,13 +126,15 @@
regdata to; // = regdata(*this);
int nmeasured = 0;
for (int i = 0; i < nids; i++)
+ {
if (masked_data[i] == 0)
nmeasured++;
- to.nids = nmeasured;
- //cout << to.nids << " in get_unmasked_data\n";
- to.ncov = ncov;
- to.ngpreds = ngpreds;
- to.noutcomes = noutcomes;
+ }
+
+ to.nids = nmeasured;
+ to.ncov = ncov;
+ to.ngpreds = ngpreds;
+ to.noutcomes = noutcomes;
to.is_interaction_excluded = is_interaction_excluded;
int dim2Y = Y.ncol;
int dim2X = X.ncol;
@@ -150,7 +162,9 @@
//delete [] to.masked_data;
to.masked_data = new unsigned short int[to.nids];
for (int i = 0; i < to.nids; i++)
+ {
to.masked_data[i] = 0;
+ }
// std::cout << "get_unmasked: " << to.nids << " "
// << dim2X << " " << dim2Y << "\n";
return (to);
@@ -161,7 +175,11 @@
mematrix<double> out;
out.reinit(X.nrow, ngpreds);
for (int i = 0; i < X.nrow; i++)
+ {
for (int j = 0; j < ngpreds; j++)
+ {
out[i * ngpreds + j] = X.get(i, (ncov - ngpreds + 1 + j));
+ }
+ }
return out;
}
Modified: pkg/ProbABEL/src/usage.cpp
===================================================================
--- pkg/ProbABEL/src/usage.cpp 2012-11-28 10:58:40 UTC (rev 1032)
+++ pkg/ProbABEL/src/usage.cpp 2012-11-28 12:34:05 UTC (rev 1033)
@@ -30,7 +30,8 @@
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)"
More information about the Genabel-commits
mailing list