[Genabel-commits] r880 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 1 14:50:56 CEST 2012
Author: maartenk
Date: 2012-04-01 14:50:56 +0200 (Sun, 01 Apr 2012)
New Revision: 880
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp
branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h
branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/data.h
branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.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/phedata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/phedata.h
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
branches/ProbABEL-refactoring/ProbABEL/src/regdata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/regdata.h
branches/ProbABEL-refactoring/ProbABEL/src/usage.cpp
branches/ProbABEL-refactoring/ProbABEL/src/utilities.cpp
Log:
convert to BSD code style
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-04-01 12:50:56 UTC (rev 880)
@@ -1,4 +1,3 @@
-
/*
* cholesky.cpp
*
@@ -12,7 +11,6 @@
#include <cstdio>
#include <cstdlib>
-
/* SCCS @(#)cholesky2.c 5.2 10/27/98
** subroutine to do Cholesky decompostion on a matrix: C = FDF'
** where F is lower triangular with 1's on the diagonal, and D is diagonal
@@ -34,10 +32,11 @@
** Terry Therneau
*/
-
// modified Yurii Aulchenko 2008-05-20
-int cholesky2_mm(mematrix<double> &matrix, double toler) {
- if (matrix.ncol != matrix.nrow) {
+int cholesky2_mm(mematrix<double> &matrix, double toler)
+{
+ if (matrix.ncol != matrix.nrow)
+ {
fprintf(stderr, "cholesky2_mm: matrix should be square\n");
exit(1);
}
@@ -50,7 +49,8 @@
nonneg = 1;
eps = 0;
- for (i = 0; i < n; i++) {
+ for (i = 0; i < n; i++)
+ {
if (matrix[i * n + i] > eps)
eps = matrix[i * n + i];
for (j = (i + 1); j < n; j++)
@@ -59,15 +59,20 @@
eps *= toler;
rank = 0;
- for (i = 0; i < n; i++) {
+ for (i = 0; i < n; i++)
+ {
pivot = matrix[i * n + i];
- if (pivot < eps) {
+ if (pivot < eps)
+ {
matrix[i * n + i] = 0;
if (pivot < -8 * eps)
nonneg = -1;
- } else {
+ }
+ else
+ {
rank++;
- for (j = (i + 1); j < n; j++) {
+ for (j = (i + 1); j < n; j++)
+ {
temp = matrix[j * n + i] / pivot;
matrix[j * n + i] = temp;
matrix[j * n + j] -= temp * temp * pivot;
@@ -79,8 +84,10 @@
return (rank * nonneg);
}
-void chinv2_mm(mematrix<double> &matrix) {
- if (matrix.ncol != matrix.nrow) {
+void chinv2_mm(mematrix<double> &matrix)
+{
+ if (matrix.ncol != matrix.nrow)
+ {
fprintf(stderr, "cholesky2_mm: matrix should be square\n");
exit(1);
}
@@ -93,10 +100,13 @@
** invert the cholesky in the lower triangle
** take full advantage of the cholesky's diagonal of 1's
*/
- for (i = 0; i < n; i++) {
- if (matrix[i * n + i] > 0) {
+ for (i = 0; i < n; i++)
+ {
+ if (matrix[i * n + i] > 0)
+ {
matrix[i * n + i] = 1 / matrix[i * n + i]; /*this line inverts D */
- for (j = (i + 1); j < n; j++) {
+ for (j = (i + 1); j < n; j++)
+ {
matrix[j * n + i] = -matrix[j * n + i];
for (k = 0; k < i; k++) /*sweep operator */
matrix[j * n + k] += matrix[j * n + i] * matrix[i * n + k];
@@ -109,14 +119,19 @@
** calculate F'DF (inverse of cholesky decomp process) to get inverse
** of original matrix
*/
- for (i = 0; i < n; i++) {
- if (matrix[i * n + i] == 0) { /* singular row */
+ for (i = 0; i < n; i++)
+ {
+ if (matrix[i * n + i] == 0)
+ { /* singular row */
for (j = 0; j < i; j++)
matrix[j * n + i] = 0;
for (j = i; j < n; j++)
matrix[i * n + j] = 0;
- } else {
- for (j = (i + 1); j < n; j++) {
+ }
+ else
+ {
+ for (j = (i + 1); j < n; j++)
+ {
temp = matrix[j * n + i] * matrix[j * n + j];
if (j != i)
matrix[i * n + j] = temp;
@@ -131,4 +146,3 @@
matrix[col * n + row] = matrix[row * n + col];
}
-
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-04-01 12:50:56 UTC (rev 880)
@@ -9,9 +9,6 @@
#define CHOLESKY_H_
#include "mematrix.h"
-
-
-
int cholesky2_mm(mematrix<double> &matrix, double toler);
void chinv2_mm(mematrix<double> &matrix);
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-04-01 12:50:56 UTC (rev 880)
@@ -5,7 +5,6 @@
* Author: mkooyman
*/
-
#include "fvlib/AbstractMatrix.h"
#include "fvlib/CastUtils.h"
#include "fvlib/const.h"
@@ -17,7 +16,8 @@
#include "fvlib/Transposer.h"
// compare for sort of times
-int cmpfun(const void *a, const void *b) {
+int cmpfun(const void *a, const void *b)
+{
double el1 = *(double*) a;
double el2 = *(double*) b;
if (el1 > el2)
@@ -28,7 +28,8 @@
return 0;
}
-coxph_data::coxph_data(const coxph_data &obj) {
+coxph_data::coxph_data(const coxph_data &obj)
+{
nids = obj.nids;
ncov = obj.ncov;
ngpreds = obj.ngpreds;
@@ -43,7 +44,8 @@
for (int i = 0; i < nids; i++)
masked_data[i] = 0;
}
-coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum) {
+coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum)
+{
nids = gend.nids;
masked_data = new unsigned short int[nids];
for (int i = 0; i < nids; i++)
@@ -53,7 +55,8 @@
ncov = phed.ncov + ngpreds;
else
ncov = phed.ncov;
- if (phed.noutcomes != 2) {
+ if (phed.noutcomes != 2)
+ {
fprintf(stderr,
"coxph_data: number of outcomes should be 2 (now: %d)\n",
phed.noutcomes);
@@ -67,11 +70,13 @@
offset.reinit(nids, 1);
strata.reinit(nids, 1);
order.reinit(nids, 1);
- for (int i = 0; i < nids; i++) {
+ for (int i = 0; i < nids; i++)
+ {
// X.put(1.,i,0);
stime[i] = (phed.Y).get(i, 0);
sstat[i] = int((phed.Y).get(i, 1));
- if (sstat[i] != 1 & sstat[i] != 0) {
+ if (sstat[i] != 1 & sstat[i] != 0)
+ {
fprintf(stderr,
"coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",
phed.noutcomes);
@@ -83,7 +88,8 @@
X.put((phed.X).get(i, j), i, j);
if (snpnum > 0)
- for (int j = 0; j < ngpreds; j++) {
+ for (int j = 0; j < ngpreds; j++)
+ {
float snpdata[nids];
gend.get_var(snpnum * ngpreds + j, snpdata);
for (int i = 0; i < nids; i++)
@@ -94,7 +100,8 @@
// 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 i = 0; i < nids; i++)
+ {
weights[i] = 1.0;
offset[i] = 0.0;
strata[i] = 0;
@@ -102,22 +109,26 @@
// sort by time
double tmptime[nids];
int passed_sorted[nids];
- for (int i = 0; i < nids; i++) {
+ for (int i = 0; i < nids; i++)
+ {
tmptime[i] = stime[i];
passed_sorted[i] = 0;
}
qsort(tmptime, nids, sizeof(double), cmpfun);
- for (int i = 0; i < nids; i++) {
+ for (int i = 0; i < nids; i++)
+ {
int passed = 0;
for (int j = 0; j < nids; j++)
if (tmptime[j] == stime[i])
- if (!passed_sorted[j]) {
+ if (!passed_sorted[j])
+ {
order[i] = j;
passed_sorted[j] = 1;
passed = 1;
break;
}
- if (passed != 1) {
+ if (passed != 1)
+ {
fprintf(stderr, "cannot recover element %d\n", i);
exit(1);
}
@@ -135,15 +146,18 @@
// stime.print();
// sstat.print();
}
-void coxph_data::update_snp(gendata &gend, int snpnum) {
+void coxph_data::update_snp(gendata &gend, int snpnum)
+{
// note this sorts by "order"!!!
- for (int j = 0; j < ngpreds; j++) {
+ for (int j = 0; j < ngpreds; j++)
+ {
float snpdata[nids];
for (int i = 0; i < nids; i++)
masked_data[i] = 0;
gend.get_var(snpnum * ngpreds + j, snpdata);
- for (int i = 0; i < nids; i++) {
+ for (int i = 0; i < nids; i++)
+ {
X.put(snpdata[i], (ncov - ngpreds + j), order[i]);
if (isnan(snpdata[i]))
masked_data[order[i]] = 1;
@@ -153,7 +167,8 @@
// for (int j=0;j<ngpreds;j++)
// X.put((gend.G).get(i,(snpnum*ngpreds+j)),(ncov-ngpreds+j),order[i]);
}
-coxph_data::~coxph_data() {
+coxph_data::~coxph_data()
+{
delete[] coxph_data::masked_data;
// delete X;
// delete sstat;
@@ -164,7 +179,8 @@
// delete order;
}
-coxph_data coxph_data::get_unmasked_data() {
+coxph_data coxph_data::get_unmasked_data()
+{
// std::cout << " !!! in get_unmasked_data !!! ";
coxph_data to; // = coxph_data(*this);
// filter missing data
@@ -193,9 +209,11 @@
// std::cout << "(to.X).nrow=" << (to.X).nrow << "\n";
// std::cout << " !!! just before cycle !!! ";
int j = 0;
- for (int i = 0; i < nids; i++) {
+ for (int i = 0; i < nids; i++)
+ {
// std::cout << nids << " " << i << " " << masked_data[i] << "\n";
- if (masked_data[i] == 0) {
+ if (masked_data[i] == 0)
+ {
(to.weights).put(weights.get(i, 1), j, 1);
(to.stime).put(stime.get(i, 1), j, 1);
(to.sstat).put(sstat.get(i, 1), j, 1);
@@ -217,5 +235,3 @@
return (to);
}
-
-
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-04-01 12:50:56 UTC (rev 880)
@@ -8,8 +8,8 @@
#ifndef COXPH_DATA_H_
#define COXPH_DATA_H_
-
-class coxph_data {
+class coxph_data
+{
public:
int nids;
int ncov;
@@ -23,7 +23,8 @@
mematrix<int> order;
unsigned short int * masked_data;
coxph_data get_unmasked_data();
- coxph_data() {
+ coxph_data()
+ {
}
coxph_data(const coxph_data &obj);
coxph_data(phedata &phed, gendata &gend, int snpnum);
@@ -32,6 +33,4 @@
};
-
-
#endif /* COXPH_DATA_H_ */
Modified: branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-04-01 12:50:56 UTC (rev 880)
@@ -22,27 +22,32 @@
extern bool is_interaction_excluded;
-unsigned int Nmeasured(char * fname, int nphenocols, int npeople) {
+unsigned int Nmeasured(char * fname, int nphenocols, int npeople)
+{
int ncov = nphenocols - 2;
int nids_all = npeople;
// first pass -- find unmeasured people
std::ifstream infile(fname);
- if (!infile) {
+ if (!infile)
+ {
std::cerr << "Nmeasured: cannot open file " << fname << endl;
}
char tmp[100];
- for (int i = 0; i < nphenocols; i++) {
+ for (int i = 0; i < nphenocols; i++)
+ {
infile >> tmp;
}
unsigned short int * allmeasured = new unsigned short int[npeople];
int nids = 0;
- for (int i = 0; i < npeople; i++) {
+ for (int i = 0; i < npeople; i++)
+ {
allmeasured[i] = 1;
infile >> tmp;
- for (int j = 1; j < nphenocols; j++) {
+ for (int j = 1; j < nphenocols; j++)
+ {
infile >> tmp;
if (tmp[0] == 'N' || tmp[0] == 'n')
allmeasured[i] = 0;
@@ -57,26 +62,30 @@
return (nids);
}
+mlinfo::mlinfo(char * filename, char * mapname)
+{
-
-mlinfo::mlinfo(char * filename, char * mapname) {
-
char tmp[100];
unsigned int nlin = 0;
std::ifstream infile(filename);
- if (infile.is_open()) {
- while (infile.good()) {
+ if (infile.is_open())
+ {
+ while (infile.good())
+ {
infile >> tmp;
nlin++;
}
nlin--; // Subtract one, the previous loop added 1 too much
- } else {
+ }
+ else
+ {
std::cerr << "mlinfo: cannot open file " << filename << endl;
exit(1);
}
infile.close();
- if (nlin % 7) {
+ if (nlin % 7)
+ {
std::cerr << "mlinfo: number of columns != 7 in " << filename << endl;
exit(1);
}
@@ -92,7 +101,8 @@
map = new std::string[nsnps];
infile.open(filename);
- if (!infile) { // file couldn't be opened
+ if (!infile)
+ { // file couldn't be opened
std::cerr << "mlinfo: cannot open file " << filename << endl;
exit(1);
}
@@ -100,7 +110,8 @@
for (int i = 0; i < 7; i++)
infile >> tmp;
- for (int i = 0; i < nsnps; i++) {
+ for (int i = 0; i < nsnps; i++)
+ {
infile >> tmp;
name[i] = tmp;
infile >> tmp;
@@ -119,16 +130,19 @@
}
infile.close();
- if (mapname != NULL) {
+ if (mapname != NULL)
+ {
std::ifstream instr(mapname);
int BFS = 1000;
char line[BFS], tmp[BFS];
- if (!instr.is_open()) {
+ if (!instr.is_open())
+ {
std::cerr << "mlinfo: cannot open file " << mapname << endl;
exit(1);
}
instr.getline(line, BFS);
- for (int i = 0; i < nsnps; i++) {
+ for (int i = 0; i < nsnps; i++)
+ {
instr.getline(line, BFS);
std::stringstream line_stream(line);
line_stream >> tmp >> map[i];
@@ -136,7 +150,8 @@
instr.close();
}
}
-mlinfo::~mlinfo() {
+mlinfo::~mlinfo()
+{
delete[] mlinfo::name;
delete[] mlinfo::A1;
delete[] mlinfo::A2;
@@ -147,65 +162,75 @@
delete[] mlinfo::map;
}
-
//_________________________________________Maksim_start
-InvSigma::InvSigma(const char * filename_, phedata * phe) {
-filename = filename_;
-npeople = phe->nids;
-std::ifstream myfile(filename_);
-char * line = new char[MAXIMUM_PEOPLE_AMOUNT];
-double val;
-std::string id;
-unsigned row = 0, col = 0;
+InvSigma::InvSigma(const char * filename_, phedata * phe)
+{
+ filename = filename_;
+ npeople = phe->nids;
+ std::ifstream myfile(filename_);
+ char * line = new char[MAXIMUM_PEOPLE_AMOUNT];
+ double val;
+ std::string id;
+ unsigned row = 0, col = 0;
-matrix.reinit(npeople, npeople);
+ matrix.reinit(npeople, npeople);
//idnames[k], if (allmeasured[i]==1)
-if (myfile.is_open()) {
- while (myfile.getline(line, MAXIMUM_PEOPLE_AMOUNT)) {
+ if (myfile.is_open())
+ {
+ while (myfile.getline(line, MAXIMUM_PEOPLE_AMOUNT))
+ {
- std::stringstream line_stream(line);
- line_stream >> id;
+ std::stringstream line_stream(line);
+ line_stream >> id;
- if (phe->idnames[row] != id) {
- 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";
- exit(1);
- }
+ if (phe->idnames[row] != id)
+ {
+ 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";
+ exit(1);
+ }
- while (line_stream >> val) {
- matrix.put(val, row, col);
- col++;
- }
+ while (line_stream >> val)
+ {
+ matrix.put(val, row, col);
+ col++;
+ }
- if (col != npeople) {
- fprintf(stderr,
- "error: inv file: Number of columns in row %d equals to %d but people amount is %d\n",
- row, col, npeople);
- myfile.close();
- exit(1);
+ if (col != npeople)
+ {
+ fprintf(stderr,
+ "error: inv file: Number of columns in row %d equals to %d but people amount is %d\n",
+ row, col, npeople);
+ myfile.close();
+ exit(1);
+ }
+ col = 0;
+ row++;
}
- col = 0;
- row++;
+ myfile.close();
}
- myfile.close();
-} else {
- fprintf(stderr, "error: inv file: cannot open file '%s'\n", filename_);
-}
+ else
+ {
+ fprintf(stderr, "error: inv file: cannot open file '%s'\n", filename_);
+ }
-delete[] line;
+ delete[] line;
}
;
-InvSigma::~InvSigma() {
+InvSigma::~InvSigma()
+{
//af af
}
-mematrix<double> & InvSigma::get_matrix(void) {
-return matrix;
+mematrix<double> & InvSigma::get_matrix(void)
+{
+ return matrix;
}
//________________________________________Maksim_end
Modified: branches/ProbABEL-refactoring/ProbABEL/src/data.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/data.h 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/data.h 2012-04-01 12:50:56 UTC (rev 880)
@@ -14,7 +14,8 @@
#include "phedata.h"
#include "gendata.h"
-class mlinfo {
+class mlinfo
+{
public:
int nsnps;
std::string * name;
@@ -25,13 +26,15 @@
double * Quality;
double * Rsq;
std::string * map;
- mlinfo() {
+ mlinfo()
+ {
}
mlinfo(char * filename, char * mapname);
~mlinfo();
};
-class InvSigma {
+class InvSigma
+{
private:
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-04-01 12:50:56 UTC (rev 880)
@@ -11,11 +11,13 @@
#include "mematri1.h"
#include "utilities.h"
-void gendata::get_var(int var, float * data) {
+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);
- else if (DAG != NULL) {
+ else if (DAG != NULL)
+ {
float tmpdata[DAG->getNumObservations()];
DAG->readVariableAs((unsigned long int) var, tmpdata);
unsigned int j = 0;
@@ -23,17 +25,19 @@
if (!DAGmask[i])
data[j++] = tmpdata[i];
//fprintf(stdout,"%i %i %i\n",j,DAG->get_nobservations(),nids);
- } else
+ }
+ else
report_error("cannot get gendata");
}
-gendata::gendata() {
+gendata::gendata()
+{
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) {
+ int nmeasured, unsigned short int * allmeasured, std::string * idnames)
+{
nsnps = insnps;
ngpreds = ingpreds;
DAG = new FileVector(filename, 128, true);
@@ -43,10 +47,12 @@
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++) {
+ for (unsigned int i = 0; i < npeople; i++)
+ {
if (allmeasured[i] == 0)
DAGmask[i] = 1;
- else {
+ else
+ {
DAGmask[i] = 0;
j++;
}
@@ -73,7 +79,8 @@
}
void gendata::re_gendata(char * fname, int insnps, int ingpreds, int npeople,
int nmeasured, unsigned short int * allmeasured, int skipd,
- std::string * idnames) {
+ std::string * idnames)
+{
nids = nmeasured;
nsnps = insnps;
ngpreds = ingpreds;
@@ -85,7 +92,8 @@
std::ifstream infile;
infile.open(fname);
- if (!infile) {
+ if (!infile)
+ {
std::cerr << "gendata: cannot open file " << fname << endl;
}
@@ -94,8 +102,10 @@
int k = 0;
for (int i = 0; i < npeople; i++)
- if (allmeasured[i] == 1) {
- if (skipd > 0) {
+ if (allmeasured[i] == 1)
+ {
+ if (skipd > 0)
+ {
// int ttt;
char ttt[100];
infile >> tmp;
@@ -104,14 +114,18 @@
// sscanf(tmp,"%s->%s",&ttt, tmpn);
// sscanf(tmp,"%[^->]->%[^->]",&ttt, tmpn);
tmpstr = tmp;
- if (tmpstr.find("->") != string::npos) {
+ if (tmpstr.find("->") != string::npos)
+ {
sscanf(tmp, "%[^->]->%s", ttt, tmpn);
tmpid = tmpn;
- } else {
+ }
+ 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]) {
+ if (tmpid != idnames[k])
+ {
fprintf(stderr,
"phenofile and dosefile did not match at line %d ",
i + 2);
@@ -120,13 +134,18 @@
exit(1);
}
}
- for (int j = 1; j < skipd; j++) {
+ for (int j = 1; j < skipd; j++)
+ {
infile >> tmp;
}
- for (int j = 0; j < (nsnps * ngpreds); j++) {
- if (infile.good()) {
+ for (int j = 0; j < (nsnps * ngpreds); j++)
+ {
+ if (infile.good())
+ {
infile >> tmp;
- } else {
+ }
+ else
+ {
std::cerr
<< "cannot read dose-file: check skipd and ngpreds parameters\n";
infile.close();
@@ -135,7 +154,9 @@
G.put(atof(tmp), k, j);
}
k++;
- } else {
+ }
+ else
+ {
for (int j = 0; j < skipd; j++)
infile >> tmp;
for (int j = 0; j < (nsnps * ngpreds); j++)
@@ -144,9 +165,11 @@
infile.close();
}
// HERE NEED A NEW CONSTRUCTOR BASED ON DATABELBASECPP OBJECT
-gendata::~gendata() {
+gendata::~gendata()
+{
- if (DAG != NULL) {
+ if (DAG != NULL)
+ {
delete DAG;
delete[] DAGmask;
}
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-04-01 12:50:56 UTC (rev 880)
@@ -11,34 +11,30 @@
#include "fvlib/FileVector.h"
#include "mematrix.h"
-
class gendata
{
public:
- int nsnps;
- int nids;
- 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 get_var(int var, float * data);
- ~gendata();
+ int nsnps;
+ int nids;
+ 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 get_var(int var, float * data);
+ ~gendata();
- // MAKE THAT PRIVATE, ACCESS THROUGH GET_SNP
- // ANOTHER PRIVATE OBJECT IS A POINTER TO DATABELBASECPP
- // UPDATE SNP, ALL REGRESSION METHODS: ACCOUNT FOR MISSING
+ // MAKE THAT PRIVATE, ACCESS THROUGH GET_SNP
+ // ANOTHER PRIVATE OBJECT IS A POINTER TO DATABELBASECPP
+ // UPDATE SNP, ALL REGRESSION METHODS: ACCOUNT FOR MISSING
private:
- mematrix <float> G;
- AbstractMatrix * DAG;
- unsigned short int * DAGmask;
- // mematrix<double> G;
+ mematrix<float> G;
+ AbstractMatrix * DAG;
+ unsigned short int * DAGmask;
+ // mematrix<double> G;
};
#endif /* GENDATA_H_ */
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-03-30 22:32:41 UTC (rev 879)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-04-01 12:50:56 UTC (rev 880)
@@ -50,22 +50,36 @@
bool is_interaction_excluded = false; //Oh Holy Matrix, forgive me for this!
-int main(int argc, char * argv[]) {
+int main(int argc, char * argv[])
+{
int next_option;
const char * const short_options = "p:i:d:m:n:c:o:s:t:g:a:erlh:b:vu";
//b - interaction parameter
// ADD --fv FLAG (FILEVECTOR), IN WHICH CASE USE ALTERNATIVE
// CONSTRUCTOR FOR GENDATA
- const struct option long_options[] = { { "pheno", 1, NULL, 'p' }, { "info",
- 1, NULL, 'i' }, { "dose", 1, NULL, 'd' }, { "map", 1, NULL, 'm' }, {
- "nids", 1, NULL, 'n' }, { "chrom", 1, NULL, 'c' }, { "out", 1, NULL,
- 'o' }, { "skipd", 1, NULL, 's' }, { "ntraits", 1, NULL, 't' }, {
- "ngpreds", 1, NULL, 'g' }, { "separat", 1, NULL, 'a' }, { "score",
- 0, NULL, 'r' }, { "no-head", 0, NULL, 'e' }, { "allcov", 0, NULL,
- 'l' }, { "help", 0, NULL, 'h' }, { "interaction", 1, NULL, 'b' }, {
- "interaction_only", 1, NULL, 'k' }, { "mmscore", 1, NULL, 'v' }, {
- "robust", 0, NULL, 'u' }, { NULL, 0, NULL, 0 } };
+ const struct option long_options[] =
+ {
+ { "pheno", 1, NULL, 'p' },
+ { "info", 1, NULL, 'i' },
+ { "dose", 1, NULL, 'd' },
+ { "map", 1, NULL, 'm' },
+ { "nids", 1, NULL, 'n' },
+ { "chrom", 1, NULL, 'c' },
+ { "out", 1, NULL, 'o' },
+ { "skipd", 1, NULL, 's' },
+ { "ntraits", 1, NULL, 't' },
+ { "ngpreds", 1, NULL, 'g' },
+ { "separat", 1, NULL, 'a' },
+ { "score", 0, NULL, 'r' },
+ { "no-head", 0, NULL, 'e' },
+ { "allcov", 0, NULL, 'l' },
+ { "help", 0, NULL, 'h' },
+ { "interaction", 1, NULL, 'b' },
+ { "interaction_only", 1, NULL, 'k' },
+ { "mmscore", 1, NULL, 'v' },
+ { "robust", 0, NULL, 'u' },
+ { NULL, 0, NULL, 0 } };
char * program_name = argv[0];
char *phefilename = NULL;
@@ -83,7 +97,8 @@
int interaction_excluded = 0;
int robust = 0;
string chrom = "-1";
- int neco[] = { 0, 0, 0 };
+ int neco[] =
+ { 0, 0, 0 };
bool iscox = false;
int isFVF = 0;
#if COXPH
@@ -94,10 +109,12 @@
#endif
int skipd = 2;
int allcov = 0;
- do {
+ do
+ {
next_option = getopt_long(argc, argv, short_options, long_options,
NULL);
- switch (next_option) {
+ switch (next_option)
+ {
case 'h':
print_help(program_name, 0);
case 'p':
@@ -171,11 +188,13 @@
"%s v. %s (C) Yurii Aulchenko, Lennart C. Karssen, Maksim Struchalin, EMCR\n\n",
PACKAGE, PACKAGE_VERSION);
- if (neco[0] != 1 || neco[1] != 1 || neco[2] != 1) {
+ if (neco[0] != 1 || neco[1] != 1 || neco[2] != 1)
+ {
print_usage(program_name, 1);
}
- if (score) {
+ if (score)
+ {
cout << "option --score suppressed from v 0.1-6\n";
exit(1);
}
@@ -238,13 +257,15 @@
else
fprintf(stdout, "\t --robust = OFF\n");
- if (ngpreds != 1 && ngpreds != 2) {
+ if (ngpreds != 1 && ngpreds != 2)
+ {
fprintf(stderr,
"\n\n--ngpreds should be 1 for MLDOSE or 2 for MLPROB\n");
exit(1);
}
- if (interaction_excluded != 0) {
+ if (interaction_excluded != 0)
+ {
interaction = interaction_excluded; //ups
is_interaction_excluded = true;
}
@@ -272,8 +293,8 @@
#if COXPH
interaction_cox--;
#endif
- if (interaction < 0 || interaction > phd.ncov
- || interaction_cox > phd.ncov) {
+ if (interaction < 0 || interaction > phd.ncov || interaction_cox > phd.ncov)
+ {
std::cerr << "error: Interaction parameter is out of range (ineraction="
<< interaction << ") \n";
exit(1);
@@ -315,13 +336,15 @@
}
#endif
- if (inverse_filename != NULL) {
+ if (inverse_filename != NULL)
+ {
std::cout << "you are running mmscore...\n";
}
std::cout << "Reading data ...";
- if (inverse_filename != NULL) {
+ if (inverse_filename != NULL)
+ {
InvSigma inv(inverse_filename, &phd);
invvarmatrix = inv.get_matrix();
double par = 1.; //var(phd.Y)*phd.nids/(phd.nids-phd.ncov-1);
@@ -390,21 +413,24 @@
//________________________________________________________________
//Maksim, 9 Jan, 2009
- if (outfilename == NULL) {
+ if (outfilename == NULL)
+ {
outfilename = (char *) string("regression").c_str();
}
std::string outfilename_str(outfilename);
std::vector<std::ofstream*> outfile;
- if (nohead != 1) {
+ if (nohead != 1)
+ {
if (ngpreds == 2) //All models output. One file per each model
- {
+ {
// open a file for output
//_____________________
- for (int i = 0; i < 5; i++) {
+ for (int i = 0; i < 5; i++)
+ {
outfile.push_back(new std::ofstream());
}
@@ -414,27 +440,32 @@
outfile[3]->open((outfilename_str + "_recess.out.txt").c_str());
outfile[4]->open((outfilename_str + "_over_domin.out.txt").c_str());
- if (!outfile[0]->is_open()) {
+ if (!outfile[0]->is_open())
+ {
std::cerr << "Cannot open file for writing: "
<< outfilename_str + "_2df.out.txt" << "\n";
exit(1);
}
- if (!outfile[1]->is_open()) {
+ if (!outfile[1]->is_open())
+ {
std::cerr << "Cannot open file for writing: "
<< outfilename_str + "_add.out.txt" << "\n";
exit(1);
}
- if (!outfile[2]->is_open()) {
+ if (!outfile[2]->is_open())
+ {
std::cerr << "Cannot open file for writing: "
<< outfilename_str + "_domin.out.txt" << "\n";
exit(1);
}
- if (!outfile[3]->is_open()) {
+ if (!outfile[3]->is_open())
+ {
std::cerr << "Cannot open file for writing: "
<< outfilename_str + "_recess.out.txt" << "\n";
exit(1);
}
- if (!outfile[4]->is_open()) {
+ if (!outfile[4]->is_open())
+ {
std::cerr << "Cannot open file for writing: "
<< outfilename_str + "_over_domin.out.txt" << "\n";
exit(1);
@@ -443,7 +474,8 @@
//Header
//_____________________
- for (int i = 0; i < outfile.size(); i++) {
+ for (int i = 0; i < outfile.size(); i++)
+ {
(*outfile[i]) << "name" << sep << "A1" << sep << "A2" << sep
<< "Freq1" << sep << "MAF" << sep << "Quality" << sep
<< "Rsq" << sep << "n" << sep
@@ -470,7 +502,8 @@
*outfile[3] << sep << "beta_SNP_recA1" << sep << "sebeta_SNP_recA1";
*outfile[4] << sep << "beta_SNP_odom" << sep << "sebeta_SNP_odom";
- if (interaction != 0) {
+ if (interaction != 0)
+ {
//Han Chen
*outfile[0] << sep << "beta_SNP_A1A2_"
<< phd.model_terms[interaction_cox] << sep
@@ -480,7 +513,8 @@
<< sep << "sebeta_SNP_A1A1_"
<< phd.model_terms[interaction_cox];
#if !COXPH
- if (inverse_filename == NULL && !allcov) {
+ if (inverse_filename == NULL && !allcov)
+ {
*outfile[0] << sep << "cov_SNP_A1A2_int_SNP_"
<< phd.model_terms[interaction_cox] << sep
<< "cov_SNP_A1A1_int_SNP_"
@@ -488,14 +522,16 @@
}
#endif
//Oct 26, 2009
- for (int file = 1; file < outfile.size(); file++) {
+ for (int file = 1; file < outfile.size(); file++)
+ {
*outfile[file] << sep << "beta_SNP_"
<< phd.model_terms[interaction_cox] << sep
<< "sebeta_SNP_"
<< phd.model_terms[interaction_cox];
//Han Chen
#if !COXPH
- if (inverse_filename == NULL && !allcov) {
+ if (inverse_filename == NULL && !allcov)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 880
More information about the Genabel-commits
mailing list