[Genabel-commits] r981 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 23 18:52:51 CEST 2012
Author: maartenk
Date: 2012-10-23 18:52:51 +0200 (Tue, 23 Oct 2012)
New Revision: 981
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h
branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp
branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
branches/ProbABEL-refactoring/ProbABEL/src/maskedmatrix.cpp
branches/ProbABEL-refactoring/ProbABEL/src/maskedmatrix.h
branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
Log:
-added const to some functions : e.g foo(int &bar) to foo(const int &bar)
-removed extra whitespace
Modified: branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp 2012-10-23 16:52:51 UTC (rev 981)
@@ -21,7 +21,6 @@
{
return str_genfilename;
}
-;
int cmdvars::getAllcov() const
{
@@ -336,20 +335,19 @@
#if COXPH
if (score)
{
- fprintf(stderr,"\n\nOption --score is implemented for linear and logistic models only\n");
+ fprintf(stderr, "\n\nOption --score is implemented for linear and logistic models only\n");
exit(1);
}
- if(inverse_filename != NULL)
+ if (inverse_filename != NULL)
{
- std::cerr<<"ERROR: mmscore is forbidden for cox regression\n";
+ std::cerr << "ERROR: mmscore is forbidden for cox regression\n";
exit(1);
}
if (robust)
{
- std::cerr<<"ERROR: robust standard errors not implemented for Cox regression\n";
+ std::cerr << "ERROR: robust standard errors not implemented for Cox regression\n";
exit(1);
}
#endif
-
}
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-10-23 16:52:51 UTC (rev 981)
@@ -38,7 +38,7 @@
return -9;
}
-coxph_data::coxph_data(const coxph_data &obj)
+coxph_data::coxph_data( coxph_data &obj)
{
nids = obj.nids;
ncov = obj.ncov;
@@ -54,7 +54,7 @@
for (int i = 0; i < nids; i++)
masked_data[i] = 0;
}
-coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum)
+coxph_data::coxph_data( phedata &phed, gendata &gend, int snpnum)
{
nids = gend.nids;
masked_data = new unsigned short int[nids];
@@ -156,7 +156,7 @@
// stime.print();
// sstat.print();
}
-void coxph_data::update_snp(gendata &gend, int snpnum)
+void coxph_data::update_snp( gendata &gend, int snpnum)
{
// note this sorts by "order"!!!
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-10-23 16:52:51 UTC (rev 981)
@@ -37,9 +37,9 @@
coxph_data()
{
}
- coxph_data(const coxph_data &obj);
- coxph_data(phedata &phed, gendata &gend, int snpnum);
- void update_snp(gendata &gend, int snpnum);
+ coxph_data( coxph_data &obj);
+ coxph_data( phedata &phed, gendata &gend, int snpnum);
+ void update_snp( gendata &gend, int snpnum);
~coxph_data();
};
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp 2012-10-23 16:52:51 UTC (rev 981)
@@ -6,7 +6,7 @@
#include "survproto.h"
}
-coxph_reg::coxph_reg(coxph_data &cdatain) {
+coxph_reg::coxph_reg(const coxph_data &cdatain) {
coxph_data cdata = cdatain.get_unmasked_data();
beta.reinit(cdata.X.nrow, 1);
sebeta.reinit(cdata.X.nrow, 1);
@@ -19,7 +19,7 @@
// delete beta;
// delete sebeta;
}
-void coxph_reg::estimate(coxph_data &cdatain, int verbose, int maxiter,
+void coxph_reg::estimate( coxph_data &cdatain, int verbose, int maxiter,
double eps, double tol_chol, int model, int interaction, int ngpreds,
bool iscox, int nullmodel = 0) {
// cout << "model = " << model << "\n";
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h 2012-10-23 16:52:51 UTC (rev 981)
@@ -25,9 +25,9 @@
double chi2_score;
int niter;
- coxph_reg(coxph_data &cdatain);
+ coxph_reg(const coxph_data &cdatain);
virtual ~coxph_reg();
- void coxph_reg::estimate(coxph_data &cdatain, int verbose, int maxiter,
+ void coxph_reg::estimate( coxph_data &cdatain, int verbose, int maxiter,
double eps, double tol_chol, int model, int interaction, int ngpreds,
bool iscox, int nullmodel = 0);
};
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-23 15:42:20 UTC (rev 980)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-23 16:52:51 UTC (rev 981)
@@ -51,907 +51,816 @@
#define EPS 1.e-8
#define CHOLTOL 1.5e-12
-void update_progress_to_cmd_line(int csnp, int nsnps)
-{
- if (csnp % 1000 == 0)
- {
- if (csnp == 0)
- {
- fprintf(stdout, "Analysis: %6.2f ...",
- 100. * static_cast<double>(csnp) / static_cast<double>(nsnps));
- }
- else
- {
- fprintf(stdout, "\b\b\b\b\b\b\b\b\b\b%6.2f ...",
- 100. * static_cast<double>(csnp) / static_cast<double>(nsnps));
- }
- std::cout.flush();
- }
+void update_progress_to_cmd_line(int csnp, int nsnps) {
+ if (csnp % 1000 == 0) {
+ if (csnp == 0) {
+ fprintf(stdout, "Analysis: %6.2f ...",
+ 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps));
+ } else {
+ fprintf(stdout, "\b\b\b\b\b\b\b\b\b\b%6.2f ...",
+ 100. * static_cast<double>(csnp)
+ / static_cast<double>(nsnps));
+ }
+ std::cout.flush();
+ }
}
void open_files_for_output(std::vector<std::ofstream*>& outfile,
- std::string& outfilename_str)
-{
- //create a list of filenames
- const int amount_of_files = 5;
- std::string filenames[amount_of_files] =
- { outfilename_str + "_2df.out.txt", outfilename_str
- + "_add.out.txt", outfilename_str + "_domin.out.txt",
- outfilename_str + "_recess.out.txt", outfilename_str
- + "_over_domin.out.txt" };
+ std::string& outfilename_str) {
+ //create a list of filenames
+ const int amount_of_files = 5;
+ std::string filenames[amount_of_files] = { outfilename_str + "_2df.out.txt",
+ outfilename_str + "_add.out.txt", outfilename_str
+ + "_domin.out.txt", outfilename_str + "_recess.out.txt",
+ outfilename_str + "_over_domin.out.txt" };
- for (int i = 0; i < amount_of_files; i++)
- {
- outfile.push_back(new std::ofstream());
- outfile[i]->open((filenames[i]).c_str());
- if (!outfile[i]->is_open())
- {
- std::cerr << "Can not open file for writing: " << filenames[i]
- << "\n";
- exit(1);
- }
- }
+ for (int i = 0; i < amount_of_files; i++) {
+ outfile.push_back(new std::ofstream());
+ outfile[i]->open((filenames[i]).c_str());
+ if (!outfile[i]->is_open()) {
+ std::cerr << "Can not open file for writing: " << filenames[i]
+ << "\n";
+ exit(1);
+ }
+ }
}
-int create_phenoytype(phedata& phd, cmdvars& input_var)
-{
- phd.set_is_interaction_excluded(input_var.isIsInteractionExcluded());
- phd.setphedata(input_var.getPhefilename(), input_var.getNoutcomes(),
- input_var.getNpeople(), input_var.getInteraction(),
- input_var.isIscox());
- int interaction_cox = input_var.getInteraction();
+int create_phenoytype(phedata& phd, cmdvars& input_var) {
+ phd.set_is_interaction_excluded(input_var.isIsInteractionExcluded());
+ phd.setphedata(input_var.getPhefilename(), input_var.getNoutcomes(),
+ input_var.getNpeople(), input_var.getInteraction(),
+ input_var.isIscox());
+ int interaction_cox = input_var.getInteraction();
#if COXPH
- interaction_cox--;
+ interaction_cox--;
#endif
- if (input_var.getInteraction() < 0 || input_var.getInteraction() > phd.ncov
- || interaction_cox > phd.ncov)
- {
- std::cerr << "error: Interaction parameter is out of range (ineraction="
- << input_var.getInteraction() << ") \n";
- exit(1);
- }
+ if (input_var.getInteraction() < 0 || input_var.getInteraction() > phd.ncov
+ || interaction_cox > phd.ncov) {
+ std::cerr << "error: Interaction parameter is out of range (ineraction="
+ << input_var.getInteraction() << ") \n";
+ exit(1);
+ }
- return interaction_cox;
+ return interaction_cox;
}
void loadInvSigma(cmdvars& input_var, phedata& phd,
- masked_matrix& invvarmatrix)
-{
- std::cout << "you are running mmscore...\n";
- InvSigma inv(input_var.getInverseFilename(), &phd);
- // invvarmatrix = inv.get_matrix();
- //double par = 1.; //var(phd.Y)*phd.nids/(phd.nids-phd.ncov-1);
- invvarmatrix.set_matrix(inv.get_matrix());// = invvarmatrix * par;
- std::cout << " loaded InvSigma ...";
+ masked_matrix& invvarmatrix) {
+ std::cout << "you are running mmscore...\n";
+ InvSigma inv(input_var.getInverseFilename(), &phd);
+ // invvarmatrix = inv.get_matrix();
+ //double par = 1.; //var(phd.Y)*phd.nids/(phd.nids-phd.ncov-1);
+ invvarmatrix.set_matrix(inv.get_matrix()); // = invvarmatrix * par;
+ std::cout << " loaded InvSigma ...";
}
void create_start_of_header(std::vector<std::ofstream*>& outfile,
- cmdvars& input_var, phedata& phd)
-{
- for (unsigned int i = 0; i < outfile.size(); i++)
- {
- (*outfile[i]) << "name" << input_var.getSep() << "A1"
- << input_var.getSep() << "A2" << input_var.getSep() << "Freq1"
- << input_var.getSep() << "MAF" << input_var.getSep()
- << "Quality" << input_var.getSep() << "Rsq"
- << input_var.getSep() << "n" << input_var.getSep()
- << "Mean_predictor_allele";
- if (input_var.getChrom() != "-1")
- (*outfile[i]) << input_var.getSep() << "chrom";
- if (input_var.getMapfilename() != NULL)
- (*outfile[i]) << input_var.getSep() << "position";
- }
+ cmdvars& input_var, phedata& phd) {
+ for (unsigned int i = 0; i < outfile.size(); i++) {
+ (*outfile[i]) << "name" << input_var.getSep() << "A1"
+ << input_var.getSep() << "A2" << input_var.getSep() << "Freq1"
+ << input_var.getSep() << "MAF" << input_var.getSep()
+ << "Quality" << input_var.getSep() << "Rsq"
+ << input_var.getSep() << "n" << input_var.getSep()
+ << "Mean_predictor_allele";
+ if (input_var.getChrom() != "-1")
+ (*outfile[i]) << input_var.getSep() << "chrom";
+ if (input_var.getMapfilename() != NULL)
+ (*outfile[i]) << input_var.getSep() << "position";
+ }
- if (input_var.getAllcov()) //All covariates in output
- {
- for (unsigned int file = 0; file < outfile.size(); file++)
- for (int i = 0; i < phd.n_model_terms - 1; i++)
- *outfile[file] << input_var.getSep() << "beta_"
- << phd.model_terms[i] << input_var.getSep() << "sebeta_"
- << phd.model_terms[i];
- }
+ if (input_var.getAllcov()) //All covariates in output
+ {
+ for (unsigned int file = 0; file < outfile.size(); file++)
+ for (int i = 0; i < phd.n_model_terms - 1; i++)
+ *outfile[file] << input_var.getSep() << "beta_"
+ << phd.model_terms[i] << input_var.getSep() << "sebeta_"
+ << phd.model_terms[i];
+ }
}
void create_header_1(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- phedata& phd, int& interaction_cox)
-{
- create_start_of_header(outfile, input_var, phd);
+ phedata& phd, int& interaction_cox) {
+ create_start_of_header(outfile, input_var, phd);
- *outfile[0] << input_var.getSep() << "beta_SNP_A1A2" << input_var.getSep()
- << "beta_SNP_A1A1" << input_var.getSep() << "sebeta_SNP_A1A2"
- << input_var.getSep() << "sebeta_SNP_A1A1";
+ *outfile[0] << input_var.getSep() << "beta_SNP_A1A2" << input_var.getSep()
+ << "beta_SNP_A1A1" << input_var.getSep() << "sebeta_SNP_A1A2"
+ << input_var.getSep() << "sebeta_SNP_A1A1";
- *outfile[1] << input_var.getSep() << "beta_SNP_addA1" << input_var.getSep()
- << "sebeta_SNP_addA1";
- *outfile[2] << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep()
- << "sebeta_SNP_domA1";
- *outfile[3] << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep()
- << "sebeta_SNP_recA1";
- *outfile[4] << input_var.getSep() << "beta_SNP_odom" << input_var.getSep()
- << "sebeta_SNP_odom";
-//TODO(unknown): compare in create_header_1 and create_header_2 the next lines.
- if (input_var.getInteraction() != 0)
- {
- //Han Chen
- *outfile[0] << input_var.getSep() << "beta_SNP_A1A2_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_A1A2_" << phd.model_terms[interaction_cox]
- << input_var.getSep() << "beta_SNP_A1A1_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_A1A1_" << phd.model_terms[interaction_cox];
+ *outfile[1] << input_var.getSep() << "beta_SNP_addA1" << input_var.getSep()
+ << "sebeta_SNP_addA1";
+ *outfile[2] << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep()
+ << "sebeta_SNP_domA1";
+ *outfile[3] << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep()
+ << "sebeta_SNP_recA1";
+ *outfile[4] << input_var.getSep() << "beta_SNP_odom" << input_var.getSep()
+ << "sebeta_SNP_odom";
+ if (input_var.getInteraction() != 0) {
+ //Han Chen
+ *outfile[0] << input_var.getSep() << "beta_SNP_A1A2_"
+ << phd.model_terms[interaction_cox] << input_var.getSep()
+ << "sebeta_SNP_A1A2_" << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "beta_SNP_A1A1_"
+ << phd.model_terms[interaction_cox] << input_var.getSep()
+ << "sebeta_SNP_A1A1_" << phd.model_terms[interaction_cox];
#if !COXPH
- if (input_var.getInverseFilename() == NULL && !input_var.getAllcov())
- {
- *outfile[0] << input_var.getSep() << "cov_SNP_A1A2_int_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "cov_SNP_A1A1_int_SNP_"
- << phd.model_terms[interaction_cox];
- }
+ if (input_var.getInverseFilename() == NULL && !input_var.getAllcov()) {
+ *outfile[0] << input_var.getSep() << "cov_SNP_A1A2_int_SNP_"
+ << phd.model_terms[interaction_cox] << input_var.getSep()
+ << "cov_SNP_A1A1_int_SNP_"
+ << phd.model_terms[interaction_cox];
+ }
#endif
- //Oct 26, 2009
- for (unsigned int file = 1; file < outfile.size(); file++)
- {
- *outfile[file] << input_var.getSep() << "beta_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_" << phd.model_terms[interaction_cox];
- //Han Chen
+ //Oct 26, 2009
+ for (unsigned int file = 1; file < outfile.size(); file++) {
+ *outfile[file] << input_var.getSep() << "beta_SNP_"
+ << phd.model_terms[interaction_cox] << input_var.getSep()
+ << "sebeta_SNP_" << phd.model_terms[interaction_cox];
+ //Han Chen
#if !COXPH
- if (input_var.getInverseFilename() == NULL
- && !input_var.getAllcov())
- {
- *outfile[file] << input_var.getSep() << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
- }
+ if (input_var.getInverseFilename() == NULL
+ && !input_var.getAllcov()) {
+ *outfile[file] << input_var.getSep() << "cov_SNP_int_SNP_"
+ << phd.model_terms[interaction_cox];
+ }
#endif
- //Oct 26, 2009
- }
- }
- *outfile[0] << input_var.getSep() << "loglik\n"; //"chi2_SNP_2df\n";
- *outfile[1] << input_var.getSep() << "loglik\n"; //"chi2_SNP_A1\n";
- *outfile[2] << input_var.getSep() << "loglik\n"; //"chi2_SNP_domA1\n";
- *outfile[3] << input_var.getSep() << "loglik\n"; //"chi2_SNP_recA1\n";
- *outfile[4] << input_var.getSep() << "loglik\n"; //"chi2_SNP_odom\n";
+ //Oct 26, 2009
+ }
+ }
+ *outfile[0] << input_var.getSep() << "loglik\n"; //"chi2_SNP_2df\n";
+ *outfile[1] << input_var.getSep() << "loglik\n"; //"chi2_SNP_A1\n";
+ *outfile[2] << input_var.getSep() << "loglik\n"; //"chi2_SNP_domA1\n";
+ *outfile[3] << input_var.getSep() << "loglik\n"; //"chi2_SNP_recA1\n";
+ *outfile[4] << input_var.getSep() << "loglik\n"; //"chi2_SNP_odom\n";
}
void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- 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";
+ 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";
-//TODO(unknown): compare in create_header_1 and create_header_2 the next lines.
- 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];
- }
+ 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];
+ }
- if (input_var.getInverseFilename() == NULL)
- //Han Chen
- {
+ if (input_var.getInverseFilename() == NULL)
+ //Han Chen
+ {
#if !COXPH
- if (input_var.getInteraction() != 0 && !input_var.getAllcov())
- {
- *outfile[0] << input_var.getSep() << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
- }
+ if (input_var.getInteraction() != 0 && !input_var.getAllcov()) {
+ *outfile[0] << input_var.getSep() << "cov_SNP_int_SNP_"
+ << phd.model_terms[interaction_cox];
+ }
#endif
- *outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
- }
- //Oct 26, 2009
- *outfile[0] << "\n";
+ *outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
+ }
+ //Oct 26, 2009
+ *outfile[0] << "\n";
}
-int main(int argc, char * argv[])
-{
- cmdvars input_var;
- input_var.set_variables(argc, argv);
+int main(int argc, char * argv[]) {
+ cmdvars input_var;
+ input_var.set_variables(argc, argv);
- input_var.printinfo();
- // if (allcov && ngpreds>1)
- // {
- // fprintf(stdout,"\n\nWARNING: --allcov allowed only for 1 predictor (MLDOSE)\n");
- // allcov = 0;
- // }
- mlinfo mli(input_var.getMlinfofilename(), input_var.getMapfilename());
- int nsnps = mli.nsnps;
- phedata phd;
- int interaction_cox = create_phenoytype(phd, input_var);
+ input_var.printinfo();
+ // if (allcov && ngpreds>1)
+ // {
+ // fprintf(stdout,"\n\nWARNING: --allcov allowed only for 1 predictor (MLDOSE)\n");
+ // allcov = 0;
+ // }
+ mlinfo mli(input_var.getMlinfofilename(), input_var.getMapfilename());
+ int nsnps = mli.nsnps;
+ phedata phd;
+ int interaction_cox = create_phenoytype(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);
- // }
- 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);}
- #endif
- */
+ //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);
+ // }
+ 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);}
+ #endif
+ */
- std::cout << "Reading data ...";
- if (input_var.getInverseFilename() != NULL)
- {
- loadInvSigma(input_var, phd, invvarmatrix);
- // matrix.print();
- }
- std::cout.flush();
- gendata gtd;
- if (!input_var.getIsFvf())
- //use the non non filevector input format
- gtd.re_gendata(input_var.getGenfilename(), nsnps,
- input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
- input_var.getSkipd(), phd.idnames);
- else
- //use the filevector input format (missing second last skipd parameter)
- gtd.re_gendata(input_var.getStrGenfilename(), nsnps,
- input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
- phd.idnames);
+ std::cout << "Reading data ...";
+ if (input_var.getInverseFilename() != NULL) {
+ loadInvSigma(input_var, phd, invvarmatrix);
+ // matrix.print();
+ }
+ std::cout.flush();
+ gendata gtd;
+ if (!input_var.getIsFvf())
+ //use the non non filevector input format
+ gtd.re_gendata(input_var.getGenfilename(), nsnps,
+ input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
+ input_var.getSkipd(), phd.idnames);
+ else
+ //use the filevector input format (missing second last skipd parameter)
+ gtd.re_gendata(input_var.getStrGenfilename(), nsnps,
+ input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
+ phd.idnames);
- std::cout << " loaded genotypic data ...";
- /**
- 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);
- **/
- // estimate null model
- //TODO(maarten): remove this unused variable if there is not a reason to keep it
- //double null_loglik = 0.;
+ std::cout << " loaded genotypic data ...";
+ /**
+ 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);
+ **/
+ // estimate null model
+ //TODO(maarten): remove this unused variable if there is not a reason to keep it
+ //double null_loglik = 0.;
#if COXPH
- coxph_data nrgd = coxph_data(phd, gtd, -1);
+ coxph_data nrgd = coxph_data(phd, gtd, -1);
#else
- regdata nrgd = regdata(phd, gtd, -1, input_var.isIsInteractionExcluded());
+ regdata nrgd = regdata(phd, gtd, -1, input_var.isIsInteractionExcluded());
#endif
- std::cout << " loaded null data ...";
+ std::cout << " loaded null data ...";
#if LOGISTIC
- logistic_reg nrd = logistic_reg(nrgd);
- nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ logistic_reg nrd = logistic_reg(nrgd);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif LINEAR
- linear_reg nrd = linear_reg(nrgd);
+ linear_reg nrd = linear_reg(nrgd);
#if DEBUG
- std::cout << "[DEBUG] linear_reg nrd = linear_reg(nrgd); DONE.";
+ std::cout << "[DEBUG] linear_reg nrd = linear_reg(nrgd); DONE.";
#endif
- nrd.estimate(nrgd, 0, CHOLTOL, 0, input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ nrd.estimate(nrgd, 0, CHOLTOL, 0, input_var.getInteraction(),
+ input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif COXPH
- coxph_reg nrd(nrgd);
+ coxph_reg nrd(nrgd);
- nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), 1);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), 1);
#endif
- std::cout << " estimated null model ...";
- // end null
+ std::cout << " estimated null model ...";
+ // end null
#if COXPH
- coxph_data rgd(phd, gtd, 0);
+ coxph_data rgd(phd, gtd, 0);
#else
- regdata rgd(phd, gtd, 0, input_var.isIsInteractionExcluded());
+ regdata rgd(phd, gtd, 0, input_var.isIsInteractionExcluded());
#endif
- std::cout << " formed regression object ...";
+ std::cout << " formed regression object ...";
- std::cout << " done\n";
- std::cout.flush();
+ std::cout << " done\n";
+ std::cout.flush();
- //________________________________________________________________
- //Maksim, 9 Jan, 2009
+ //________________________________________________________________
+ //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)
- {
- open_files_for_output(outfile, outfilename_str);
- if (input_var.getNohead() != 1)
- {
- create_header_1(outfile, input_var, phd, interaction_cox);
- }
- }
- else //Only additive model. Only one output file
- {
- outfile.push_back(
- new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
+ 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) {
+ open_files_for_output(outfile, outfilename_str);
+ if (input_var.getNohead() != 1) {
+ create_header_1(outfile, input_var, phd, interaction_cox);
+ }
+ } else //Only additive model. Only one output file
+ {
+ outfile.push_back(
+ new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
- if (!outfile[0]->is_open())
- {
- std::cerr << "Can not open file for writing: " << outfilename_str
- << "\n";
- exit(1);
- }
- if (input_var.getNohead() != 1)
- {
- create_header2(outfile, input_var, phd, interaction_cox);
- }
- }
+ if (!outfile[0]->is_open()) {
+ std::cerr << "Can not open file for writing: " << outfilename_str
+ << "\n";
+ exit(1);
+ }
+ if (input_var.getNohead() != 1) {
+ create_header2(outfile, input_var, phd, interaction_cox);
+ }
+ }
- //________________________________________________________________
+ //________________________________________________________________
- /*
- if (input_var.getAllcov())
- {
- if (score)
- {
- outfile << input_var.getSep() << "beta_mu"; // << input_var.getSep() << "beta_SNP_A1";
- outfile << input_var.getSep() << "sebeta_mu"; // << input_var.getSep() << "sebeta_SNP_A1";
- }
- else
- {
- for (int i =0; i<phd.n_model_terms-1;i++)
- outfile << input_var.getSep() << "beta_" << phd.model_terms[i] << input_var.getSep() << "sebeta_" << phd.model_terms[i];
- }
- if(interactio != 0) outfile << input_var.getSep() << "beta_SNP_" << phd.model_terms[interaction];
- }
- 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";
- }
- else
- {
- 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
- int maxmod = 5;
- int start_pos, end_pos;
+ /*
+ if (input_var.getAllcov())
+ {
+ if (score)
+ {
+ outfile << input_var.getSep() << "beta_mu"; // << input_var.getSep() << "beta_SNP_A1";
+ outfile << input_var.getSep() << "sebeta_mu"; // << input_var.getSep() << "sebeta_SNP_A1";
+ }
+ else
+ {
+ for (int i =0; i<phd.n_model_terms-1;i++)
+ outfile << input_var.getSep() << "beta_" << phd.model_terms[i] << input_var.getSep() << "sebeta_" << phd.model_terms[i];
+ }
+ if(interactio != 0) outfile << input_var.getSep() << "beta_SNP_" << phd.model_terms[interaction];
+ }
+ 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";
+ }
+ else
+ {
+ 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
+ int maxmod = 5;
+ int start_pos, end_pos;
- std::vector<std::ostringstream *> beta_sebeta;
- //Han Chen
- std::vector<std::ostringstream *> covvalue;
- //Oct 26, 2009
- std::vector<std::ostringstream *> chi2;
+ std::vector<std::ostringstream *> beta_sebeta;
+ //Han Chen
+ std::vector<std::ostringstream *> covvalue;
+ //Oct 26, 2009
+ std::vector<std::ostringstream *> chi2;
- for (int i = 0; i < maxmod; i++)
- {
- beta_sebeta.push_back(new std::ostringstream());
- //Han Chen
- covvalue.push_back(new std::ostringstream());
- //Oct 26, 2009
- chi2.push_back(new std::ostringstream());
- }
+ for (int i = 0; i < maxmod; i++) {
+ beta_sebeta.push_back(new std::ostringstream());
+ //Han Chen
+ covvalue.push_back(new std::ostringstream());
+ //Oct 26, 2009
+ chi2.push_back(new std::ostringstream());
+ }
- for (int csnp = 0; csnp < nsnps; csnp++)
- {
- rgd.update_snp(gtd, csnp);
- double freq = 0.;
- int gcount = 0;
- float snpdata1[gtd.nids];
- float snpdata2[gtd.nids];
- if (input_var.getNgpreds() == 2)
- {
- //freq = ((gtd.G).column_mean(csnp*2)*2.+(gtd.G).column_mean(csnp*2+1))/2.;
- gtd.get_var(csnp * 2, snpdata1);
- gtd.get_var(csnp * 2 + 1, snpdata2);
- for (unsigned int ii = 0; ii < gtd.nids; ii++)
- if (!isnan(snpdata1[ii]) && !isnan(snpdata2[ii]))
- {
- gcount++;
- freq += snpdata1[ii] + snpdata2[ii] * 0.5;
- }
- }
- else
- {
- // freq = (gtd.G).column_mean(csnp)/2.;
- gtd.get_var(csnp, snpdata1);
- for (unsigned int ii = 0; ii < gtd.nids; ii++)
- if (!isnan(snpdata1[ii]))
- {
- gcount++;
- freq += snpdata1[ii] * 0.5;
- }
- }
- freq /= static_cast<double> (gcount);
- int poly = 1;
- if (fabs(freq) < 1.e-16 || fabs(1. - freq) < 1.e-16)
- poly = 0;
+ for (int csnp = 0; csnp < nsnps; csnp++) {
+ rgd.update_snp(gtd, csnp);
+ double freq = 0.;
+ int gcount = 0;
+ float snpdata1[gtd.nids];
+ float snpdata2[gtd.nids];
+ if (input_var.getNgpreds() == 2) {
+ //freq = ((gtd.G).column_mean(csnp*2)*2.+(gtd.G).column_mean(csnp*2+1))/2.;
+ gtd.get_var(csnp * 2, snpdata1);
+ gtd.get_var(csnp * 2 + 1, snpdata2);
+ for (unsigned int ii = 0; ii < gtd.nids; ii++)
+ if (!isnan(snpdata1[ii]) && !isnan(snpdata2[ii])) {
+ gcount++;
+ freq += snpdata1[ii] + snpdata2[ii] * 0.5;
+ }
+ } else {
+ // freq = (gtd.G).column_mean(csnp)/2.;
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 981
More information about the Genabel-commits
mailing list