[Genabel-commits] r973 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 6 00:28:34 CEST 2012
Author: maartenk
Date: 2012-10-06 00:28:33 +0200 (Sat, 06 Oct 2012)
New Revision: 973
Added:
branches/ProbABEL-refactoring/ProbABEL/src/maskedmatrix.cpp
branches/ProbABEL-refactoring/ProbABEL/src/maskedmatrix.h
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/mematrix.h
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
Log:
-separate invarmatrix from rest of code (target for optimalization) wile making this class I found the bug described in svn commit 965 (and fix this)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp 2012-10-05 21:46:18 UTC (rev 972)
+++ branches/ProbABEL-refactoring/ProbABEL/src/comand_line_settings.cpp 2012-10-05 22:28:33 UTC (rev 973)
@@ -244,6 +244,9 @@
fprintf(stdout,
"%s v. %s (C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, EMCR\n\n",
PACKAGE, PACKAGE_VERSION);
+#if EIGEN
+ fprintf(stdout, "using EIGEN for matrix operations\n");
+#endif
if (neco[0] != 1 || neco[1] != 1 || neco[2] != 1)
{
Modified: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h 2012-10-05 21:46:18 UTC (rev 972)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h 2012-10-05 22:28:33 UTC (rev 973)
@@ -86,7 +86,7 @@
// return temp;
//}
template<class DT>
-mematrix<DT> mematrix<DT>::operator+(mematrix<DT> &M)
+mematrix<DT> mematrix<DT>::operator+(const mematrix<DT> &M)
{
if (ncol != M.ncol || nrow != M.nrow)
{
@@ -102,14 +102,14 @@
}
template<class DT>
-mematrix<DT> mematrix<DT>::operator-(DT toadd)
+mematrix<DT> mematrix<DT>::operator-(const DT toadd)
{
mematrix<DT> temp(nrow, ncol);
temp.data = data.array() - toadd;
return temp;
}
template<class DT>
-mematrix<DT> mematrix<DT>::operator-(mematrix<DT> &M)
+mematrix<DT> mematrix<DT>::operator-(const mematrix<DT> &M)
{
if (ncol != M.ncol || nrow != M.nrow)
{
@@ -127,7 +127,7 @@
}
template<class DT>
-mematrix<DT> mematrix<DT>::operator*(DT multiplyer)
+mematrix<DT> mematrix<DT>::operator*( DT multiplyer)
{
// MatrixXd add = MatrixXd::Constant(nrow, ncol, toadd);
mematrix<DT> temp(nrow, ncol);
@@ -136,7 +136,7 @@
}
template<class DT>
-mematrix<DT> mematrix<DT>::operator*(mematrix<DT> &M)
+mematrix<DT> mematrix<DT>::operator*(const mematrix<DT> &M)
{
if (ncol != M.nrow)
{
@@ -153,6 +153,27 @@
return temp;
}
+
+template<class DT>
+mematrix<DT> mematrix<DT>::operator*(const mematrix<DT> *M)
+{
+ if (ncol != M->nrow)
+ {
+ fprintf(stderr, "mematrix*: ncol != nrow (%d,%d) and (%d,%d)", nrow,
+ ncol, M->nrow, M->ncol);
+
+ }
+ mematrix<DT> temp;
+ temp.data = data * M->data;
+ temp.ncol=temp.data.cols();
+ temp.nrow=temp.data.rows();
+ temp.nelements=temp.nrow* temp.ncol;
+
+ return temp;
+}
+
+
+
//
// operations
//
@@ -272,20 +293,21 @@
return temp;
}
-//template<class DT>
-//mematrix<DT> reorder(mematrix<DT> &M, mematrix<int> order)
-//{
-// if (M.nrow != order.nrow)
-// {
-// fprintf(stderr, "reorder: M & order have differet # of rows\n");
-// exit(1);
-// }
-// mematrix<DT> temp(M.nrow, M.ncol);
+template<class DT>
+mematrix<DT> reorder(mematrix<DT> &M, mematrix<int> order)
+{
+ if (M.nrow != order.nrow)
+ {
+ fprintf(stderr, "reorder: M & order have differet # of rows\n");
+ exit(1);
+ }
+ mematrix<DT> temp(M.nrow, M.ncol);
+//FIXME: commented out to get compilation running
// for (int i = 0; i < temp.nrow; i++)
// for (int j = 0; j < temp.ncol; j++)
// temp.data[order[i] * temp.ncol + j] = M.data[i * M.ncol + j];
-// return temp;
-//}
+ return temp;
+}
//
//
//template<class DT>
Modified: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h 2012-10-05 21:46:18 UTC (rev 972)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h 2012-10-05 22:28:33 UTC (rev 973)
@@ -33,12 +33,14 @@
mematrix & operator=(const mematrix &M);
DT & operator[]( int i);
// mematrix operator+(DT toadd);
- mematrix operator+(mematrix &M);
+ mematrix operator+(const mematrix &M);
mematrix operator-(DT toadd);
- mematrix operator-(mematrix &M);
+ mematrix operator-(const mematrix &M);
mematrix operator*(DT toadd);
- mematrix operator*(mematrix &M);
+ mematrix operator*(const mematrix &M);
+ mematrix operator*(const mematrix *M);
+
void reinit(int nr, int nc);
unsigned int getnrow(void)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-05 21:46:18 UTC (rev 972)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-10-05 22:28:33 UTC (rev 973)
@@ -42,6 +42,7 @@
#include "mematrix.h"
#include "mematri1.h"
#endif
+#include "maskedmatrix.h"
#include "data.h"
#include "reg1.h"
#include "comand_line_settings.h"
@@ -54,41 +55,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,49 +97,49 @@
{
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)
+ 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 = invvarmatrix * par;
+ // 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)
+ 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";
+ (*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
@@ -152,58 +153,58 @@
}
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
+ //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
+ {
+ *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 +214,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";
@@ -271,7 +272,7 @@
// std::cerr<<"Error: In mmscore you can use additive model without any inetractions only\n";
// exit(1);
// }
- mematrix<double> invvarmatrix;
+ masked_matrix invvarmatrix;
/*
* now should be possible... delete this part later when everything works
#if LOGISTIC
@@ -282,21 +283,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 +326,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 +354,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 +423,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);
+ 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);
+ 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 (!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:
+ 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];
- }
+ {
+ *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
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 973
More information about the Genabel-commits
mailing list