[Genabel-commits] r1686 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Apr 8 09:28:06 CEST 2014
Author: lckarssen
Date: 2014-04-08 09:28:06 +0200 (Tue, 08 Apr 2014)
New Revision: 1686
Modified:
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/gendata.cpp
pkg/ProbABEL/src/gendata.h
pkg/ProbABEL/src/reg1.cpp
pkg/ProbABEL/src/reg1.h
pkg/ProbABEL/src/regdata.cpp
Log:
Fixing some easy-to-fix cpplint errors found by Jenkins after the merging of the ProbABEL 0.5.0 branch.
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/coxph_data.cpp 2014-04-08 07:28:06 UTC (rev 1686)
@@ -81,7 +81,7 @@
freq = 0;
masked_data = new unsigned short int[nids];
- std::copy(obj.masked_data, obj.masked_data+nids,masked_data);
+ std::copy(obj.masked_data, obj.masked_data + nids, masked_data);
}
coxph_data::coxph_data(phedata &phed, gendata &gend, const int snpnum)
@@ -92,7 +92,7 @@
masked_data = new unsigned short int[nids];
- std::fill(masked_data,masked_data+nids,0);
+ std::fill(masked_data,masked_data + nids, 0);
ngpreds = gend.ngpreds;
if (snpnum >= 0)
@@ -162,7 +162,7 @@
// sort by time
double *tmptime = new double[nids];
int *passed_sorted = new int[nids];
- std::fill (passed_sorted,passed_sorted+nids,0);
+ std::fill (passed_sorted, passed_sorted + nids, 0);
for (int i = 0; i < nids; i++)
@@ -235,7 +235,7 @@
for (int j = 0; j < ngpreds; j++) {
double *snpdata = new double[nids];
- std::fill (masked_data,masked_data+nids,0);
+ std::fill (masked_data, masked_data + nids, 0);
gend->get_var(snpnum * ngpreds + j, snpdata);
@@ -311,7 +311,7 @@
// filter missing data
- int nmeasured=std::count (masked_data, masked_data+nids, 0);
+ int nmeasured = std::count(masked_data, masked_data + nids, 0);
to.nids = nmeasured;
to.ncov = ncov;
@@ -346,7 +346,7 @@
//delete [] to.masked_data;
to.masked_data = new unsigned short int[to.nids];
- std::fill(to.masked_data,to.masked_data+to.nids,0);
+ std::fill(to.masked_data, to.masked_data + to.nids, 0);
return (to);
}
Modified: pkg/ProbABEL/src/gendata.cpp
===================================================================
--- pkg/ProbABEL/src/gendata.cpp 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/gendata.cpp 2014-04-08 07:28:06 UTC (rev 1686)
@@ -41,27 +41,30 @@
#include "utilities.h"
-void gendata::mldose_line_to_matrix(int k,const char *all_numbers,int amount_of_numbers){
+void gendata::mldose_line_to_matrix(int k, const char *all_numbers,
+ int amount_of_numbers){
int j = 0;
- //check if not a null pointer
+ // Check if not a null pointer
if (!*all_numbers){
perror("Error while reading genetic data (expected pointer to char but found a null pointer)");
exit(EXIT_FAILURE);
}
- while (j<amount_of_numbers)
+
+ while (j < amount_of_numbers)
{
double result = 0;
- //skip whitespace
+ // Skip whitespace
while (*all_numbers == ' ')
{
all_numbers++;
}
- //check NaN (right now checks only first character)
- //TODO: make catching of NaN more rigid
+
+ // check NaN (right now checks only first character)
+ // TODO: make catching of NaN more rigid
if (*all_numbers == 'N')
{
result = std::numeric_limits<double>::quiet_NaN();
- //skip other characters of NaN
+ // Skip other characters of NaN
while ((*all_numbers == 'a') | (*all_numbers == 'N'))
{
all_numbers++;
@@ -70,18 +73,18 @@
else
{
int sign = 0;
- //set sign to -1 if negative: multiply by sign just before return
+ // set sign to -1 if negative: multiply by sign just before return
if (*all_numbers == '-')
{
all_numbers++;
sign = -1;
}
- //read digits before dot
+ // Read digits before dot
while (*all_numbers <= '9' && *all_numbers >= '0')
{
result = result * 10 + (*all_numbers++ - '0');
}
- //read digit after dot
+ // Read digit after dot
if (*all_numbers == '.')
{
double decimal_counter = 1.0;
@@ -92,7 +95,7 @@
result += (*all_numbers++ - '0') * decimal_counter;
}
}
- //correct for negative number
+ // Correct for negative number
if (sign == -1)
{
result = sign * result;
@@ -103,6 +106,7 @@
}
}
+
void gendata::get_var(int var, double * data)
{
// Read the genetic data for SNP 'var' and store in the array 'data'
@@ -163,10 +167,12 @@
}
}
+
gendata::gendata() : nsnps(0), nids(0), ngpreds(0), DAG(NULL), DAGmask(NULL)
{
}
+
void gendata::re_gendata(string filename, unsigned int insnps,
unsigned int ingpreds, unsigned int npeople,
unsigned int nmeasured,
@@ -188,7 +194,9 @@
for (unsigned int i = 0; i < npeople; i++)
{
if (allmeasured[i] == 0)
+ {
DAGmask[i] = 1;
+ }
else
{
DAGmask[i] = 0;
@@ -197,7 +205,9 @@
string DAGobsname = DAG->readObservationName(i).name;
if (DAGobsname.find("->") != string::npos)
+ {
DAGobsname = DAGobsname.substr(DAGobsname.find("->") + 2);
+ }
// if (allmeasured[i] && idnames[j] != DAGobsname)
// std::cerr << "names do not match for observation at phenofile "
@@ -206,10 +216,11 @@
// << 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",
+ "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";
@@ -217,6 +228,7 @@
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,
@@ -335,15 +347,21 @@
else
{
for (int j = 0; j < skipd; j++)
+ {
infile >> tmpstr;
+ }
for (unsigned int j = 0; j < (nsnps * ngpreds); j++)
+ {
infile >> tmpstr;
+ }
}
}
infile.close();
}
+
+
// HERE NEED A NEW CONSTRUCTOR BASED ON DATABELBASECPP OBJECT
gendata::~gendata()
{
Modified: pkg/ProbABEL/src/gendata.h
===================================================================
--- pkg/ProbABEL/src/gendata.h 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/gendata.h 2014-04-08 07:28:06 UTC (rev 1686)
@@ -44,7 +44,8 @@
unsigned int nids;
unsigned int ngpreds;
gendata();
- void mldose_line_to_matrix(int k,const char *all_numbers,int amount_of_numbers);
+ void mldose_line_to_matrix(int k, const char *all_numbers,
+ int amount_of_numbers);
void re_gendata(char * fname, unsigned int insnps, unsigned int ingpreds,
unsigned int npeople, unsigned int nmeasured,
Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/reg1.cpp 2014-04-08 07:28:06 UTC (rev 1686)
@@ -310,9 +310,9 @@
chi2_score = chi2[0];
}
+
void linear_reg::mmscore_regression(const mematrix<double>& X,
const masked_matrix& W_masked, LDLT<MatrixXd>& Ch) {
-
VectorXd Y = reg_data.Y.data.col(0);
/*
in ProbABEL <0.50 this calculation was performed like t(X)*W
@@ -327,9 +327,9 @@
VectorXd beta_vec = Ch.solve(tXW.transpose() * Y);
sigma2 = (Y - tXW * beta_vec).squaredNorm();
beta.data = beta_vec;
-
}
+
void linear_reg::logLikelihood(const mematrix<double>& X) {
/*
loglik = 0.;
@@ -364,7 +364,6 @@
//residuals[i] -= resid_sub;
loglik -= (residuals.data.array().square() * halfrecsig2).sum();
loglik -= static_cast<double>(reg_data.nids) * log(sqrt(sigma2));
-
#else
for (int i = 0; i < reg_data.nids; i++)
{
@@ -378,6 +377,7 @@
#endif
}
+
void linear_reg::estimate(int verbose, double tol_chol,
int model, int interaction, int ngpreds, masked_matrix& invvarmatrixin,
int robust, int nullmodel) {
@@ -453,7 +453,7 @@
// now compute residual variance
sigma2 = 0.;
//next line is: 1000+5000+= 6000 flops
- mematrix<double> sigma2_matrix = reg_data.Y - (transpose(tXW) * beta); //flops: 1000+5000
+ mematrix<double> sigma2_matrix = reg_data.Y - (transpose(tXW) * beta);
for (int i = 0; i < sigma2_matrix.nrow; i++)
{
double val = sigma2_matrix.get(i, 0);
@@ -467,9 +467,8 @@
// YSA, 2009.07.20
sigma2_internal = 1.0;
sigma2 /= N;
-
}
- else//NO mm-score regression : normal least square regression
+ else // NO mm-score regression : normal least square regression
{
#if EIGEN
int m = X.ncol;
@@ -520,7 +519,7 @@
#if EIGEN
MatrixXd tXX_inv = Ch.solve(MatrixXd(length_beta, length_beta).
- Identity(length_beta,length_beta));
+ Identity(length_beta, length_beta));
#endif
mematrix<double> robust_sigma2(X.ncol, X.ncol);
@@ -528,10 +527,10 @@
{
#if EIGEN
MatrixXd Xresiduals = X.data.array().colwise()\
- *residuals.data.col(0).array();
+ *residuals.data.col(0).array();
MatrixXd XbyR = MatrixXd(X.ncol, X.ncol).setZero()\
- .selfadjointView<Lower>().rankUpdate(Xresiduals.adjoint());
- robust_sigma2.data= tXX_inv*XbyR *tXX_inv;
+ .selfadjointView<Lower>().rankUpdate(Xresiduals.adjoint());
+ robust_sigma2.data = tXX_inv * XbyR * tXX_inv;
#else
mematrix<double> XbyR = X;
@@ -547,7 +546,6 @@
robust_sigma2 = robust_sigma2 * tXX_i;
#endif
-
}
//cout << "estimate 0\n";
#if EIGEN
@@ -565,7 +563,7 @@
//if additive and interaction and 2 predictors and more then 2 betas
if (model == 0 && interaction != 0 && ngpreds == 2 && length_beta > 2){
- offset = X.ncol - 2;
+ offset = X.ncol - 2;
}
if (robust)
@@ -639,6 +637,7 @@
#endif
}
+
void linear_reg::score(mematrix<double>& resid,
double tol_chol, int model, int interaction, int ngpreds,
const masked_matrix& invvarmatrix, int nullmodel) {
@@ -647,6 +646,7 @@
invvarmatrix, nullmodel = 0);
}
+
logistic_reg::logistic_reg(regdata& rdatain) {
reg_data = rdatain.get_unmasked_data();
int length_beta = reg_data.X.ncol;
@@ -665,6 +665,7 @@
chi2_score = -1.;
}
+
void logistic_reg::estimate(int verbose, int maxiter,
double eps, int model, int interaction, int ngpreds,
masked_matrix& invvarmatrixin, int robust, int nullmodel) {
@@ -896,6 +897,7 @@
// exit(1);
}
+
void logistic_reg::score(mematrix<double>& resid,
double tol_chol, int model, int interaction, int ngpreds,
masked_matrix& invvarmatrix, int nullmodel) {
Modified: pkg/ProbABEL/src/reg1.h
===================================================================
--- pkg/ProbABEL/src/reg1.h 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/reg1.h 2014-04-08 07:28:06 UTC (rev 1686)
@@ -100,7 +100,7 @@
double tol_chol, int model, int interaction, int ngpreds,
const masked_matrix& invvarmatrix, int nullmodel = 0);
-private:
+ private:
void mmscore_regression(const mematrix<double>& X,
const masked_matrix& W_masked, LDLT<MatrixXd>& Ch);
void logLikelihood(const mematrix<double>& X);
Modified: pkg/ProbABEL/src/regdata.cpp
===================================================================
--- pkg/ProbABEL/src/regdata.cpp 2014-04-07 20:34:25 UTC (rev 1685)
+++ pkg/ProbABEL/src/regdata.cpp 2014-04-08 07:28:06 UTC (rev 1686)
@@ -63,8 +63,7 @@
is_interaction_excluded = obj.is_interaction_excluded;
masked_data = new unsigned short int[nids];
- std::copy(obj.masked_data, obj.masked_data+nids,masked_data);
-
+ std::copy(obj.masked_data, obj.masked_data + nids, masked_data);
}
@@ -76,8 +75,9 @@
nids = gend.nids;
masked_data = new unsigned short int[nids];
- std::fill (masked_data,masked_data+nids,0);
+ std::fill(masked_data, masked_data + nids, 0);
ngpreds = gend.ngpreds;
+
if (snpnum >= 0)
{
ncov = phed.ncov + ngpreds;
@@ -86,9 +86,11 @@
{
ncov = phed.ncov;
}
+
noutcomes = phed.noutcomes;
X.reinit(nids, (ncov + 1));
Y.reinit(nids, noutcomes);
+
for (int i = 0; i < nids; i++)
{
X.put(1., i, 0);
@@ -132,7 +134,7 @@
for (int j = 0; j < ngpreds; j++)
{
double *snpdata = new double[nids];
- std::fill (masked_data,masked_data+nids,0);
+ std::fill(masked_data, masked_data + nids, 0);
gend->get_var(snpnum * ngpreds + j, snpdata);
for (int i = 0; i < nids; i++) {
@@ -198,7 +200,7 @@
regdata regdata::get_unmasked_data()
{
regdata to; // = regdata(*this);
- int nmeasured=std::count (masked_data, masked_data+nids, 0);
+ int nmeasured = std::count(masked_data, masked_data + nids, 0);
to.nids = nmeasured;
to.ncov = ncov;
to.ngpreds = ngpreds;
@@ -230,7 +232,7 @@
// delete [] to.masked_data;
const int arr_size = nids;
to.masked_data = new unsigned short int[arr_size];
- std::copy(masked_data, masked_data+arr_size,to.masked_data);
+ std::copy(masked_data, masked_data + arr_size, to.masked_data);
return (to);
}
More information about the Genabel-commits
mailing list