[GenABEL-dev] [Genabel-commits] r1700 - pkg/ProbABEL/src
L.C. Karssen
lennart at karssen.org
Sun Apr 27 22:30:31 CEST 2014
Hi Maarten,
More clean ups. Great!
Some comments below.
On 27-04-14 11:01, noreply at r-forge.r-project.org wrote:
> Author: maartenk
> Date: 2014-04-27 11:01:42 +0200 (Sun, 27 Apr 2014)
> New Revision: 1700
>
> Removed:
> pkg/ProbABEL/src/cholesky.cpp
> pkg/ProbABEL/src/cholesky.h
> Modified:
> pkg/ProbABEL/src/Makefile.am
> pkg/ProbABEL/src/fvlib
> pkg/ProbABEL/src/reg1.cpp
> pkg/ProbABEL/src/reg1.h
> Log:
> -removed dependency of reg1.* on cholesky.* since this is now done with EIGEN (remove about 150 lines of code from our codebase)
> -removed cholesky.h and cholesky.cpp
> -added some consts to functions in reg1.*
> -removed some whitespace in reg1.cpp
Happy with the consts!
>
> Modified: pkg/ProbABEL/src/fvlib
> ===================================================================
> --- pkg/ProbABEL/src/fvlib 2014-04-25 06:26:38 UTC (rev 1699)
> +++ pkg/ProbABEL/src/fvlib 2014-04-27 09:01:42 UTC (rev 1700)
> @@ -1 +1 @@
> -link ../../../tags/filevector/v.1.0.0/fvlib
> \ No newline at end of file
> +link include/filevector/fvlib
> \ No newline at end of file
This is strange. You seem to have replace the symlink for fvlib to a
place that doesn't exist in the SVN tree. Probably a local thing.
Can you revert this and point the symlink back to the v.1.0.0 tag of fvlib?
>
> Modified: pkg/ProbABEL/src/reg1.cpp
> ===================================================================
> --- pkg/ProbABEL/src/reg1.cpp 2014-04-25 06:26:38 UTC (rev 1699)
> +++ pkg/ProbABEL/src/reg1.cpp 2014-04-27 09:01:42 UTC (rev 1700)
> @@ -275,10 +275,12 @@
> reg_data.is_interaction_excluded, false, nullmodel);
> beta.reinit(X.ncol, 1);
> sebeta.reinit(X.ncol, 1);
> + int length_beta=X.ncol;
Could you please add spaces around the = sign, according to the coding
standards?
> double N = static_cast<double>(resid.nrow);
> mematrix<double> tX = transpose(X);
> - if (invvarmatrix.length_of_mask != 0)
> + if (invvarmatrix.length_of_mask != 0){
> tX = tX * invvarmatrix.masked_data;
> + }
>
> mematrix<double> u = tX * resid;
> mematrix<double> v = tX * X;
> @@ -287,12 +289,16 @@
> csum = csum * (1. / N);
> v = v - csum;
> // use cholesky to invert
> - mematrix<double> v_i = v;
> - cholesky2_mm(v_i, tol_chol);
> - chinv2_mm(v_i);
> +
> + LDLT <MatrixXd> Ch = LDLT < MatrixXd > (v.data.selfadjointView<Lower>());
I get the feeling here that you added too many spaces in this case :-).
The < and > here are not operators.
Thanks,
Lennart.
> // before was
> // mematrix<double> v_i = invert(v);
> - beta = v_i * u;
> + beta.data = Ch.solve(v.data.adjoint() * u.data);
> + //TODO(maartenk): set size of v_i directly or remove mematrix class
> + mematrix<double> v_i = v;
> + v_i.data = Ch.solve(MatrixXd(length_beta, length_beta).
> + Identity(length_beta, length_beta));
> +
> double sr = 0.;
> double srr = 0.;
> for (int i = 0; i < resid.nrow; i++)
> @@ -327,7 +333,7 @@
> sigma2 = (Y - tXW * beta_vec).squaredNorm();
> beta.data = beta_vec;
> }
> -void linear_reg::LeastSquaredRegression(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());
> @@ -368,12 +374,11 @@
> //residuals[i] -= resid_sub;
> loglik -= (residuals.data.array().square() * halfrecsig2).sum();
> loglik -= static_cast<double>(reg_data.nids) * log(sqrt(sigma2));
> -
> }
>
>
>
> -void linear_reg::RobustSEandCovariance(mematrix<double> X, mematrix<double> robust_sigma2,
> +void linear_reg::RobustSEandCovariance(const mematrix<double> &X, mematrix<double> robust_sigma2,
> MatrixXd tXX_inv, int offset) {
> MatrixXd Xresiduals = X.data.array().colwise()
> * residuals.data.col(0).array();
> @@ -386,7 +391,7 @@
> robust_sigma2.data.bottomLeftCorner(offset, offset).diagonal();
> }
>
> -void linear_reg::PlainSEandCovariance(double sigma2_internal, MatrixXd tXX_inv,
> +void linear_reg::PlainSEandCovariance(double sigma2_internal,const MatrixXd &tXX_inv,
> int offset) {
> sebeta.data = (sigma2_internal * tXX_inv.diagonal().array()).sqrt();
> covariance.data = sigma2_internal
> @@ -438,7 +443,6 @@
>
> double sigma2_internal;
>
> -
> LDLT <MatrixXd> Ch;
> if (invvarmatrixin.length_of_mask != 0)
> {
> @@ -481,10 +485,8 @@
>
> MatrixXd tXX_inv = Ch.solve(MatrixXd(length_beta, length_beta).
> Identity(length_beta, length_beta));
> -
> mematrix<double> robust_sigma2(X.ncol, X.ncol);
>
> -
> int offset = X.ncol- 1;
> //if additive and interaction and 2 predictors and more then 2 betas
> if (model == 0 && interaction != 0 && ngpreds == 2 && length_beta > 2){
> @@ -499,10 +501,8 @@
> {
> PlainSEandCovariance(sigma2_internal, tXX_inv, offset);
> }
> -
> }
>
> -
> void linear_reg::score(mematrix<double>& resid,
> double tol_chol, int model, int interaction, int ngpreds,
> const masked_matrix& invvarmatrix, int nullmodel) {
> @@ -511,7 +511,6 @@
> invvarmatrix, nullmodel = 0);
> }
>
> -
> logistic_reg::logistic_reg(regdata& rdatain) {
> reg_data = rdatain.get_unmasked_data();
> int length_beta = reg_data.X.ncol;
>
> Modified: pkg/ProbABEL/src/reg1.h
> ===================================================================
> --- pkg/ProbABEL/src/reg1.h 2014-04-25 06:26:38 UTC (rev 1699)
> +++ pkg/ProbABEL/src/reg1.h 2014-04-27 09:01:42 UTC (rev 1700)
> @@ -50,11 +50,9 @@
> #ifndef REG1_H_
> #define REG1_H_
> #include <cmath>
> -#include "cholesky.h"
> #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);
> @@ -99,11 +97,10 @@
> const masked_matrix& W_masked,
> LDLT<MatrixXd>& Ch);
> void logLikelihood(const mematrix<double>& X);
> - void LeastSquaredRegression(mematrix<double> X, LDLT<MatrixXd>& Ch);
> - void RobustSEandCovariance(mematrix<double> X,
> - mematrix<double> robust_sigma2,
> - MatrixXd tXX_inv, int offset);
> - void PlainSEandCovariance(double sigma2_internal, MatrixXd tXX_inv,
> + 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);
> };
>
>
> _______________________________________________
> Genabel-commits mailing list
> Genabel-commits at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-commits
>
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
L.C. Karssen
Utrecht
The Netherlands
lennart at karssen.org
http://blog.karssen.org
GPG key ID: A88F554A
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 213 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20140427/77d80e62/attachment.sig>
More information about the genabel-devel
mailing list