[Genabel-commits] r1714 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 28 22:47:50 CEST 2014
Author: lckarssen
Date: 2014-04-28 22:47:50 +0200 (Mon, 28 Apr 2014)
New Revision: 1714
Modified:
pkg/ProbABEL/src/reg1.cpp
pkg/ProbABEL/src/reg1.h
Log:
Added consts to the base_reg, linear_reg and logistic_reg classes.
Also includes a little bit of code layout improvements.
The biggest questions that remain are: why is there and assignment nullmodel = 0 in the calls to the score() functions of linear_reg and logistic_reg? Is this a simple copy/paste error?
Modified: pkg/ProbABEL/src/reg1.cpp
===================================================================
--- pkg/ProbABEL/src/reg1.cpp 2014-04-28 20:36:10 UTC (rev 1713)
+++ pkg/ProbABEL/src/reg1.cpp 2014-04-28 20:47:50 UTC (rev 1714)
@@ -23,17 +23,36 @@
#include "reg1.h"
-mematrix<double> apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool is_interaction_excluded, bool iscox, int nullmodel)
-// if ngpreds==1 (dose data):
-// model 0 = additive 1 df
-// if ngpreds==2 (prob data):
-// model 0 = 2 df
-// model 1 = additive 1 df
-// model 2 = dominant 1 df
-// model 3 = recessive 1 df
-// model 4 = over-dominant 1 df
- {
+
+/**
+ *
+ *
+ * if ngpreds==1 (dose data):
+ * model 0 = additive 1 df
+ * if ngpreds==2 (prob data):
+ * model 0 = 2 df
+ * model 1 = additive 1 df
+ * model 2 = dominant 1 df
+ * model 3 = recessive 1 df
+ * model 4 = over-dominant 1 df
+ * @param X Design matrix, including SNP column
+ * @param model
+ * @param interaction
+ * @param ngpreds
+ * @param is_interaction_excluded
+ * @param iscox
+ * @param nullmodel
+ *
+ * @return Matrix with the model applied to it.
+ */
+mematrix<double> apply_model(const mematrix<double>& X,
+ const int model,
+ const int interaction,
+ const int ngpreds,
+ const bool is_interaction_excluded,
+ const bool iscox,
+ const int nullmodel)
+{
if (nullmodel)
{
// No need to apply any genotypic model when calculating the
@@ -239,8 +258,14 @@
return nX;
}
-mematrix<double> t_apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool iscox, int nullmodel) {
+
+mematrix<double> t_apply_model(const mematrix<double>& X,
+ const int model,
+ const int interaction,
+ const int ngpreds,
+ const bool iscox,
+ const int nullmodel)
+{
mematrix<double> tmpX = transpose(X);
mematrix<double> nX = apply_model(tmpX, model, interaction, ngpreds,
interaction, iscox, nullmodel);
@@ -248,7 +273,8 @@
return out;
}
-linear_reg::linear_reg(regdata& rdatain) {
+
+linear_reg::linear_reg(const regdata& rdatain) {
reg_data = rdatain.get_unmasked_data();
// std::cout << "linear_reg: " << rdata.nids << " " << (rdata.X).ncol
// << " " << (rdata.Y).ncol << "\n";
@@ -267,9 +293,13 @@
chi2_score = -1.;
}
-void base_reg::base_score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- const masked_matrix& invvarmatrix, int nullmodel) {
+
+void base_reg::base_score(const mematrix<double>& resid,
+ const double tol_chol, const int model,
+ const int interaction,
+ const int ngpreds,
+ const masked_matrix& invvarmatrix,
+ const int nullmodel) {
mematrix<double> oX = reg_data.extract_genotypes();
mematrix<double> X = apply_model(oX, model, interaction, ngpreds,
reg_data.is_interaction_excluded, false, nullmodel);
@@ -316,8 +346,10 @@
chi2_score = chi2[0];
}
+
void linear_reg::mmscore_regression(const mematrix<double>& X,
- const masked_matrix& W_masked, LDLT<MatrixXd>& Ch) {
+ 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
@@ -333,7 +365,10 @@
sigma2 = (Y - tXW * beta_vec).squaredNorm();
beta.data = beta_vec;
}
-void linear_reg::LeastSquaredRegression(const mematrix<double>& X, LDLT<MatrixXd>& Ch) {
+
+
+void linear_reg::LeastSquaredRegression(const mematrix<double>& X,
+ LDLT<MatrixXd>& Ch) {
int m = X.ncol;
MatrixXd txx = MatrixXd(m, m).setZero().selfadjointView<Lower>().rankUpdate(
X.data.adjoint());
@@ -342,6 +377,7 @@
sigma2 = (reg_data.Y.data - (X.data * beta.data)).squaredNorm();
}
+
void linear_reg::logLikelihood(const mematrix<double>& X) {
/*
loglik = 0.;
@@ -377,9 +413,10 @@
}
-
-void linear_reg::RobustSEandCovariance(const mematrix<double> &X, mematrix<double> robust_sigma2,
- MatrixXd tXX_inv, int offset) {
+void linear_reg::RobustSEandCovariance(const mematrix<double> &X,
+ mematrix<double> robust_sigma2,
+ MatrixXd tXX_inv,
+ const int offset) {
MatrixXd Xresiduals = X.data.array().colwise()
* residuals.data.col(0).array();
MatrixXd XbyR =
@@ -391,16 +428,21 @@
robust_sigma2.data.bottomLeftCorner(offset, offset).diagonal();
}
-void linear_reg::PlainSEandCovariance(double sigma2_internal,const MatrixXd &tXX_inv,
- int offset) {
+
+void linear_reg::PlainSEandCovariance(double sigma2_internal,
+ const MatrixXd &tXX_inv,
+ const int offset) {
sebeta.data = (sigma2_internal * tXX_inv.diagonal().array()).sqrt();
covariance.data = sigma2_internal
* tXX_inv.bottomLeftCorner(offset, offset).diagonal().array();
}
-void linear_reg::estimate(int verbose, double tol_chol,
- int model, int interaction, int ngpreds, masked_matrix& invvarmatrixin,
- int robust, int nullmodel) {
+
+void linear_reg::estimate(const int verbose, const double tol_chol,
+ const int model, const int interaction,
+ const int ngpreds,
+ masked_matrix& invvarmatrixin,
+ const int robust, const int nullmodel) {
// suda interaction parameter
// model should come here
//regdata rdata = rdatain.get_unmasked_data();
@@ -503,15 +545,22 @@
}
}
-void linear_reg::score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- const masked_matrix& invvarmatrix, int nullmodel) {
+
+void linear_reg::score(const mematrix<double>& resid,
+ const double tol_chol, const int model,
+ const int interaction, const int ngpreds,
+ const masked_matrix& invvarmatrix,
+ int nullmodel)
+{
//regdata rdata = rdatain.get_unmasked_data();
- base_score(resid, tol_chol, model, interaction, ngpreds,
- invvarmatrix, nullmodel = 0);
+ base_score(resid, tol_chol, model, interaction, ngpreds,
+ invvarmatrix, nullmodel = 0);
+ // TODO: find out why nullmodel is assigned 0 in the call above.
}
-logistic_reg::logistic_reg(regdata& rdatain) {
+
+logistic_reg::logistic_reg(const regdata& rdatain)
+{
reg_data = rdatain.get_unmasked_data();
int length_beta = reg_data.X.ncol;
beta.reinit(length_beta, 1);
@@ -546,7 +595,8 @@
invvarmatrixin.update_mask(reg_data.masked_data);
}
mematrix<double> X = apply_model(reg_data.X, model, interaction, ngpreds,
- reg_data.is_interaction_excluded, false, nullmodel);
+ reg_data.is_interaction_excluded, false,
+ nullmodel);
int length_beta = X.ncol;
beta.reinit(length_beta, 1);
sebeta.reinit(length_beta, 1);
@@ -762,9 +812,12 @@
}
-void logistic_reg::score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- masked_matrix& invvarmatrix, int nullmodel) {
+void logistic_reg::score(const mematrix<double>& resid,
+ const double tol_chol, const int model,
+ const int interaction, const int ngpreds,
+ const masked_matrix& invvarmatrix,
+ int nullmodel) {
base_score(resid, tol_chol, model, interaction, ngpreds,
- invvarmatrix, nullmodel = 0);
+ invvarmatrix, nullmodel = 0);
+ // TODO: find out why nullmodel is assigned 0 in the call above.
}
Modified: pkg/ProbABEL/src/reg1.h
===================================================================
--- pkg/ProbABEL/src/reg1.h 2014-04-28 20:36:10 UTC (rev 1713)
+++ pkg/ProbABEL/src/reg1.h 2014-04-28 20:47:50 UTC (rev 1714)
@@ -53,13 +53,23 @@
#include "regdata.h"
#include "maskedmatrix.h"
-mematrix<double> apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool is_interaction_excluded, bool iscox = false,
- int nullmodel = 0);
-mematrix<double> t_apply_model(mematrix<double>& X, int model, int interaction,
- int ngpreds, bool iscox, int nullmodel = 0);
+mematrix<double> apply_model(const mematrix<double>& X,
+ const int model,
+ const int interaction,
+ const int ngpreds,
+ const bool is_interaction_excluded,
+ const bool iscox = false,
+ const int nullmodel = 0);
+mematrix<double> t_apply_model(const mematrix<double>& X,
+ const int model,
+ const int interaction,
+ const int ngpreds,
+ const bool iscox,
+ const int nullmodel = 0);
+
+
class base_reg {
public:
mematrix<double> beta;
@@ -73,52 +83,61 @@
double chi2_score;
regdata reg_data;
- void base_score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- const masked_matrix& invvarmatrix, int nullmodel);
+ void base_score(const mematrix<double>& resid,
+ const double tol_chol, const int model,
+ const int interaction, const int ngpreds,
+ const masked_matrix& invvarmatrix,
+ int nullmodel);
};
+
class linear_reg: public base_reg {
public:
- linear_reg(regdata& rdatain);
- // ~linear_reg();
+ linear_reg(const regdata& rdatain);
- void estimate(int verbose, double tol_chol, int model,
- int interaction, int ngpreds,
+ void estimate(const int verbose, const double tol_chol, const int model,
+ const int interaction, const int ngpreds,
masked_matrix& invvarmatrixin,
- int robust, int nullmodel = 0);
+ const int robust, const int nullmodel = 0);
- void score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- const masked_matrix& invvarmatrix, int nullmodel = 0);
+ void score(const mematrix<double>& resid,
+ const double tol_chol, const int model, const int interaction,
+ const int ngpreds, const masked_matrix& invvarmatrix,
+ int nullmodel = 0);
private:
void mmscore_regression(const mematrix<double>& X,
const masked_matrix& W_masked,
LDLT<MatrixXd>& Ch);
void logLikelihood(const mematrix<double>& X);
- void LeastSquaredRegression(const mematrix<double> & X,LDLT<MatrixXd>& Ch);
+ void LeastSquaredRegression(const mematrix<double> & X,
+ LDLT<MatrixXd>& Ch);
void RobustSEandCovariance(const mematrix<double> & X,
- mematrix <double> robust_sigma2, MatrixXd tXX_inv, int offset);
- void PlainSEandCovariance(double sigma2_internal, const MatrixXd & tXX_inv,
- int offset);
+ mematrix <double> robust_sigma2,
+ MatrixXd tXX_inv,
+ const int offset);
+ void PlainSEandCovariance(double sigma2_internal,
+ const MatrixXd & tXX_inv,
+ const int offset);
};
+
class logistic_reg: public base_reg {
public:
int niter;
- logistic_reg(regdata& rdatain);
- // ~logistic_reg();
+ logistic_reg(const regdata& rdatain);
- void estimate(int verbose, int maxiter, double eps,
- int model, int interaction, int ngpreds,
- masked_matrix& invvarmatrixin, int robust,
- int nullmodel = 0);
+ void estimate(const int verbose, const int maxiter, const double eps,
+ const int model, const int interaction, const int ngpreds,
+ masked_matrix& invvarmatrixin, const int robust,
+ const int nullmodel = 0);
+
// just a stupid copy from linear_reg
- void score(mematrix<double>& resid,
- double tol_chol, int model, int interaction, int ngpreds,
- masked_matrix& invvarmatrix, int nullmodel = 0);
+ void score(const mematrix<double>& resid,
+ const double tol_chol, const int model, const int interaction,
+ const int ngpreds, const masked_matrix& invvarmatrix,
+ int nullmodel = 0);
};
#endif // REG1_H_
More information about the Genabel-commits
mailing list