[Genabel-commits] r970 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 4 00:40:55 CEST 2012
Author: lckarssen
Date: 2012-10-04 00:40:55 +0200 (Thu, 04 Oct 2012)
New Revision: 970
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.h
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
Log:
In the ProbABEL-refactoring branch: forward-ported the changes from SVN r.968: Fixed a lot of compiler warnings. Mostly by replacing ints with unsigned ints, for example for nids, nsnps. This makes perfect sense and serves as extra sanity check on our data.
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-10-03 22:24:14 UTC (rev 969)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-10-03 22:40:55 UTC (rev 970)
@@ -33,6 +33,9 @@
{
return 0;
}
+
+ // You should never come here...
+ return -9;
}
coxph_data::coxph_data(const coxph_data &obj)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-10-03 22:24:14 UTC (rev 969)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-10-03 22:40:55 UTC (rev 970)
@@ -19,20 +19,20 @@
void gendata::get_var(int var, float * data)
{
if (DAG == NULL)
- for (int i = 0; i < G.nrow; i++)
- data[i] = G.get(i, var);
+ for (int i = 0; i < G.nrow; i++)
+ data[i] = G.get(i, var);
else if (DAG != NULL)
{
- float tmpdata[DAG->getNumObservations()];
- DAG->readVariableAs((unsigned long int) var, tmpdata);
- unsigned int j = 0;
- for (unsigned int i = 0; i < DAG->getNumObservations(); i++)
- if (!DAGmask[i])
- data[j++] = tmpdata[i];
- //fprintf(stdout,"%i %i %i\n",j,DAG->get_nobservations(),nids);
+ float tmpdata[DAG->getNumObservations()];
+ DAG->readVariableAs((unsigned long int) var, tmpdata);
+ unsigned int j = 0;
+ for (unsigned int i = 0; i < DAG->getNumObservations(); i++)
+ if (!DAGmask[i])
+ data[j++] = tmpdata[i];
+ //fprintf(stdout,"%i %i %i\n",j,DAG->get_nobservations(),nids);
}
else
- report_error("cannot get gendata");
+ report_error("cannot get gendata");
}
gendata::gendata()
@@ -40,51 +40,57 @@
nsnps = nids = ngpreds = 0;
}
-void gendata::re_gendata(string filename, int insnps, int ingpreds, int npeople,
- int nmeasured, unsigned short int * allmeasured, std::string * idnames)
+void gendata::re_gendata(string filename, unsigned int insnps,
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured,
+ std::string * idnames)
{
nsnps = insnps;
ngpreds = ingpreds;
DAG = new FileVector(filename, 128, true);
DAGmask = new unsigned short int[DAG->getNumObservations()];
if (DAG->getNumObservations() != npeople)
- report_error("dimension of fvf-data and phenotype data do not match\n");
+ report_error("dimension of fvf-data and phenotype data do not match\n");
if (DAG->getNumVariables() != insnps * ingpreds)
- report_error("dimension of fvf-data and mlinfo data do not match\n");
+ report_error("dimension of fvf-data and mlinfo data do not match\n");
long int j = -1;
for (unsigned int i = 0; i < npeople; i++)
{
- if (allmeasured[i] == 0)
- DAGmask[i] = 1;
- else
- {
- DAGmask[i] = 0;
- j++;
- }
- string DAGobsname = DAG->readObservationName(i).name;
+ if (allmeasured[i] == 0)
+ DAGmask[i] = 1;
+ else
+ {
+ DAGmask[i] = 0;
+ j++;
+ }
+ string DAGobsname = DAG->readObservationName(i).name;
- if (DAGobsname.find("->") != string::npos)
- DAGobsname = DAGobsname.substr(DAGobsname.find("->") + 2);
+ if (DAGobsname.find("->") != string::npos)
+ DAGobsname = DAGobsname.substr(DAGobsname.find("->") + 2);
//if (allmeasured[i] && idnames[j] != DAGobsname)
- // error("names do not match for observation at phenofile line (phe/geno) %i/+1 (%s/%s)\n",
- // i+1,idnames[i].c_str(),DAGobsname.c_str());
- // fix thanks to Vadym Pinchuk
- if (allmeasured[i] && idnames[j] != DAGobsname)
- report_error(
- "names do not match for observation at phenofile line(phe/geno) %i/+1 (%s/%s)\n",
- i + 1, idnames[j].c_str(), DAGobsname.c_str());
+ // error("names do not match for observation at phenofile line (phe/geno) %i/+1 (%s/%s)\n",
+ // i+1,idnames[i].c_str(),DAGobsname.c_str());
+ // fix thanks to Vadym Pinchuk
+ if (allmeasured[i] && idnames[j] != DAGobsname)
+ report_error(
+ "names do not match for observation at phenofile line(phe/geno) %i/+1 (%s/%s)\n",
+ i + 1, idnames[j].c_str(), DAGobsname.c_str());
}
nids = j + 1;
//fprintf(stdout,"in INI: %i %i\n",nids,npeople);
if (nids != nmeasured)
- report_error("nids != mneasured (%i != %i)\n", nids, nmeasured);
+ report_error("nids != mneasured (%i != %i)\n", nids, nmeasured);
}
-void gendata::re_gendata(char * fname, int insnps, int ingpreds, int npeople,
- int nmeasured, unsigned short int * allmeasured, int skipd,
- std::string * idnames)
+void gendata::re_gendata(char * fname, unsigned int insnps,
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured,
+ int skipd,
+ std::string * idnames)
{
nids = nmeasured;
nsnps = insnps;
@@ -99,74 +105,74 @@
infile.open(fname);
if (!infile)
{
- std::cerr << "gendata: cannot open file " << fname << endl;
+ std::cerr << "gendata: cannot open file " << fname << endl;
}
char tmp[100], tmpn[100];
std::string tmpid, tmpstr;
int k = 0;
- for (int i = 0; i < npeople; i++)
- if (allmeasured[i] == 1)
- {
- if (skipd > 0)
- {
- // int ttt;
- char ttt[100];
- infile >> tmp;
- // sscanf(tmp,"%d->%s",&ttt, tmpn);
- // these changes are thanks to BMM & BP :)
- // sscanf(tmp,"%s->%s",&ttt, tmpn);
- // sscanf(tmp,"%[^->]->%[^->]",&ttt, tmpn);
- tmpstr = tmp;
- if (tmpstr.find("->") != string::npos)
- {
- sscanf(tmp, "%[^->]->%s", ttt, tmpn);
- tmpid = tmpn;
- }
- else
- {
- tmpid = tmpstr;
- //fprintf(stdout,"%s;%s;%s;%s;%s\n",tmp,ttt,tmpn,tmpid.c_str(),idnames[k].c_str());
- }
- if (tmpid != idnames[k])
- {
- fprintf(stderr,
- "phenofile and dosefile did not match at line %d ",
- i + 2);
- cerr << "(" << tmpid << " != " << idnames[k] << ")\n";
- infile.close();
- exit(1);
- }
- }
- for (int j = 1; j < skipd; j++)
- {
- infile >> tmp;
- }
- for (int j = 0; j < (nsnps * ngpreds); j++)
- {
- if (infile.good())
- {
- infile >> tmp;
- }
- else
- {
- std::cerr
- << "cannot read dose-file: check skipd and ngpreds parameters\n";
- infile.close();
- exit(1);
- }
- G.put(atof(tmp), k, j);
- }
- k++;
- }
- else
- {
- for (int j = 0; j < skipd; j++)
- infile >> tmp;
- for (int j = 0; j < (nsnps * ngpreds); j++)
- infile >> tmp;
- }
+ for (unsigned int i = 0; i < npeople; i++)
+ if (allmeasured[i] == 1)
+ {
+ if (skipd > 0)
+ {
+ // int ttt;
+ char ttt[100];
+ infile >> tmp;
+ // sscanf(tmp,"%d->%s",&ttt, tmpn);
+ // these changes are thanks to BMM & BP :)
+ // sscanf(tmp,"%s->%s",&ttt, tmpn);
+ // sscanf(tmp,"%[^->]->%[^->]",&ttt, tmpn);
+ tmpstr = tmp;
+ if (tmpstr.find("->") != string::npos)
+ {
+ sscanf(tmp, "%[^->]->%s", ttt, tmpn);
+ tmpid = tmpn;
+ }
+ else
+ {
+ tmpid = tmpstr;
+ //fprintf(stdout,"%s;%s;%s;%s;%s\n",tmp,ttt,tmpn,tmpid.c_str(),idnames[k].c_str());
+ }
+ if (tmpid != idnames[k])
+ {
+ fprintf(stderr,
+ "phenofile and dosefile did not match at line %d ",
+ i + 2);
+ cerr << "(" << tmpid << " != " << idnames[k] << ")\n";
+ infile.close();
+ exit(1);
+ }
+ }
+ for (int j = 1; j < skipd; j++)
+ {
+ infile >> tmp;
+ }
+ for (int j = 0; j < (nsnps * ngpreds); j++)
+ {
+ if (infile.good())
+ {
+ infile >> tmp;
+ }
+ else
+ {
+ std::cerr
+ << "cannot read dose-file: check skipd and ngpreds parameters\n";
+ infile.close();
+ exit(1);
+ }
+ G.put(atof(tmp), k, j);
+ }
+ k++;
+ }
+ else
+ {
+ for (int j = 0; j < skipd; j++)
+ infile >> tmp;
+ for (int j = 0; j < (nsnps * ngpreds); j++)
+ infile >> tmp;
+ }
infile.close();
}
// HERE NEED A NEW CONSTRUCTOR BASED ON DATABELBASECPP OBJECT
@@ -175,10 +181,9 @@
if (DAG != NULL)
{
- delete DAG;
- delete[] DAGmask;
+ delete DAG;
+ delete[] DAGmask;
}
// delete G;
}
-
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-10-03 22:24:14 UTC (rev 969)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-10-03 22:40:55 UTC (rev 970)
@@ -19,16 +19,21 @@
class gendata
{
public:
- int nsnps;
- int nids;
- int ngpreds;
+ unsigned int nsnps;
+ unsigned int nids;
+ unsigned int ngpreds;
gendata();
- void re_gendata(char * fname, int insnps, int ingpreds, int npeople,
- int nmeasured, unsigned short int * allmeasured, int skipd,
- std::string * idnames);
- void re_gendata(string filename, int insnps, int ingpreds, int npeople,
- int nmeasured, unsigned short int * allmeasured,
- std::string * idnames);
+ void re_gendata(char * fname, unsigned int insnps,
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured,
+ int skipd,
+ std::string * idnames);
+ void re_gendata(string filename, unsigned int insnps,
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured,
+ std::string * idnames);
void get_var(int var, float * data);
~gendata();
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-03 22:24:14 UTC (rev 969)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-03 22:40:55 UTC (rev 970)
@@ -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 << "Can not 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++)
+ 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";
+ (*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 (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)
+ 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 (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
- }
+ //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 << "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);
+ }
}
//________________________________________________________________
@@ -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 (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;
- 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 (unsigned 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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 970
More information about the Genabel-commits
mailing list