[Genabel-commits] r1032 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 28 11:58:40 CET 2012
Author: lckarssen
Date: 2012-11-28 11:58:40 +0100 (Wed, 28 Nov 2012)
New Revision: 1032
Modified:
pkg/ProbABEL/src/command_line_settings.cpp
pkg/ProbABEL/src/command_line_settings.h
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/coxph_data.h
pkg/ProbABEL/src/data.cpp
pkg/ProbABEL/src/data.h
pkg/ProbABEL/src/eigen_mematrix.h
pkg/ProbABEL/src/extract-snp.cpp
pkg/ProbABEL/src/gendata.cpp
pkg/ProbABEL/src/gendata.h
pkg/ProbABEL/src/maskedmatrix.cpp
pkg/ProbABEL/src/mematri1.h
pkg/ProbABEL/src/reg1.cpp
pkg/ProbABEL/src/reg1.h
pkg/ProbABEL/src/regdata.h
pkg/ProbABEL/src/testchol.cpp
pkg/ProbABEL/src/usage.cpp
Log:
ProbABEL: improving code layout. Should not result in any functional changes.
Modified: pkg/ProbABEL/src/command_line_settings.cpp
===================================================================
--- pkg/ProbABEL/src/command_line_settings.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/command_line_settings.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -241,8 +241,9 @@
void cmdvars::printinfo()
{
cout << PACKAGE
- << " v. " << PACKAGE_VERSION
- << "(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, EMCR\n\n";
+ << " v. " << PACKAGE_VERSION
+ << "\n(C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, "
+ << "EMCR\n\n";
#if EIGEN
cout << "Using EIGEN for matrix operations\n";
#endif
@@ -334,7 +335,8 @@
#if COXPH
if (score)
{
- cerr << "\n\nOption --score is implemented for linear and logistic models only\n"
+ cerr << "\n\nOption --score is implemented for "
+ << "linear and logistic models only\n"
<< endl;
exit(1);
}
@@ -347,7 +349,8 @@
}
if (robust)
{
- cerr << "ERROR: robust standard errors not implemented for Cox regression"
+ cerr << "ERROR: robust standard errors not implemented "
+ << "for Cox regression"
<< endl;
exit(1);
}
Modified: pkg/ProbABEL/src/command_line_settings.h
===================================================================
--- pkg/ProbABEL/src/command_line_settings.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/command_line_settings.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -7,13 +7,13 @@
#ifndef COMMAND_LINE_SETTINGS_H_
#define COMMAND_LINE_SETTINGS_H_
+#include <string>
using namespace std;
class cmdvars
{
private:
-
char * program_name;
char *phefilename;
@@ -48,16 +48,16 @@
program_name = NULL;
std::fill_n(neco, 3, 0);
- phefilename = NULL;
- mlinfofilename = NULL;
- genfilename = NULL;
- mapfilename = NULL;
- outfilename = NULL;
+ phefilename = NULL;
+ mlinfofilename = NULL;
+ genfilename = NULL;
+ mapfilename = NULL;
+ outfilename = NULL;
inverse_filename = NULL;
- sep = " ";
- nohead = 0;
- score = 0;
+ sep = " ";
+ nohead = 0;
+ score = 0;
npeople = -1;
ngpreds = 1;
interaction = 0;
@@ -66,20 +66,18 @@
robust = 0;
chrom = "-1";
str_genfilename = "";
-
- isFVF = 0;
-
- skipd = 2;
+ isFVF = 0;
+ skipd = 2;
allcov = 0;
#if COXPH
noutcomes = 2;
- iscox=true;
+ iscox = true;
#else
noutcomes = 1;
- iscox = false;
+ iscox = false;
#endif
-
}
+
void set_variables(int, char *[]);
char* getPhefilename() const;
int getAllcov() const;
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/coxph_data.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -95,8 +95,9 @@
sstat[i] = int((phed.Y).get(i, 1));
if (sstat[i] != 1 && sstat[i] != 0)
{
- std::cerr << "coxph_data: status not 0/1 (correct order: id, fuptime, status ...)"
- << endl;
+ std::cerr << "coxph_data: status not 0/1 "
+ <<"(correct order: id, fuptime, status ...)"
+ << endl;
exit(1);
}
}
@@ -114,9 +115,9 @@
X.put(snpdata[i], i, (ncov - ngpreds + j));
}
- // for (int i=0;i<nids;i++)
- // for (int j=0;j<ngpreds;j++)
- // X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+j));
+ // for (int i=0;i<nids;i++)
+ // for (int j=0;j<ngpreds;j++)
+ // X.put((gend.G).get(i,(snpnum*ngpreds+j)),i,(ncov-ngpreds+j));
for (int i = 0; i < nids; i++)
{
@@ -127,6 +128,7 @@
// sort by time
double tmptime[nids];
int passed_sorted[nids];
+
for (int i = 0; i < nids; i++)
{
tmptime[i] = stime[i];
@@ -305,7 +307,8 @@
int flag;
double sctest = 1.0;
- //TODO(maarten): this function works only in combination with eigen remove .data() from arguments to make is available under old matrix
+ //TODO(maarten): this function works only in combination with eigen
+ // remove .data() from arguments to make is available under old matrix
coxfit2(&maxiter, &cdata.nids, &X.nrow, cdata.stime.data.data(),
cdata.sstat.data.data(), X.data.data(), newoffset.data.data(),
cdata.weights.data.data(), cdata.strata.data.data(),
Modified: pkg/ProbABEL/src/coxph_data.h
===================================================================
--- pkg/ProbABEL/src/coxph_data.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/coxph_data.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -22,7 +22,6 @@
class coxph_data {
public:
-
coxph_data get_unmasked_data();
coxph_data()
{
@@ -48,6 +47,7 @@
mematrix<int> order;
unsigned short int * masked_data;
};
+
class coxph_reg {
public:
mematrix<double> beta;
Modified: pkg/ProbABEL/src/data.cpp
===================================================================
--- pkg/ProbABEL/src/data.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/data.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -69,7 +69,6 @@
mlinfo::mlinfo(char * filename, char * mapname)
{
-
char tmp[100];
unsigned int nlin = 0;
std::ifstream infile(filename);
@@ -95,14 +94,14 @@
}
nsnps = int(nlin / 7) - 1;
std::cout << "Number of SNPs = " << nsnps << endl;
- name = new std::string[nsnps];
- A1 = new std::string[nsnps];
- A2 = new std::string[nsnps];
- Freq1 = new double[nsnps];
- MAF = new double[nsnps];
+ name = new std::string[nsnps];
+ A1 = new std::string[nsnps];
+ A2 = new std::string[nsnps];
+ Freq1 = new double[nsnps];
+ MAF = new double[nsnps];
Quality = new double[nsnps];
- Rsq = new double[nsnps];
- map = new std::string[nsnps];
+ Rsq = new double[nsnps];
+ map = new std::string[nsnps];
infile.open(filename);
if (!infile)
@@ -185,7 +184,6 @@
unsigned row = 0;
while (myfile.getline(line, MAXIMUM_PEOPLE_AMOUNT))
{
-
std::stringstream line_stream(line);
line_stream >> id;
@@ -194,7 +192,8 @@
std::cerr << "error:in row " << row << " id="
<< phe->idnames[row]
<< " in inverse variance matrix but id=" << id
- << " must be there. Wrong inverse variance matrix (only measured id must be there)\n";
+ << " must be there. Wrong inverse variance matrix"
+ << " (only measured id must be there)\n";
exit(1);
}
unsigned col = 0;
@@ -224,11 +223,10 @@
delete[] line;
}
-;
+
InvSigma::~InvSigma()
{
-//af af
}
mematrix<double> & InvSigma::get_matrix(void)
Modified: pkg/ProbABEL/src/data.h
===================================================================
--- pkg/ProbABEL/src/data.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/data.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -7,6 +7,7 @@
#ifndef DATA_H_
#define DATA_H_
+#include <string>
extern bool is_interaction_excluded;
@@ -42,12 +43,12 @@
};
class InvSigma {
-
private:
static const unsigned MAXIMUM_PEOPLE_AMOUNT = 1000000;
unsigned npeople; //amount of people
std::string filename;
mematrix<double> matrix; //file is stored here
+
public:
InvSigma(const char * filename_, phedata * phe);
mematrix<double> & get_matrix();
Modified: pkg/ProbABEL/src/eigen_mematrix.h
===================================================================
--- pkg/ProbABEL/src/eigen_mematrix.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/eigen_mematrix.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -1,6 +1,7 @@
#ifndef __EIGEN_MEMATRIX_H__
#define __EIGEN_MEMATRIX_H__
#include <Eigen/Dense>
+#include <Eigen/LU>
#include <iostream>
using namespace Eigen;
Modified: pkg/ProbABEL/src/extract-snp.cpp
===================================================================
--- pkg/ProbABEL/src/extract-snp.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/extract-snp.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -21,11 +21,11 @@
******************************************************************************/
#include <stdio.h>
+#include <getopt.h>
#include <string>
#include <iostream>
#include <fstream>
#include <map>
-#include <getopt.h>
using namespace std;
@@ -37,16 +37,16 @@
void info(char *program_name)
{
cout << program_name
- << " extracts the dosage for a SNP for one or more individuals from a"
- << "file in filevector format (.fvi/.fvd)." << endl;
+ << " extracts the dosage for a SNP for one or more individuals from a"
+ << "file in filevector format (.fvi/.fvd)." << endl;
cout << endl;
cout << "Usage: " << program_name << " --file <fv file> --snp <snpname>"
- << endl;
+ << endl;
cout << " or: " << program_name << " -f <fv file> -s <snpname>" << endl;
cout << endl;
cout << "Additional options:" << endl;
- cout << "\t--id (-i) <idname>: only print the dosage for the given SNP and ID"
- << endl;
+ cout << "\t--id (-i) <idname>: "
+ << "only print the dosage for the given SNP and ID" << endl;
cout << "\t--debug (-d): show debugging output" << endl;
cout << "\t--help (-h): show this information" << endl;
}
@@ -56,21 +56,21 @@
{
if (argc < 2)
{
- info(argv[0]);
- exit(3);
+ info(argv[0]);
+ exit(3);
}
int next_option;
const char * const short_options = "df:hi:s:";
const struct option long_options [] =
- {
- {"debug", 0, NULL, 'd'},
- {"file", 1, NULL, 'f'},
- {"help", 0, NULL, 'h'},
- {"id", 1, NULL, 'i'},
- {"snp", 1, NULL, 's'},
- {NULL , 0, NULL, 0 }
- };
+ {
+ {"debug", 0, NULL, 'd'},
+ {"file", 1, NULL, 'f'},
+ {"help", 0, NULL, 'h'},
+ {"id", 1, NULL, 'i'},
+ {"snp", 1, NULL, 's'},
+ {NULL , 0, NULL, 0 }
+ };
char *program_name = argv[0];
bool debug = false;
string inputFileName = "";
@@ -78,35 +78,35 @@
string snpname = "";
do
{
- next_option = getopt_long(argc, argv, short_options, long_options, NULL);
- switch (next_option)
- {
- case 'h':
- info(program_name);
- exit(0);
- case 'd':
- debug = true;
- break;
- case 'f':
- inputFileName = string(optarg);
- break;
- case 'i':
- idname = string(optarg);
- break;
- case 's':
- snpname = string(optarg);
- break;
- case '?': break;
- case -1 : break;
- }
+ next_option = getopt_long(argc, argv, short_options, long_options, NULL);
+ switch (next_option)
+ {
+ case 'h':
+ info(program_name);
+ exit(0);
+ case 'd':
+ debug = true;
+ break;
+ case 'f':
+ inputFileName = string(optarg);
+ break;
+ case 'i':
+ idname = string(optarg);
+ break;
+ case 's':
+ snpname = string(optarg);
+ break;
+ case '?': break;
+ case -1 : break;
+ }
}
while (next_option != -1);
if (snpname.compare("") == 0)
{
- cerr << "Error: No SNP name given" << endl << endl;
- info(program_name);
- exit(2);
+ cerr << "Error: No SNP name given" << endl << endl;
+ info(program_name);
+ exit(2);
}
unsigned long int snprow = 0;
@@ -119,90 +119,90 @@
if (debug)
{
- for (unsigned long int col=0; col<fv.getNumObservations(); col++)
- {
- cout << fv.readObservationName(col).name << " ";
- }
+ for (unsigned long int col = 0; col < fv.getNumObservations(); col++)
+ {
+ cout << fv.readObservationName(col).name << " ";
+ }
- cout << endl;
- cout << "----------------------" << endl;
+ cout << endl;
+ cout << "----------------------" << endl;
}
// Look at the SNPs (rows) first
- for (unsigned long int row=0; row<fv.getNumVariables(); row++)
+ for (unsigned long int row = 0; row < fv.getNumVariables(); row++)
{
- string current_snpname = fv.readVariableName(row).name;
- if (current_snpname.compare(snpname) == 0)
- {
- snpfound = true;
- snprow = row;
- if (debug)
- {
- cout << "*" << current_snpname << "* ";
- }
- } else {
- if (debug)
- {
- cout << current_snpname << " ";
- }
- }
+ string current_snpname = fv.readVariableName(row).name;
+ if (current_snpname.compare(snpname) == 0)
+ {
+ snpfound = true;
+ snprow = row;
+ if (debug)
+ {
+ cout << "*" << current_snpname << "* ";
+ }
+ } else {
+ if (debug)
+ {
+ cout << current_snpname << " ";
+ }
+ }
}
if (debug)
{
- cout << endl;
- cout << "----------------------" << endl;
+ cout << endl;
+ cout << "----------------------" << endl;
- cout << "N_obs = " << fv.getNumObservations() << "\telement size: "
- << fv.getElementSize() << "\tDataType: "
- << dataTypeToString(fv.getElementType()) << endl;
+ cout << "N_obs = " << fv.getNumObservations() << "\telement size: "
+ << fv.getElementSize() << "\tDataType: "
+ << dataTypeToString(fv.getElementType()) << endl;
}
- if(!snpfound)
+ if (!snpfound)
{
- cerr << "SNP name '" << snpname << "' not found in data file "
- << inputFileName << endl;
- exit(1);
+ cerr << "SNP name '" << snpname << "' not found in data file "
+ << inputFileName << endl;
+ exit(1);
}
char * data = new (nothrow) char[fv.getNumObservations() *
- fv.getElementSize()];
+ fv.getElementSize()];
if (!data)
{
- cerr << "Cannot allocate memory for data vector. Exiting..." << endl;
- exit(2);
+ cerr << "Cannot allocate memory for data vector. Exiting..." << endl;
+ exit(2);
}
fv.readVariable(snprow, data);
- for (unsigned long int col=0; col<fv.getNumObservations(); col++)
+ for (unsigned long int col = 0; col < fv.getNumObservations(); col++)
{
- if (idname.compare("") != 0)
- {
- // An ID name has been set, only print the dosage for that ID
- // and SNP combination
- if(idname.compare(fv.readObservationName(col).name) == 0)
- {
- idfound = true;
- cout << fv.readObservationName(col).name << "\t"
- << bufToString(fv.getElementType(),
- &data[col*fv.getElementSize()],
- string("NA"))
- << endl;
- }
- } else {
- // Print the dosages for all IDs
- cout << fv.readObservationName(col).name << "\t"
- << bufToString(fv.getElementType(),
- &data[col*fv.getElementSize()],
- string("NA"))
- << endl;
- }
+ if (idname.compare("") != 0)
+ {
+ // An ID name has been set, only print the dosage for that ID
+ // and SNP combination
+ if (idname.compare(fv.readObservationName(col).name) == 0)
+ {
+ idfound = true;
+ cout << fv.readObservationName(col).name << "\t"
+ << bufToString(fv.getElementType(),
+ &data[col*fv.getElementSize()],
+ string("NA"))
+ << endl;
+ }
+ } else {
+ // Print the dosages for all IDs
+ cout << fv.readObservationName(col).name << "\t"
+ << bufToString(fv.getElementType(),
+ &data[col*fv.getElementSize()],
+ string("NA"))
+ << endl;
+ }
}
if (idfound == false && idname.compare("") != 0)
{
- cerr << "Id '" << idname << "' not found in data file "
- << inputFileName << endl;
- exit(1);
+ cerr << "Id '" << idname << "' not found in data file "
+ << inputFileName << endl;
+ exit(1);
}
delete [] data;
Modified: pkg/ProbABEL/src/gendata.cpp
===================================================================
--- pkg/ProbABEL/src/gendata.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/gendata.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -31,8 +31,11 @@
data[j++] = tmpdata[i];
// std::cout << j << " " << DAG->get_nobservations() << " "
// << nids << "\n";
- } else
+ }
+ else
+ {
report_error("cannot get gendata");
+ }
}
gendata::gendata() : nsnps(0), nids(0), ngpreds(0), DAG(NULL), DAGmask(NULL)
@@ -40,8 +43,10 @@
}
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)
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured,
+ std::string * idnames)
{
nsnps = insnps;
ngpreds = ingpreds;
@@ -49,9 +54,12 @@
DAGmask = new unsigned short int[DAG->getNumObservations()];
if (DAG->getNumObservations() != npeople)
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");
+
long int j = -1;
+
for (unsigned int i = 0; i < npeople; i++)
{
if (allmeasured[i] == 0)
@@ -66,31 +74,34 @@
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());
+ // if (allmeasured[i] && idnames[j] != DAGobsname)
+ // std::cerr << "names do not match for observation at phenofile "
+ // << "line (phe/geno) " << i+1 << "/+1 ("
+ // << idnames[i].c_str() << "/"
+ // << DAGobsname.c_str() << ")\n";
// 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());
+ "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;
// std::cout << "in INI: " << nids << " " << npeople << "\n";
if (nids != nmeasured)
report_error("nids != mneasured (%i != %i)\n", nids, nmeasured);
-
}
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)
+ unsigned int ingpreds, unsigned int npeople,
+ unsigned int nmeasured,
+ unsigned short int * allmeasured, int skipd,
+ std::string * idnames)
{
- nids = nmeasured;
- nsnps = insnps;
+ nids = nmeasured;
+ nsnps = insnps;
ngpreds = ingpreds;
- DAG = NULL;
+ DAG = NULL;
// int nids_all = npeople;
G.reinit(nids, (nsnps * ngpreds));
@@ -116,11 +127,12 @@
// arrow of MaCH/minimac. If found only use the part
// after the arrow as ID.
infile >> tmpstr;
- size_t strpos = tmpstr.find("->") ;
+ size_t strpos = tmpstr.find("->");
if (strpos != string::npos)
{
tmpid = tmpstr.substr(strpos+2, string::npos);
- } else
+ }
+ else
{
tmpid = tmpstr;
}
@@ -150,7 +162,8 @@
// size 8" error messages here. A bug in Valgrind?!
float dosage = strtod(tmpstr.c_str(), (char **) NULL);
G.put(dosage, k, j);
- } else
+ }
+ else
{
std::cerr << "cannot read dose-file: "
<< "check skipd and ngpreds parameters\n";
@@ -159,7 +172,8 @@
}
}
k++;
- } else
+ }
+ else
{
for (int j = 0; j < skipd; j++)
infile >> tmpstr;
@@ -172,7 +186,6 @@
// HERE NEED A NEW CONSTRUCTOR BASED ON DATABELBASECPP OBJECT
gendata::~gendata()
{
-
if (DAG != NULL)
{
delete DAG;
Modified: pkg/ProbABEL/src/gendata.h
===================================================================
--- pkg/ProbABEL/src/gendata.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/gendata.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -13,7 +13,6 @@
#if EIGEN
#include "eigen_mematrix.h"
#include "eigen_mematrix.cpp"
-
#else
#include "mematrix.h"
#endif
Modified: pkg/ProbABEL/src/maskedmatrix.cpp
===================================================================
--- pkg/ProbABEL/src/maskedmatrix.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/maskedmatrix.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -49,17 +49,16 @@
length_of_mask = M.nrow;
//TODO:set type (row,column,symmetric)
//type="symmetric";
-
}
masked_matrix::~masked_matrix()
{
delete[] mask_of_old;
}
+
void masked_matrix::update_mask(short unsigned int *newmask)
{
-
-//find length of masked matrix
+ //find length of masked matrix
int nmeasured = 0;
for (int i = 0; i < length_of_mask; i++)
{
@@ -68,22 +67,25 @@
nmeasured++;
}
}
+
//Check update mask is the same as original matrix
if (nmeasured == length_of_mask)
{
//masked matrix is the same as original matrix
masked_data = &matrix_original;
- } else
+ }
+ else
{
//Check update mask is the same as old matrix
if (is_equal_array(newmask, mask_of_old, length_of_mask))
{
//new mask is the same as old matrix
masked_data = &matrix_masked_data;
- } else
+ }
+ else
{
- //new mask differs from old matrix and create new.
-// //mask_of_old = newmask;
+ // new mask differs from old matrix and create new.
+ // mask_of_old = newmask;
//TODO(maarten): there must be a smarter way to copy these values
for (int i = 0; i < length_of_mask; i++)
{
@@ -97,10 +99,11 @@
void masked_matrix::mask_symmetric(int nmeasured)
{
-//mask a symmetric matrix: this matrix is always a square matrix and will mask
-//rows and columns. The result is always a square matrix
+ // Mask a symmetric matrix: this matrix is always a square matrix and will
+ // mask rows and columns. The result is always a square matrix
matrix_masked_data.reinit(nmeasured, nmeasured);
int i1 = 0, j1 = 0;
+
for (int i = 0; i < length_of_mask; i++)
if (mask_of_old[i] == 0)
{
@@ -117,7 +120,7 @@
}
bool masked_matrix::is_equal_array(unsigned short int *a, unsigned short int *b,
- int size)
+ int size)
{
for (int i = 0; i < size; i++)
{
Modified: pkg/ProbABEL/src/mematri1.h
===================================================================
--- pkg/ProbABEL/src/mematri1.h 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/mematri1.h 2012-11-28 10:58:40 UTC (rev 1032)
@@ -6,7 +6,7 @@
#include <cstdarg>
#include <cstdio>
#include <cstdlib>
-//
+//
// constructors
//
@@ -33,8 +33,8 @@
nrow, ncol);
exit(1);
}
- // fprintf(stderr,"mematrix(nr,nc): can allocate memory (%d,%d)\n",nrow,ncol);
}
+
template<class DT>
mematrix<DT>::mematrix(const mematrix<DT> & M)
{
@@ -53,7 +53,8 @@
for (int i = 0; i < M.ncol * M.nrow; i++)
data[i] = M.data[i];
}
-//
+
+//
// operators
//
template<class DT>
@@ -75,11 +76,13 @@
nrow = M.nrow;
nelements = M.nelements;
for (int i = 0; i < M.ncol * M.nrow; i++)
+ {
data[i] = M.data[i];
- // fprintf(stderr,"mematrix=: can allocate memory (%d,%d)\n",M.nrow,M.ncol);
+ }
}
return *this;
}
+
template<class DT>
DT &mematrix<DT>::operator[](int i)
{
@@ -91,6 +94,7 @@
}
return data[i];
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator+(DT toadd)
{
@@ -99,6 +103,7 @@
temp.data[i] = data[i] + toadd;
return temp;
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator+(mematrix<DT> &M)
{
@@ -114,6 +119,7 @@
temp.data[i] = data[i] + M.data[i];
return temp;
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator-(DT toadd)
{
@@ -122,6 +128,7 @@
temp.data[i] = data[i] - toadd;
return temp;
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator-(mematrix<DT> &M)
{
@@ -137,6 +144,7 @@
temp.data[i] = data[i] - M.data[i];
return temp;
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator*(DT toadd)
{
@@ -146,6 +154,7 @@
temp.data[i] = data[i] * toadd;
return temp;
}
+
template<class DT>
mematrix<DT> mematrix<DT>::operator*(mematrix<DT> &M)
{
@@ -192,7 +201,7 @@
return temp;
}
-//
+//
// operations
//
template<class DT>
@@ -221,6 +230,7 @@
exit(1);
}
}
+
template<class DT>
DT mematrix<DT>::get(int nr, int nc)
{
@@ -239,6 +249,7 @@
DT temp = data[nr * ncol + nc];
return temp;
}
+
template<class DT>
void mematrix<DT>::put(DT value, int nr, int nc)
{
@@ -256,6 +267,7 @@
}
data[nr * ncol + nc] = value;
}
+
template<class DT>
DT mematrix<DT>::column_mean(int nc)
{
@@ -270,6 +282,7 @@
out /= DT(nrow);
return out;
}
+
template<class DT>
void mematrix<DT>::print(void)
{
@@ -283,6 +296,7 @@
cout << "\n";
}
}
+
template<class DT>
void mematrix<DT>::delete_column(int delcol)
{
@@ -313,8 +327,8 @@
if (nc != delcol)
data[nr * ncol + (newcol++)] = temp[nr * temp.ncol + nc];
}
+}
-}
template<class DT>
void mematrix<DT>::delete_row(int delrow)
{
@@ -345,10 +359,9 @@
if (nr != delrow)
data[nr * ncol + (newrow++)] = temp[nr * temp.ncol + nc];
}
-
}
-//
+//
// other functions
//
template<class DT>
@@ -365,6 +378,7 @@
}
return out;
}
+
template<class DT>
mematrix<DT> column_mean(mematrix<DT> &M)
{
@@ -380,6 +394,7 @@
}
return out;
}
+
template<class DT>
mematrix<DT> transpose(mematrix<DT> &M)
{
Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp 2012-11-27 23:59:11 UTC (rev 1031)
+++ pkg/ProbABEL/src/reg1.cpp 2012-11-28 10:58:40 UTC (rev 1032)
@@ -1,7 +1,8 @@
#include "reg1.h"
mematrix<double> apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool is_interaction_excluded, bool iscox, int nullmodel)
+ int ngpreds, bool is_interaction_excluded,
+ bool iscox, int nullmodel)
// model 0 = 2 df
// model 1 = additive 1 df
// model 2 = dominant 1 df
@@ -32,7 +33,8 @@
* X[i * X.ncol + interaction - 1];
nX[i * nX.ncol + c2] = X[i * X.ncol + csnp_p2]
* X[i * X.ncol + interaction - 1];
- } else
+ }
+ else
{
//Maksim: interaction with SNP;;
nX[i * nX.ncol + c1] = X[i * X.ncol + csnp_p1]
@@ -90,7 +92,8 @@
{
nX[i * nX.ncol + c1] = X[i * X.ncol + csnp_p1]
* X[i * X.ncol + interaction - 1]; //Maksim: interaction with SNP;;
- } else
+ }
+ else
{
nX[i * nX.ncol + c1] = X[i * X.ncol + csnp_p1]
* X[i * X.ncol + interaction]; //Maksim: interaction with SNP;;
@@ -111,15 +114,17 @@
{
col_new++;
nX_without_interact_phe[row
- * nX_without_interact_phe.ncol + col_new] =
- nX[row * nX.ncol + col];
+ * nX_without_interact_phe.ncol
+ + col_new] =
+ nX[row * nX.ncol + col];
}
if (col != interaction - 1 && iscox)
{
col_new++;
nX_without_interact_phe[row
- * nX_without_interact_phe.ncol + col_new] =
- nX[row * nX.ncol + col];
+ * nX_without_interact_phe.ncol
+ + col_new] =
+ nX[row * nX.ncol + col];
}
}
}
@@ -128,24 +133,29 @@
//________________________
return (nX);
}
- } else
+ }
+ else
{
return (X);
}
}
+
mematrix<double> nX;
if (interaction != 0)
{
nX.reinit(X.nrow, (X.ncol));
- } else
+ }
+ else
{
nX.reinit(X.nrow, (X.ncol - 1));
}
+
int c1 = X.ncol - 2;
int c2 = X.ncol - 1;
for (int i = 0; i < X.nrow; i++)
for (int j = 0; j < (X.ncol - 2); j++)
nX[i * nX.ncol + j] = X[i * X.ncol + j];
+
for (int i = 0; i < nX.nrow; i++)
{
if (model == 1)
@@ -192,11 +202,11 @@
}
mematrix<double> t_apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool iscox, int nullmodel)
+ int ngpreds, bool iscox, int nullmodel)
{
mematrix<double> tmpX = transpose(X);
mematrix<double> nX = apply_model(tmpX, model, interaction, ngpreds, iscox,
- nullmodel);
+ nullmodel);
mematrix<double> out = transpose(nX);
return out;
}
@@ -222,12 +232,14 @@
}
void base_reg::base_score(mematrix<double>& resid, regdata& rdata, int verbose,
- double tol_chol, int model, int interaction, int ngpreds,
- const masked_matrix& invvarmatrix, int nullmodel)
+ double tol_chol, int model, int interaction,
+ int ngpreds, const masked_matrix& invvarmatrix,
+ int nullmodel)
{
mematrix<double> oX = rdata.extract_genotypes();
mematrix<double> X = apply_model(oX, model, interaction, ngpreds,
- rdata.is_interaction_excluded, false, nullmodel);
+ rdata.is_interaction_excluded, false,
+ nullmodel);
beta.reinit(X.ncol, 1);
sebeta.reinit(X.ncol, 1);
double N = static_cast<double>(resid.nrow);
@@ -259,21 +271,22 @@
double sigma2_internal = (srr - N * mean_r * mean_r) / (N - beta.nrow);
for (int i = 0; i < beta.nrow; i++)
sebeta[i] = sqrt(v_i.get(i, i) * sigma2_internal);
+
mematrix<double> chi2 = transpose(u) * v_i * u;
chi2 = chi2 * (1. / sigma2_internal);
chi2_score = chi2[0];
}
void linear_reg::estimate(regdata& rdatain, int verbose, double tol_chol,
- int model, int interaction, int ngpreds, masked_matrix& invvarmatrixin,
- int robust, int nullmodel)
+ int model, int interaction, int ngpreds,
+ masked_matrix& invvarmatrixin, int robust,
+ int nullmodel)
{
//suda ineraction parameter
// model should come here
regdata rdata = rdatain.get_unmasked_data();
if (invvarmatrixin.length_of_mask != 0)
{
-
invvarmatrixin.update_mask(rdatain.masked_data);
// invvarmatrixin.masked_data->print();
}
@@ -287,7 +300,8 @@
rdata.X.print();
}
mematrix<double> X = apply_model(rdata.X, model, interaction, ngpreds,
- rdata.is_interaction_excluded, false, nullmodel);
+ rdata.is_interaction_excluded, false,
+ nullmodel);
if (verbose)
{
std::cout << "X:\n";
@@ -302,11 +316,13 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1032
More information about the Genabel-commits
mailing list