[Genabel-commits] r1695 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 23 20:54:26 CEST 2014
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
More information about the Genabel-commits
mailing list