[Genabel-commits] r922 - branches/ProbABEL-refactoring/ProbABEL/src pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 14 08:47:58 CEST 2012
Author: lckarssen
Date: 2012-06-14 08:47:57 +0200 (Thu, 14 Jun 2012)
New Revision: 922
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
pkg/ProbABEL/src/main.cpp
Log:
Undo previous commit, because it contained more changes than expected.
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-06-14 06:37:57 UTC (rev 921)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-06-14 06:47:57 UTC (rev 922)
@@ -54,41 +54,41 @@
{
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();
+ 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)
+ 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" };
+ { 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 << "Cannot open file for writing: " << filenames[i]
- << "\n";
- exit(1);
- }
+ 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);
+ }
}
}
@@ -96,25 +96,25 @@
{
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());
+ input_var.getNpeople(), input_var.getInteraction(),
+ input_var.isIscox());
int interaction_cox = input_var.getInteraction();
#if COXPH
interaction_cox--;
#endif
if (input_var.getInteraction() < 0 || input_var.getInteraction() > phd.ncov
- || interaction_cox > phd.ncov)
+ || interaction_cox > phd.ncov)
{
- std::cerr << "error: Interaction parameter is out of range (ineraction="
- << input_var.getInteraction() << ") \n";
- exit(1);
+ std::cerr << "error: Interaction parameter is out of range (ineraction="
+ << input_var.getInteraction() << ") \n";
+ exit(1);
}
return interaction_cox;
}
void loadInvSigma(cmdvars& input_var, phedata& phd,
- mematrix<double>& invvarmatrix)
+ mematrix<double>& invvarmatrix)
{
std::cout << "you are running mmscore...\n";
InvSigma inv(input_var.getInverseFilename(), &phd);
@@ -125,85 +125,85 @@
}
void create_start_of_header(std::vector<std::ofstream*>& outfile,
- cmdvars& input_var, phedata& phd)
+ cmdvars& input_var, phedata& phd)
{
for (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";
+ (*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 (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];
+ for (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)
+ 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";
+ << "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";
+ << "sebeta_SNP_addA1";
*outfile[2] << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep()
- << "sebeta_SNP_domA1";
+ << "sebeta_SNP_domA1";
*outfile[3] << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep()
- << "sebeta_SNP_recA1";
+ << "sebeta_SNP_recA1";
*outfile[4] << input_var.getSep() << "beta_SNP_odom" << input_var.getSep()
- << "sebeta_SNP_odom";
+ << "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];
+ //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 (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 (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
- }
+ //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";
@@ -213,32 +213,32 @@
}
void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- phedata phd, int interaction_cox)
+ 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";
+ << "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];
+ *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 !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";
+ *outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
}
//Oct 26, 2009
*outfile[0] << "\n";
@@ -282,21 +282,21 @@
std::cout << "Reading data ...";
if (input_var.getInverseFilename() != NULL)
{
- loadInvSigma(input_var, phd, invvarmatrix);
- // matrix.print();
+ 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);
+ //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);
+ //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 ...";
/**
@@ -325,7 +325,7 @@
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);
+ input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif COXPH
coxph_reg nrd(nrgd);
@@ -353,27 +353,27 @@
//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);
- }
+ 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()));
+ outfile.push_back(
+ new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
- if (!outfile[0]->is_open())
- {
- std::cerr << "Cannot 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);
+ }
}
//________________________________________________________________
@@ -422,499 +422,499 @@
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());
+ 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 (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 (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;
+ 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 (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 (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;
- if (fabs(mli.Rsq[csnp]) < 1.e-16)
- poly = 0;
- //All models output. One file per each model
- if (input_var.getNgpreds() == 2)
- {
- //Write mlinfo to output:
- for (int file = 0; file < outfile.size(); file++)
- {
- *outfile[file] << mli.name[csnp] << input_var.getSep()
- << mli.A1[csnp] << input_var.getSep() << mli.A2[csnp]
- << input_var.getSep() << mli.Freq1[csnp]
- << input_var.getSep() << mli.MAF[csnp]
- << input_var.getSep() << mli.Quality[csnp]
- << input_var.getSep() << mli.Rsq[csnp]
- << input_var.getSep() << gcount << input_var.getSep()
- << freq;
- if (input_var.getChrom() != "-1")
- *outfile[file] << input_var.getSep()
- << input_var.getChrom();
- if (input_var.getMapfilename() != NULL)
- *outfile[file] << input_var.getSep() << mli.map[csnp];
- }
+ if (fabs(mli.Rsq[csnp]) < 1.e-16)
+ poly = 0;
+ //All models output. One file per each model
+ if (input_var.getNgpreds() == 2)
+ {
+ //Write mlinfo to output:
+ for (int file = 0; file < outfile.size(); file++)
+ {
+ *outfile[file] << mli.name[csnp] << input_var.getSep()
+ << mli.A1[csnp] << input_var.getSep() << mli.A2[csnp]
+ << input_var.getSep() << mli.Freq1[csnp]
+ << input_var.getSep() << mli.MAF[csnp]
+ << input_var.getSep() << mli.Quality[csnp]
+ << input_var.getSep() << mli.Rsq[csnp]
+ << input_var.getSep() << gcount << input_var.getSep()
+ << freq;
+ if (input_var.getChrom() != "-1")
+ *outfile[file] << input_var.getSep()
+ << input_var.getChrom();
+ if (input_var.getMapfilename() != NULL)
+ *outfile[file] << input_var.getSep() << mli.map[csnp];
+ }
- for (int model = 0; model < maxmod; model++)
- {
- if (poly) //allel freq is not to rare
- {
+ for (int model = 0; model < maxmod; model++)
+ {
+ if (poly) //allel freq is not to rare
+ {
#if LOGISTIC
- logistic_reg rd(rgd);
- if (input_var.getScore())
- rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
- else
- rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
+ logistic_reg rd(rgd);
+ if (input_var.getScore())
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ else
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
#elif LINEAR
- linear_reg rd(rgd);
- if (input_var.getScore())
- rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
- input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix);
- else
- {
- // rd.mmscore(rgd,0,CHOLTOL,model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
- rd.estimate(rgd, 0, CHOLTOL, model,
- input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix,
- input_var.getRobust());
- }
+ linear_reg rd(rgd);
+ if (input_var.getScore())
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(), invvarmatrix);
+ else
+ {
+ // rd.mmscore(rgd,0,CHOLTOL,model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ rd.estimate(rgd, 0, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(), invvarmatrix,
+ input_var.getRobust());
+ }
#elif COXPH
- coxph_reg rd(rgd);
- rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), true, input_var.getNgpreds());
+ coxph_reg rd(rgd);
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), true, input_var.getNgpreds());
#endif
- if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() == 0)
- start_pos = rd.beta.nrow - 2;
- else if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() != 0)
- start_pos = rd.beta.nrow - 4;
- else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() == 0)
- start_pos = rd.beta.nrow - 1;
- else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() != 0)
- start_pos = rd.beta.nrow - 2;
- else
- start_pos = 0;
+ if (!input_var.getAllcov() && model == 0
+ && input_var.getInteraction() == 0)
+ start_pos = rd.beta.nrow - 2;
+ else if (!input_var.getAllcov() && model == 0
+ && input_var.getInteraction() != 0)
+ start_pos = rd.beta.nrow - 4;
+ else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() == 0)
+ start_pos = rd.beta.nrow - 1;
+ else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() != 0)
+ start_pos = rd.beta.nrow - 2;
+ else
+ start_pos = 0;
- for (int pos = start_pos; pos < rd.beta.nrow; pos++)
- {
- *beta_sebeta[model] << input_var.getSep()
- << rd.beta[pos] << input_var.getSep()
- << rd.sebeta[pos];
- //Han Chen
+ for (int pos = start_pos; pos < rd.beta.nrow; pos++)
+ {
+ *beta_sebeta[model] << input_var.getSep()
+ << rd.beta[pos] << input_var.getSep()
+ << rd.sebeta[pos];
+ //Han Chen
#if !COXPH
- if (input_var.getInverseFilename() == NULL
- && !input_var.getAllcov()
- && input_var.getInteraction() != 0)
- {
- if (pos > start_pos)
- {
- if (model == 0)
- {
- if (pos > start_pos + 2)
- {
- *covvalue[model]
- << rd.covariance[pos - 3]
- << input_var.getSep()
- << rd.covariance[pos - 2];
- }
- }
- else
- {
- *covvalue[model] << rd.covariance[pos - 1];
- }
- }
- }
+ if (input_var.getInverseFilename() == NULL
+ && !input_var.getAllcov()
+ && input_var.getInteraction() != 0)
+ {
+ if (pos > start_pos)
+ {
+ if (model == 0)
+ {
+ if (pos > start_pos + 2)
+ {
+ *covvalue[model]
+ << rd.covariance[pos - 3]
+ << input_var.getSep()
+ << rd.covariance[pos - 2];
+ }
+ }
+ else
+ {
+ *covvalue[model] << rd.covariance[pos - 1];
+ }
+ }
+ }
#endif
- //Oct 26, 2009
- }
+ //Oct 26, 2009
+ }
- //calculate chi2
- //________________________________
- if (input_var.getScore() == 0)
- {
- //*chi2[model] << 2.*(rd.loglik-null_loglik);
- *chi2[model] << rd.loglik;
- }
- else
- {
- //*chi2[model] << rd.chi2_score;
- *chi2[model] << "nan";
- }
- //________________________________
- }
- else //beta, sebeta = nan
- {
- if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() == 0)
- start_pos = rgd.X.ncol - 2;
- else if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() != 0)
- start_pos = rgd.X.ncol - 4;
- else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() == 0)
- start_pos = rgd.X.ncol - 1;
- else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() != 0)
- start_pos = rgd.X.ncol - 2;
- else
- start_pos = 0;
+ //calculate chi2
+ //________________________________
+ if (input_var.getScore() == 0)
+ {
+ //*chi2[model] << 2.*(rd.loglik-null_loglik);
+ *chi2[model] << rd.loglik;
+ }
+ else
+ {
+ //*chi2[model] << rd.chi2_score;
+ *chi2[model] << "nan";
+ }
+ //________________________________
+ }
+ else //beta, sebeta = nan
+ {
+ if (!input_var.getAllcov() && model == 0
+ && input_var.getInteraction() == 0)
+ start_pos = rgd.X.ncol - 2;
+ else if (!input_var.getAllcov() && model == 0
+ && input_var.getInteraction() != 0)
+ start_pos = rgd.X.ncol - 4;
+ else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() == 0)
+ start_pos = rgd.X.ncol - 1;
+ else if (!input_var.getAllcov() && model != 0
+ && input_var.getInteraction() != 0)
+ start_pos = rgd.X.ncol - 2;
+ else
+ start_pos = 0;
- if (model == 0)
- {
- end_pos = rgd.X.ncol;
- }
- else
- {
- end_pos = rgd.X.ncol - 1;
- }
+ if (model == 0)
+ {
+ end_pos = rgd.X.ncol;
+ }
+ else
+ {
+ end_pos = rgd.X.ncol - 1;
+ }
- if (input_var.getInteraction() != 0)
- end_pos++;
+ if (input_var.getInteraction() != 0)
+ end_pos++;
- for (int pos = start_pos; pos < end_pos; pos++)
- {
- *beta_sebeta[model] << input_var.getSep() << "nan"
- << input_var.getSep() << "nan";
- }
- //Han Chen
+ for (int pos = start_pos; pos < end_pos; pos++)
+ {
+ *beta_sebeta[model] << input_var.getSep() << "nan"
+ << input_var.getSep() << "nan";
+ }
+ //Han Chen
#if !COXPH
- if (!input_var.getAllcov()
- && input_var.getInteraction() != 0)
- {
- if (model == 0)
- {
- *covvalue[model] << "nan" << input_var.getSep()
- << "nan";
- }
- else
- {
- *covvalue[model] << "nan";
- }
- }
+ if (!input_var.getAllcov()
+ && input_var.getInteraction() != 0)
+ {
+ if (model == 0)
+ {
+ *covvalue[model] << "nan" << input_var.getSep()
+ << "nan";
+ }
+ else
+ {
+ *covvalue[model] << "nan";
+ }
+ }
#endif
- //Oct 26, 2009
- *chi2[model] << "nan";
- }
- } //end of model cycle
+ //Oct 26, 2009
+ *chi2[model] << "nan";
+ }
+ } //end of model cycle
- //Han Chen
- *outfile[0] << beta_sebeta[0]->str() << input_var.getSep();
+ //Han Chen
+ *outfile[0] << beta_sebeta[0]->str() << input_var.getSep();
#if !COXPH
- if (!input_var.getAllcov() && input_var.getInteraction() != 0)
- {
- *outfile[0] << covvalue[0]->str() << input_var.getSep();
- }
+ if (!input_var.getAllcov() && input_var.getInteraction() != 0)
+ {
+ *outfile[0] << covvalue[0]->str() << input_var.getSep();
+ }
#endif
- *outfile[0] << chi2[0]->str() << "\n";
+ *outfile[0] << chi2[0]->str() << "\n";
- *outfile[1] << beta_sebeta[1]->str() << input_var.getSep();
+ *outfile[1] << beta_sebeta[1]->str() << input_var.getSep();
#if !COXPH
- if (!input_var.getAllcov() && input_var.getInteraction() != 0)
- {
- *outfile[1] << covvalue[1]->str() << input_var.getSep();
- }
+ if (!input_var.getAllcov() && input_var.getInteraction() != 0)
+ {
+ *outfile[1] << covvalue[1]->str() << input_var.getSep();
+ }
#endif
- *outfile[1] << chi2[1]->str() << "\n";
+ *outfile[1] << chi2[1]->str() << "\n";
- *outfile[2] << beta_sebeta[2]->str() << input_var.getSep();
+ *outfile[2] << beta_sebeta[2]->str() << input_var.getSep();
#if !COXPH
- if (!input_var.getAllcov() && input_var.getInteraction() != 0)
- {
- *outfile[2] << covvalue[2]->str() << input_var.getSep();
- }
+ if (!input_var.getAllcov() && input_var.getInteraction() != 0)
+ {
+ *outfile[2] << covvalue[2]->str() << input_var.getSep();
+ }
#endif
- *outfile[2] << chi2[2]->str() << "\n";
+ *outfile[2] << chi2[2]->str() << "\n";
- *outfile[3] << beta_sebeta[3]->str() << input_var.getSep();
+ *outfile[3] << beta_sebeta[3]->str() << input_var.getSep();
#if !COXPH
- if (!input_var.getAllcov() && input_var.getInteraction() != 0)
- {
- *outfile[3] << covvalue[3]->str() << input_var.getSep();
- }
+ if (!input_var.getAllcov() && input_var.getInteraction() != 0)
+ {
+ *outfile[3] << covvalue[3]->str() << input_var.getSep();
+ }
#endif
- *outfile[3] << chi2[3]->str() << "\n";
+ *outfile[3] << chi2[3]->str() << "\n";
- *outfile[4] << beta_sebeta[4]->str() << input_var.getSep();
+ *outfile[4] << beta_sebeta[4]->str() << input_var.getSep();
#if !COXPH
- if (!input_var.getAllcov() && input_var.getInteraction() != 0)
- {
- *outfile[4] << covvalue[4]->str() << input_var.getSep();
- }
+ if (!input_var.getAllcov() && input_var.getInteraction() != 0)
+ {
+ *outfile[4] << covvalue[4]->str() << input_var.getSep();
+ }
#endif
- *outfile[4] << chi2[4]->str() << "\n";
- //Oct 26, 2009
- }
- else //Only additive model. Only one output file
- {
- //Write mlinfo to output:
- *outfile[0] << mli.name[csnp] << input_var.getSep() << mli.A1[csnp]
- << input_var.getSep() << mli.A2[csnp] << input_var.getSep();
- *outfile[0] << mli.Freq1[csnp] << input_var.getSep()
- << mli.MAF[csnp] << input_var.getSep() << mli.Quality[csnp]
- << input_var.getSep() << mli.Rsq[csnp]
- << input_var.getSep();
- *outfile[0] << gcount << input_var.getSep() << freq;
- if (input_var.getChrom() != "-1")
- *outfile[0] << input_var.getSep() << input_var.getChrom();
- if (input_var.getMapfilename() != NULL)
- *outfile[0] << input_var.getSep() << mli.map[csnp];
- int model = 0;
- if (poly) //allel freq is not to rare
- {
+ *outfile[4] << chi2[4]->str() << "\n";
+ //Oct 26, 2009
+ }
+ else //Only additive model. Only one output file
+ {
+ //Write mlinfo to output:
+ *outfile[0] << mli.name[csnp] << input_var.getSep() << mli.A1[csnp]
+ << input_var.getSep() << mli.A2[csnp] << input_var.getSep();
+ *outfile[0] << mli.Freq1[csnp] << input_var.getSep()
+ << mli.MAF[csnp] << input_var.getSep() << mli.Quality[csnp]
+ << input_var.getSep() << mli.Rsq[csnp]
+ << input_var.getSep();
+ *outfile[0] << gcount << input_var.getSep() << freq;
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 922
More information about the Genabel-commits
mailing list