[GenABEL-dev] [Genabel-commits] r1695 - pkg/ProbABEL/src
L.C. Karssen
lennart at karssen.org
Thu Apr 24 11:45:28 CEST 2014
Wow, that's a lot of code removal!
Glad to see this clean-up. Thanks Maarten! Now we can slowly start
thinking of using Eigen directly instead of going through the
eigen_mematrix "wrapper". However, I don't think this has a high
priority now. Moreover, it is quite a large task, so definitely should
be done in a branch and extensively tested.
Lennart.
On 23-04-14 20:54, noreply at r-forge.r-project.org wrote:
> Author: maartenk
> Date: 2014-04-23 20:54:26 +0200 (Wed, 23 Apr 2014)
> New Revision: 1695
>
> Removed:
> pkg/ProbABEL/src/mematri1.h
> pkg/ProbABEL/src/mematrix.h
> Modified:
> pkg/ProbABEL/src/cholesky.cpp
> pkg/ProbABEL/src/cholesky.h
> pkg/ProbABEL/src/command_line_settings.cpp
> pkg/ProbABEL/src/coxph_data.cpp
> pkg/ProbABEL/src/coxph_data.h
> pkg/ProbABEL/src/data.cpp
> pkg/ProbABEL/src/gendata.cpp
> pkg/ProbABEL/src/gendata.h
> pkg/ProbABEL/src/main.cpp
> pkg/ProbABEL/src/maskedmatrix.cpp
> pkg/ProbABEL/src/maskedmatrix.h
> pkg/ProbABEL/src/phedata.h
> pkg/ProbABEL/src/reg1.cpp
> pkg/ProbABEL/src/regdata.h
> pkg/ProbABEL/src/testchol.cpp
> pkg/ProbABEL/src/usage.cpp
> Log:
> Removed mematri1.h mematrix.h specific code. This remove about 797 lines of code.
>
> Modified: pkg/ProbABEL/src/cholesky.cpp
> ===================================================================
> --- pkg/ProbABEL/src/cholesky.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/cholesky.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -9,13 +9,8 @@
> #include <cstdio>
> #include <cstdlib>
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
>
>
> /* SCCS @(#)cholesky2.c 5.2 10/27/98
>
> Modified: pkg/ProbABEL/src/cholesky.h
> ===================================================================
> --- pkg/ProbABEL/src/cholesky.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/cholesky.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -8,12 +8,8 @@
> #ifndef CHOLESKY_H_
> #define CHOLESKY_H_
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#endif
>
> int cholesky2_mm(mematrix<double> &matrix, double toler);
> void chinv2_mm(mematrix<double> &matrix);
>
> Modified: pkg/ProbABEL/src/command_line_settings.cpp
> ===================================================================
> --- pkg/ProbABEL/src/command_line_settings.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/command_line_settings.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -31,9 +31,7 @@
> #include <iostream>
> #include "usage.h"
> #include "command_line_settings.h"
> -#if EIGEN
> #include "eigen_mematrix.h"
> -#endif
>
> // config.h and fvlib/FileVector.h are included for the upper case variables
> #if HAVE_CONFIG_H
>
> Modified: pkg/ProbABEL/src/coxph_data.cpp
> ===================================================================
> --- pkg/ProbABEL/src/coxph_data.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/coxph_data.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -405,22 +405,14 @@
>
> // When using Eigen coxfit2 needs to be called in a slightly
> // different way (i.e. the .data()-part needs to be added).
> -#if EIGEN
> 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(),
> means.data.data(), beta.data.data(), u.data.data(),
> imat.data.data(), loglik_int, &flag, work, &eps, &tol_chol,
> &sctest);
> -#else
> - coxfit2(&maxiter, &cdata.nids, &X.nrow, cdata.stime.data,
> - cdata.sstat.data, X.data, newoffset.data,
> - cdata.weights.data, cdata.strata.data,
> - means.data, beta.data, u.data,
> - imat.data, loglik_int, &flag, work, &eps, &tol_chol,
> - &sctest);
> -#endif
>
> +
> niter = maxiter;
>
> // Check the results of the Cox fit; mirrored from the same checks
> @@ -449,7 +441,6 @@
> << " setting beta and se to 'nan'\n";
> setToZero = true;
> } else {
> -#if EIGEN
> VectorXd ueigen = u.data;
> MatrixXd imateigen = imat.data;
> VectorXd infs = ueigen.transpose() * imateigen;
> @@ -463,12 +454,7 @@
>
> setToZero = true;
> }
> -#else
> - cerr << "Warning for " << snpinfo.name[cursnp]
> - << ": can't check for infinite betas."
> - << " Please compile ProbABEL with Eigen support to fix this."
> - << endl;
> -#endif
> +
> }
>
> for (int i = 0; i < X.nrow; i++)
>
> Modified: pkg/ProbABEL/src/coxph_data.h
> ===================================================================
> --- pkg/ProbABEL/src/coxph_data.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/coxph_data.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -29,13 +29,8 @@
> #ifndef COXPH_DATA_H_
> #define COXPH_DATA_H_
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
>
> #include "data.h"
> #include "reg1.h"
>
> Modified: pkg/ProbABEL/src/data.cpp
> ===================================================================
> --- pkg/ProbABEL/src/data.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/data.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -38,13 +38,8 @@
> #include "gendata.h"
> #include "data.h"
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
> #include "utilities.h"
>
>
>
> Modified: pkg/ProbABEL/src/gendata.cpp
> ===================================================================
> --- pkg/ProbABEL/src/gendata.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/gendata.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -31,13 +31,8 @@
> #include <limits>
> #include "gendata.h"
> #include "fvlib/FileVector.h"
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
> #include "utilities.h"
>
>
>
> Modified: pkg/ProbABEL/src/gendata.h
> ===================================================================
> --- pkg/ProbABEL/src/gendata.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/gendata.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -31,13 +31,10 @@
> #include <string>
> #include "fvlib/FileVector.h"
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#endif
>
> +
> class gendata {
> public:
> unsigned int nsnps;
>
> Modified: pkg/ProbABEL/src/main.cpp
> ===================================================================
> --- pkg/ProbABEL/src/main.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/main.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -67,15 +67,8 @@
>
> #include <ctime> //needed for timing loading non file vector format
>
> -
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
> -
> #include "maskedmatrix.h"
> #include "data.h"
> #include "reg1.h"
>
> Modified: pkg/ProbABEL/src/maskedmatrix.cpp
> ===================================================================
> --- pkg/ProbABEL/src/maskedmatrix.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/maskedmatrix.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -30,13 +30,8 @@
>
> #include <algorithm>
> #include "maskedmatrix.h"
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
>
> masked_matrix::masked_matrix()
> {
>
> Modified: pkg/ProbABEL/src/maskedmatrix.h
> ===================================================================
> --- pkg/ProbABEL/src/maskedmatrix.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/maskedmatrix.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -29,13 +29,8 @@
> #ifndef MASKEDMATRIX_H_
> #define MASKEDMATRIX_H_
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
>
> class masked_matrix {
> public:
>
> Deleted: pkg/ProbABEL/src/mematri1.h
> ===================================================================
> --- pkg/ProbABEL/src/mematri1.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/mematri1.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -1,636 +0,0 @@
> -/*
> - *
> - * Copyright (C) 2009--2014 Various members of the GenABEL team. See
> - * the SVN commit logs for more details.
> - *
> - * This program is free software; you can redistribute it and/or
> - * modify it under the terms of the GNU General Public License
> - * as published by the Free Software Foundation; either version 2
> - * of the License, or (at your option) any later version.
> - *
> - * This program is distributed in the hope that it will be useful,
> - * but WITHOUT ANY WARRANTY; without even the implied warranty of
> - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> - * GNU General Public License for more details.
> - *
> - * You should have received a copy of the GNU General Public License
> - * along with this program; if not, write to the Free Software
> - * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> - * MA 02110-1301, USA.
> - *
> - */
> -
> -
> -#ifndef MEMATRI1_H
> -#define MEMATRI1_H
> -
> -#include <cstdlib>
> -#include <string>
> -#include <cstdarg>
> -#include <cstdio>
> -
> -//
> -// constructors
> -//
> -
> -template<class DT>
> -mematrix<DT>::mematrix(int nr, int nc)
> -{
> - if (nr <= 0)
> - {
> - fprintf(stderr, "mematrix(): nr <= 0\n");
> - exit(1);
> - }
> - if (nc <= 0)
> - {
> - fprintf(stderr, "mematrix(): nc <= 0\n");
> - exit(1);
> - }
> - nrow = nr;
> - ncol = nc;
> - nelements = nr * nc;
> - data = new (nothrow) DT[ncol * nrow];
> - if (!data)
> - {
> - fprintf(stderr, "mematrix(nr,nc): cannot allocate memory (%d,%d)\n",
> - nrow, ncol);
> - exit(1);
> - }
> -}
> -
> -template<class DT>
> -mematrix<DT>::mematrix(const mematrix<DT> & M)
> -{
> - ncol = M.ncol;
> - nrow = M.nrow;
> - nelements = M.nelements;
> - data = new (nothrow) DT[M.ncol * M.nrow];
> - if (!data)
> - {
> - fprintf(stderr,
> - "mematrix const(mematrix): cannot allocate memory (%d,%d)\n",
> - M.nrow, M.ncol);
> - exit(1);
> - }
> - // std::cerr << "mematrix const(mematrix): can allocate memory ("
> - // << M.nrow << "," << M.ncol << ")\n";
> - for (int i = 0; i < M.ncol * M.nrow; i++)
> - data[i] = M.data[i];
> -}
> -
> -//
> -// operators
> -//
> -template<class DT>
> -mematrix<DT> &mematrix<DT>::operator=(const mematrix<DT> &M)
> -{
> - if (this != &M)
> - {
> - if (data != NULL)
> - delete[] data;
> - data = new (nothrow) DT[M.ncol * M.nrow];
> - if (!data)
> - {
> - fprintf(stderr, "mematrix=: cannot allocate memory (%d,%d)\n",
> - M.nrow, M.ncol);
> - delete[] data;
> - exit(1);
> - }
> - ncol = M.ncol;
> - nrow = M.nrow;
> - nelements = M.nelements;
> - for (int i = 0; i < M.ncol * M.nrow; i++)
> - {
> - data[i] = M.data[i];
> - }
> - }
> - return *this;
> -}
> -
> -template<class DT>
> -DT &mematrix<DT>::operator[](int i)
> -{
> - if (i < 0 || i >= (ncol * nrow))
> - {
> - fprintf(stderr, "mematrix[]: %d out of bounds (0,%d)\n", i,
> - nrow * ncol - 1);
> - exit(1);
> - }
> - return data[i];
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator+(DT toadd)
> -{
> - mematrix<DT> temp(nrow, ncol);
> - for (int i = 0; i < nelements; i++)
> - temp.data[i] = data[i] + toadd;
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator+(mematrix<DT> &M)
> -{
> - if (ncol != M.ncol || nrow != M.nrow)
> - {
> - fprintf(stderr,
> - "mematrix+: matrices not equal in size (%d,%d) and (%d,%d)",
> - nrow, ncol, M.nrow, M.ncol);
> - exit(1);
> - }
> - mematrix<DT> temp(nrow, ncol);
> - for (int i = 0; i < nelements; i++)
> - temp.data[i] = data[i] + M.data[i];
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator-(DT toadd)
> -{
> - mematrix<DT> temp(nrow, ncol);
> - for (int i = 0; i < nelements; i++)
> - temp.data[i] = data[i] - toadd;
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator-(mematrix<DT> &M)
> -{
> - if (ncol != M.ncol || nrow != M.nrow)
> - {
> - fprintf(stderr,
> - "mematrix-: matrices not equal in size (%d,%d) and (%d,%d)",
> - nrow, ncol, M.nrow, M.ncol);
> - exit(1);
> - }
> - mematrix<DT> temp(nrow, ncol);
> - for (int i = 0; i < nelements; i++)
> - temp.data[i] = data[i] - M.data[i];
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator*(DT toadd)
> -{
> - // A che naschet std::string vmesto DT? Maksim.
> - mematrix<DT> temp(nrow, ncol);
> - for (int i = 0; i < nelements; i++)
> - temp.data[i] = data[i] * toadd;
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator*(mematrix<DT> &M)
> -{
> - if (ncol != M.nrow)
> - {
> - fprintf(stderr, "mematrix*: ncol != nrow (%d,%d) and (%d,%d)", nrow,
> - ncol, M.nrow, M.ncol);
> - exit(1);
> - }
> - mematrix<DT> temp(nrow, M.ncol);
> - for (int j = 0; j < temp.nrow; j++)
> - {
> - for (int i = 0; i < temp.ncol; i++)
> - {
> - DT sum = 0;
> - for (int j1 = 0; j1 < ncol; j1++)
> - sum += data[j * ncol + j1] * M.data[j1 * M.ncol + i];
> - temp[j * temp.ncol + i] = sum;
> - }
> - }
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> mematrix<DT>::operator*(mematrix<DT> *M)
> -{
> - if (ncol != M->nrow)
> - {
> - fprintf(stderr, "mematrix*: ncol != nrow (%d,%d) and (%d,%d)", nrow,
> - ncol, M->nrow, M->ncol);
> - exit(1);
> - }
> - mematrix<DT> temp(nrow, M->ncol);
> - for (int j = 0; j < temp.nrow; j++)
> - {
> - for (int i = 0; i < temp.ncol; i++)
> - {
> - DT sum = 0;
> - for (int j1 = 0; j1 < ncol; j1++)
> - sum += data[j * ncol + j1] * M->data[j1 * M->ncol + i];
> - temp[j * temp.ncol + i] = sum;
> - }
> - }
> - return temp;
> -}
> -
> -//
> -// operations
> -//
> -template<class DT>
> -void mematrix<DT>::reinit(int nr, int nc)
> -{
> - if (nelements > 0)
> - delete[] data;
> - if (nr <= 0)
> - {
> - fprintf(stderr, "mematrix(): number of rows smaller then 1\n");
> - exit(1);
> - }
> - if (nc <= 0)
> - {
> - fprintf(stderr, "mematrix(): number of columns smaller then 1\n");
> - exit(1);
> - }
> - nrow = nr;
> - ncol = nc;
> - nelements = nr * nc;
> - data = new (nothrow) DT[ncol * nrow];
> - if (!data)
> - {
> - fprintf(stderr, "mematrix(nr,nc): cannot allocate memory (%d,%d)\n",
> - nrow, ncol);
> - exit(1);
> - }
> -}
> -
> -template<class DT>
> -DT mematrix<DT>::get(int nr, int nc)
> -{
> - if (nc < 0 || nc > ncol -1)
> - {
> - std::cerr << "mematrix::get: column out of range: " << nc + 1
> - << " not between (1," << ncol << ")\n" << std::flush;
> - exit(1);
> - }
> - if (nr < 0 || nr > nrow -1)
> - {
> - std::cerr << "mematrix::get: row out of range: " << nr + 1
> - << " not between (1," << nrow << ")\n" << std::flush;
> - exit(1);
> - }
> - DT temp = data[nr * ncol + nc];
> - return temp;
> -}
> -
> -template<class DT>
> -void mematrix<DT>::put(DT value, int nr, int nc)
> -{
> - if (nc < 0 || nc > ncol -1)
> - {
> - std::cerr << "mematrix::put: column out of range: " << nc + 1
> - << " not between (1," << ncol << ")\n" << std::flush;
> - exit(1);
> - }
> - if (nr < 0 || nr > nrow -1)
> - {
> - std::cerr << "mematrix::put: row out of range: " << nr + 1
> - << " not between (1," << nrow << ")\n" << std::flush;
> - exit(1);
> - }
> - data[nr * ncol + nc] = value;
> -}
> -
> -template<class DT>
> -DT mematrix<DT>::column_mean(int nc)
> -{
> - if (nc >= ncol || nc < 0)
> - {
> - fprintf(stderr, "colmM bad column\n");
> - exit(1);
> - }
> - DT out = 0.0;
> - for (int i = 0; i < nrow; i++)
> - out += DT(data[i * ncol + nc]);
> - out /= DT(nrow);
> - return out;
> -}
> -
> -template<class DT>
> -void mematrix<DT>::print(void)
> -{
> - cout << "nrow=" << nrow << "; ncol=" << ncol << "; nelements=" << nelements
> - << "\n";
> - for (int i = 0; i < nrow; i++)
> - {
> - cout << "nr=" << i << ":\t";
> - for (int j = 0; j < ncol; j++)
> - {
> - printf("%e\t", data[i * ncol + j]);
> - }
> - cout << "\n";
> - }
> -}
> -
> -template<class DT>
> -void mematrix<DT>::delete_column(const int delcol)
> -{
> - if (delcol > ncol || delcol < 0)
> - {
> - fprintf(stderr, "mematrix::delete_column: column out of range\n");
> - exit(1);
> - }
> - mematrix<DT> temp = *this;
> - if (nelements > 0)
> - delete[] data;
> - ncol--;
> - nelements = ncol * nrow;
> - data = new (nothrow) DT[ncol * nrow];
> - if (!data)
> - {
> - fprintf(stderr,
> - "mematrix::delete_column: cannot allocate memory (%d,%d)\n",
> - nrow, ncol);
> - delete[] data;
> - exit(1);
> - }
> - int newcol = 0;
> - for (int nr = 0; nr < temp.nrow; nr++)
> - {
> - newcol = 0;
> - for (int nc = 0; nc < temp.ncol; nc++)
> - if (nc != delcol)
> - data[nr * ncol + (newcol++)] = temp[nr * temp.ncol + nc];
> - }
> -}
> -
> -template<class DT>
> -void mematrix<DT>::delete_row(const int delrow)
> -{
> - if (delrow > nrow || delrow < 0)
> - {
> - fprintf(stderr, "mematrix::delete_row: row out of range\n");
> - exit(1);
> - }
> - mematrix<DT> temp = *this;
> - if (nelements > 0)
> - delete[] data;
> - nrow--;
> - nelements = ncol * nrow;
> - data = new (nothrow) DT[ncol * nrow];
> - if (!data)
> - {
> - fprintf(stderr,
> - "mematrix::delete_row: cannot allocate memory (%d,%d)\n", nrow,
> - ncol);
> - delete[] data;
> - exit(1);
> - }
> - int newrow = 0;
> - for (int nc = 0; nc < temp.ncol; nc++)
> - {
> - newrow = 0;
> - for (int nr = 0; nr < temp.nrow; nr++)
> - if (nr != delrow)
> - data[nr * ncol + (newrow++)] = temp[nr * temp.ncol + nc];
> - }
> -}
> -
> -//
> -// other functions
> -//
> -template<class DT>
> -mematrix<DT> column_sum(mematrix<DT> &M)
> -{
> - mematrix<DT> out;
> - out.reinit(1, M.ncol);
> - for (int j = 0; j < M.ncol; j++)
> - {
> - DT sum = 0;
> - for (int i = 0; i < M.nrow; i++)
> - sum = sum + DT(M.data[i * M.ncol + j]);
> - out.put(sum, 0, j);
> - }
> - return out;
> -}
> -
> -template<class DT>
> -mematrix<DT> column_mean(mematrix<DT> &M)
> -{
> - mematrix<DT> out;
> - out.reinit(1, M.ncol);
> - for (int j = 0; j < M.ncol; j++)
> - {
> - DT sum = 0;
> - for (int i = 0; i < M.nrow; i++)
> - sum = sum + DT(M.data[i * M.ncol + j]);
> - sum /= DT(M.nrow);
> - out.put(sum, 0, j);
> - }
> - return out;
> -}
> -
> -template<class DT>
> -mematrix<DT> transpose(mematrix<DT> &M)
> -{
> - mematrix<DT> temp(M.ncol, M.nrow);
> - for (int i = 0; i < temp.nrow; i++)
> - for (int j = 0; j < temp.ncol; j++)
> - temp.data[i * temp.ncol + j] = M.data[j * M.ncol + i];
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> reorder(mematrix<DT> &M, mematrix<int> order)
> -{
> - if (M.nrow != order.nrow)
> - {
> - std::cerr << "reorder: M & order have different # of rows\n";
> - exit(1);
> - }
> - mematrix<DT> temp(M.nrow, M.ncol);
> - for (int i = 0; i < temp.nrow; i++)
> - for (int j = 0; j < temp.ncol; j++)
> - temp.data[order[i] * temp.ncol + j] = M.data[i * M.ncol + j];
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> productMatrDiag(mematrix<DT> &M, mematrix<DT> &D)
> -{
> - //multiply all rows of M by value of first row of D
> - if (M.ncol != D.nrow)
> - {
> - fprintf(stderr, "productMatrDiag: wrong dimenstions");
> - exit(1);
> - }
> - mematrix<DT> temp(M.nrow, M.ncol);
> - for (int i = 0; i < temp.nrow; i++){
> - for (int j = 0; j < temp.ncol; j++){
> - temp.data[i * temp.ncol + j] = M.data[i * M.ncol + j] * D.data[j];
> - // temp.put(M.get(i,j)*D.get(j,0),i,j);
> - }
> - }
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<double> todouble(mematrix<DT> &M)
> -{
> - mematrix<double> temp(M.nrow, M.ncol);
> - for (int i = 0; i < temp.nelements; i++)
> - temp.data[i] = double(M.data[i]);
> - return temp;
> -}
> -
> -template<class DT>
> -mematrix<DT> productXbySymM(mematrix<DT> &X, mematrix<DT> &M)
> -{
> - if (M.ncol < 1 || M.nrow < 1 || X.ncol < 1 || X.nrow < 1)
> - {
> - fprintf(stderr,
> - "productXbySymM: M.ncol<1 || M.nrow<1 || X.ncol<1 || X.nrow < 1\n");
> - exit(1);
> - }
> - if (M.ncol != M.nrow)
> - {
> - fprintf(stderr, "productXbySymM: M.ncol != M.nrow\n");
> - exit(1);
> - }
> - if (M.ncol != X.ncol)
> - {
> - fprintf(stderr, "productXbySymM: M.ncol != X.ncol\n");
> - exit(1);
> - }
> - if (M.ncol != X.ncol)
> - {
> - fprintf(stderr, "productXbySymM: M.ncol != X.ncol\n");
> - exit(1);
> - }
> -
> - mematrix<DT> out(X.nrow, X.ncol);
> - int i, j, k;
> -
> - double temp1, temp2, value1, value2; // not good should be of <DT>!
> - for (k = 0; k < X.nrow; k++)
> - {
> - temp1 = 0.;
> - for (i = 0; i < X.ncol; i++)
> - {
> - temp1 = X.get(k, i);
> - temp2 = 0.;
> - for (j = (i + 1); j < X.ncol; j++)
> - {
> - value1 = out.get(k, j) + temp1 * M.get(i, j);
> - out.put(value1, k, j);
> - temp2 += M.get(i, j) * X.get(k, j);
> - }
> - value2 = out.get(k, i) + temp2 + M.get(i, i) * X.get(k, i);
> - out.put(value2, k, i);
> - }
> - }
> -
> - return out;
> -}
> -
> -// written by Mike Dinolfo 12/98
> -// modified Yurii Aulchenko 2008-04-22
> -template<class DT>
> -mematrix<DT> invert(mematrix<DT> &M)
> -{
> - if (M.ncol != M.nrow)
> - {
> - fprintf(stderr, "invert: only square matrices possible\n");
> - exit(1);
> - }
> - if (M.ncol == 1)
> - {
> - mematrix<DT> temp(1, 1);
> - temp[0] = 1. / M[0];
> - }
> - /*
> - for (int i=0;i<M.ncol;i++)
> - if (M.data[i*M.ncol+i]==0)
> - {
> - fprintf(stderr,"invert: zero elements in diagonal\n");
> - mematrix<DT> temp = M;
> - for (int i = 0; i < M.ncol; i++)
> - for (int j = 0; j < M.ncol; j++)
> - temp.put(NAN,i,j);
> - return temp;
> - //exit(1);
> - }
> - */
> - int actualsize = M.ncol;
> - int maxsize = M.ncol;
> - mematrix<DT> temp = M;
> - for (int i = 1; i < actualsize; i++)
> - temp.data[i] /= temp.data[0]; // normalize row 0
> - for (int i = 1; i < actualsize; i++)
> - {
> - for (int j = i; j < actualsize; j++)
> - { // do a column of L
> - DT sum = 0.0;
> - for (int k = 0; k < i; k++)
> - sum += temp.data[j * maxsize + k] * temp.data[k * maxsize + i];
> - temp.data[j * maxsize + i] -= sum;
> - }
> - if (i == actualsize - 1)
> - continue;
> - for (int j = i + 1; j < actualsize; j++)
> - { // do a row of U
> - DT sum = 0.0;
> - for (int k = 0; k < i; k++)
> - sum += temp.data[i * maxsize + k] * temp.data[k * maxsize + j];
> - temp.data[i * maxsize + j] = (temp.data[i * maxsize + j] - sum)
> - / temp.data[i * maxsize + i];
> - }
> - }
> - for (int i = 0; i < actualsize; i++) // invert L
> - for (int j = i; j < actualsize; j++)
> - {
> - DT x = 1.0;
> - if (i != j)
> - {
> - x = 0.0;
> - for (int k = i; k < j; k++)
> - x -= temp.data[j * maxsize + k]
> - * temp.data[k * maxsize + i];
> - }
> - temp.data[j * maxsize + i] = x / temp.data[j * maxsize + j];
> - }
> - for (int i = 0; i < actualsize; i++) // invert U
> - for (int j = i; j < actualsize; j++)
> - {
> - if (i == j)
> - continue;
> - DT sum = 0.0;
> - for (int k = i; k < j; k++)
> - sum += temp.data[k * maxsize + j]
> - * ((i == k) ? 1.0 : temp.data[i * maxsize + k]);
> - temp.data[i * maxsize + j] = -sum;
> - }
> - for (int i = 0; i < actualsize; i++) // final inversion
> - for (int j = 0; j < actualsize; j++)
> - {
> - DT sum = 0.0;
> - for (int k = ((i > j) ? i : j); k < actualsize; k++)
> - sum += ((j == k) ? 1.0 : temp.data[j * maxsize + k])
> - * temp.data[k * maxsize + i];
> - temp.data[j * maxsize + i] = sum;
> - }
> - return temp;
> -}
> -
> -//_________Maksim____________
> -template<class DT>
> -DT var(mematrix<DT> &M)
> -{
> - DT sum = 0;
> - for (int i = 0; i < M.nelements; i++)
> - {
> - sum += M.data[i];
> - }
> - DT mean = sum / M.nelements;
> -
> - DT sum2 = 0;
> - for (int i = 0; i < M.nelements; i++)
> - {
> - sum2 += pow(M.data[i] - mean, 2);
> - }
> -
> - return sum2 / (M.nelements - 1);
> -}
> -//_________Maksim____________
> -#endif /* MEMATRI1_H */
>
> Deleted: pkg/ProbABEL/src/mematrix.h
> ===================================================================
> --- pkg/ProbABEL/src/mematrix.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/mematrix.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -1,82 +0,0 @@
> -/*
> - *
> - * Copyright (C) 2009--2014 Various members of the GenABEL team. See
> - * the SVN commit logs for more details.
> - *
> - * This program is free software; you can redistribute it and/or
> - * modify it under the terms of the GNU General Public License
> - * as published by the Free Software Foundation; either version 2
> - * of the License, or (at your option) any later version.
> - *
> - * This program is distributed in the hope that it will be useful,
> - * but WITHOUT ANY WARRANTY; without even the implied warranty of
> - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
> - * GNU General Public License for more details.
> - *
> - * You should have received a copy of the GNU General Public License
> - * along with this program; if not, write to the Free Software
> - * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
> - * MA 02110-1301, USA.
> - *
> - */
> -
> -
> -#ifndef __MEMATRIX_H__
> -#define __MEMATRIX_H__
> -#include <iostream>
> -using namespace std;
> -
> -template<class DT> class mematrix
> -{
> - public:
> - int nrow;
> - int ncol;
> - int nelements;
> - DT * data;
> -
> - mematrix()
> - {
> - nrow = ncol = nelements = 0;
> - data = NULL;
> - }
> - mematrix(int nr, int nc);
> - mematrix(const mematrix &M);
> - ~mematrix()
> - {
> - if (nelements > 0)
> - delete[] data;
> - }
> -
> - mematrix & operator=(const mematrix &M);
> - DT & operator[](int i);
> - mematrix operator+(DT toadd);
> - mematrix operator+(mematrix &M);
> - mematrix operator-(DT toadd);
> - mematrix operator-(mematrix &M);
> - mematrix operator*(DT toadd);
> - mematrix operator*(mematrix &M);
> - mematrix operator*(mematrix *M);
> -
> - void reinit(int nr, int nc);
> -
> - unsigned int getnrow(void)
> - {
> - return nrow;
> - }
> - unsigned int getncol(void)
> - {
> - return ncol;
> - }
> - DT get(int nr, int nc);
> - void put(DT value, int nr, int nc);
> - DT column_mean(int nc);
> - void print(void);
> - void delete_column(const int delcol);
> - void delete_row(const int delrow);
> -
> -};
> -
> -// mematrix transpose(mematrix M);
> -// mematrix invert(mematrix M);
> -
> -#endif
>
> Modified: pkg/ProbABEL/src/phedata.h
> ===================================================================
> --- pkg/ProbABEL/src/phedata.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/phedata.h 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -29,13 +29,8 @@
> #ifndef PHEDATA_H_
> #define PHEDATA_H_
>
> -#if EIGEN
> #include "eigen_mematrix.h"
> #include "eigen_mematrix.cpp"
> -#else
> -#include "mematrix.h"
> -#include "mematri1.h"
> -#endif
>
> class phedata {
> public:
>
> Modified: pkg/ProbABEL/src/reg1.cpp
> ===================================================================
> --- pkg/ProbABEL/src/reg1.cpp 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/reg1.cpp 2014-04-23 18:54:26 UTC (rev 1695)
> @@ -310,7 +310,6 @@
> 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);
> @@ -329,7 +328,6 @@
> beta.data = beta_vec;
> }
>
> -
> void linear_reg::logLikelihood(const mematrix<double>& X) {
> /*
> loglik = 0.;
> @@ -348,10 +346,8 @@
> //cout << endl;
> loglik = 0.;
> double halfrecsig2 = .5 / sigma2;
> -#if EIGEN
> //loglik -= halfrecsig2 * residuals[i] * residuals[i];
>
> -
> double intercept = beta.get(0, 0);
> residuals.data = reg_data.Y.data.array() - intercept;
> //matrix.
> @@ -364,17 +360,7 @@
> //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++)
> - {
> - double resid = reg_data.Y[i] - beta.get(0, 0); // intercept
> - for (int j = 1; j < beta.nrow; j++){
> - resid -= beta.get(j, 0) * X.get(i, j);
> - }
> - residuals[i] = resid;
> - loglik -= halfrecsig2 * resid * resid;
> - }
> -#endif
> +
> }
>
>
> @@ -423,12 +409,8 @@
>
> double sigma2_internal;
>
> -#if EIGEN
>
> LDLT <MatrixXd> Ch;
> -#else
> - mematrix<double> tXX_i;
> -#endif
> if (invvarmatrixin.length_of_mask != 0)
> {
> //retrieve masked data W
> @@ -440,26 +422,7 @@
> //flops=mp(2n-1) (when n is big enough flops=mpn2)
> //Oct 26, 2009
>
> -#if EIGEN
> mmscore_regression(X, invvarmatrixin, Ch);
> -#else
> - // next line is 5997000 flops
> - mematrix<double> tXW = transpose(X) * invvarmatrixin.masked_data;
> - tXX_i = tXW * X; // 17991 flops
> - // use cholesky to invert
> - cholesky2_mm(tXX_i, tol_chol);
> - chinv2_mm(tXX_i);
> - beta = tXX_i * (tXW * reg_data.Y); // flops 15+5997
> - // now compute residual variance
> - sigma2 = 0.;
> - //next line is: 1000+5000+= 6000 flops
> - 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);
> - sigma2 += val * val; // flops: 3000 (iterations counted)
> - }
> -#endif
> double N = X.nrow;
> //sigma2_internal = sigma2 / (N - static_cast<double>(length_beta));
> // Ugly fix to the fact that if we do mmscore, sigma2 is already
> @@ -470,7 +433,7 @@
> }
> else // NO mm-score regression : normal least square regression
> {
> -#if EIGEN
> +
> int m = X.ncol;
> MatrixXd txx = MatrixXd(m, m).setZero().selfadjointView<Lower>().\
> rankUpdate(X.data.adjoint());
> @@ -478,23 +441,7 @@
> beta.data = Ch.solve(X.data.adjoint() * reg_data.Y.data);
> sigma2 = (reg_data.Y.data - (X.data * beta.data)).squaredNorm();
>
> -#else
> - mematrix<double> tX = transpose(X);
> - // use cholesky to invert
> - tXX_i = tX * X;
> - cholesky2_mm(tXX_i, tol_chol);
> - chinv2_mm(tXX_i);
> - beta = tXX_i * (tX * (reg_data.Y));
>
> - // now compute residual variance
> - sigma2 = 0.;
> - mematrix<double> sigma2_matrix = reg_data.Y - (X * beta);
> - for (int i = 0; i < sigma2_matrix.nrow; i++)
> - {
> - double val = sigma2_matrix.get(i, 0);
> - sigma2 += val * val;
> - }
> -#endif
> double N = static_cast<double>(X.nrow);
> double P = static_cast<double>(length_beta);
> sigma2_internal = sigma2 / (N - P);
> @@ -517,38 +464,19 @@
> //cout << endl;
> logLikelihood(X);
>
> -#if EIGEN
> MatrixXd tXX_inv = Ch.solve(MatrixXd(length_beta, length_beta).
> Identity(length_beta, length_beta));
> -#endif
>
> mematrix<double> robust_sigma2(X.ncol, X.ncol);
> if (robust)
> {
> -#if EIGEN
> MatrixXd Xresiduals = X.data.array().colwise()\
> *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;
> -#else
> -
> - mematrix<double> XbyR = X;
> - for (int i = 0; i < X.nrow; i++){
> - for (int j = 0; j < X.ncol; j++)
> - {
> - double tmpval = XbyR.get(i, j) * residuals[i];
> - XbyR.put(tmpval, i, j);
> - }
> - }
> - XbyR = transpose(XbyR) * XbyR;
> - robust_sigma2 = tXX_i * XbyR;
> - robust_sigma2 = robust_sigma2 * tXX_i;
> -
> -#endif
> }
> //cout << "estimate 0\n";
> -#if EIGEN
> if (robust)
> {
> sebeta.data = robust_sigma2.data.diagonal().array().sqrt();
> @@ -578,63 +506,6 @@
> offset).diagonal().array();
> }
>
> -#else
> -
> - //cout << "estimate 0\n";
> - for (int i = 0; i < (length_beta); i++)
> - {
> - if (robust)
> - {
> - // cout << "estimate :robust\n";
> - double value = sqrt(robust_sigma2.get(i, i));
> - sebeta.put(value, i, 0);
> - //Han Chen
> - if (i > 0)
> - {
> - if (model == 0 && interaction != 0 && ngpreds == 2
> - && length_beta > 2)
> - {
> - if (i > 1)
> - {
> - double covval = robust_sigma2.get(i, i - 2);
> - covariance.put(covval, i - 2, 0);
> - }
> - }
> - else
> - {
> - double covval = robust_sigma2.get(i, i - 1);
> - covariance.put(covval, i - 1, 0);
> - }
> - }
> - //Oct 26, 2009
> - }
> - else
> - {
> - // cout << "estimate :non-robust\n";
> - double value = sqrt(sigma2_internal * tXX_i.get(i, i));
> - sebeta.put(value, i, 0);
> - //Han Chen
> - if (i > 0)
> - {
> - if (model == 0 && interaction != 0 && ngpreds == 2
> - && length_beta > 2)
> - {
> - if (i > 1)
> - {
> - double covval = sigma2_internal * tXX_i.get(i, i - 2);
> - covariance.put(covval, i - 2, 0);
> - }
> - }
> - else
> - {
> - double covval = sigma2_internal * tXX_i.get(i, i - 1);
> - covariance.put(covval, i - 1, 0);
> - }
> - }
> - //Oct 26, 2009
> - }
> - }
> -#endif
> }
>
>
>
> Modified: pkg/ProbABEL/src/regdata.h
> ===================================================================
> --- pkg/ProbABEL/src/regdata.h 2014-04-23 09:52:41 UTC (rev 1694)
> +++ pkg/ProbABEL/src/regdata.h 2014-04-23 18:54:26 UTC (rev 1695)
> [TRUNCATED]
>
> To get the complete diff run:
> svnlook diff /svnroot/genabel -r 1695
> _______________________________________________
> 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/20140424/a13eec3f/attachment-0001.sig>
More information about the genabel-devel
mailing list