[Genabel-commits] r906 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 14 21:16:26 CEST 2012
Author: maartenk
Date: 2012-05-14 21:16:26 +0200 (Mon, 14 May 2012)
New Revision: 906
Added:
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp
branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h
branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.h
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/phedata.h
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
branches/ProbABEL-refactoring/ProbABEL/src/regdata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/regdata.h
branches/ProbABEL-refactoring/ProbABEL/src/testchol.cpp
Log:
-added initial use of EIGEN math library for palinear by building with -DEIGEN flag and -DARMA_NO_DEBUG (turn debug of) -I /path/to/eigen/
-commented out coxfix2 in reg1.h: cox regression is broken!
- Fixed bugs introduced by refactoring
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -4,8 +4,14 @@
* Created on: Mar 15, 2012
* Author: mkooyman
*/
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#else
#include "mematrix.h"
#include "mematri1.h"
+#endif
+
#include <string>
#include <cstdarg>
#include <cstdio>
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -7,7 +7,12 @@
#ifndef CHOLESKY_H_
#define CHOLESKY_H_
+
+#if EIGEN
+#include "eigen_mematrix.h"
+#else
#include "mematrix.h"
+#endif
int cholesky2_mm(mematrix<double> &matrix, double toler);
void chinv2_mm(mematrix<double> &matrix);
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -81,10 +81,10 @@
// X.put(1.,i,0);
stime[i] = (phed.Y).get(i, 0);
sstat[i] = int((phed.Y).get(i, 1));
- if (sstat[i] != 1 & sstat[i] != 0)
+ if (sstat[i] != 1 && sstat[i] != 0)
{
fprintf(stderr,
- "coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",
+ "coxph_data: status not 0/1 (right order: id, fuptime, status ...) %d \n",
phed.noutcomes);
exit(1);
}
Modified: branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -15,8 +15,13 @@
#include "gendata.h"
#include "data.h"
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#else
#include "mematrix.h"
#include "mematri1.h"
+#endif
#include "utilities.h"
using namespace std;
Added: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h (rev 0)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -0,0 +1,338 @@
+#ifndef EIGEN_MEMATRI1_H
+#define EIGEN_MEMATRI1_H
+#include "eigen_mematrix.h"
+#include <cstdlib>
+#include <string>
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <Eigen/Dense>
+#include <Eigen/LU>
+
+using namespace Eigen;
+//
+// 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;
+ Matrix<DT, Dynamic, Dynamic,RowMajor> data;
+
+
+}
+template<class DT>
+mematrix<DT>::mematrix(const mematrix<DT> & M)
+{
+ ncol = M.ncol;
+ nrow = M.nrow;
+ nelements = M.nelements;
+ data = M.data;
+
+}
+//
+// operators
+//
+template<class DT>
+mematrix<DT> &mematrix<DT>::operator=(const mematrix<DT> &M)
+{
+ if (this != &M)
+ {
+
+ ncol = M.ncol;
+ nrow = M.nrow;
+ nelements = M.nelements;
+
+ data = M.data;
+ // fprintf(stderr,"mematrix=: can allocate memory (%d,%d)\n",M.nrow,M.ncol);
+ }
+ return *this;
+}
+
+template<class DT>
+DT & mematrix<DT>::operator[](const int i)
+{
+ if (i < 0 || i >= (ncol * nrow))
+ {
+ fprintf(stderr, "mematrix[]: %d out of bounds (0,%d)\n", i,
+ nrow * ncol - 1);
+ exit(1);
+ }
+ int column = i % ncol;
+ int row = (int) floor((double) i / ncol);
+
+ return data(row, column);
+}
+
+//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;
+ temp.data = data + M.data;
+
+ return temp;
+}
+
+template<class DT>
+mematrix<DT> mematrix<DT>::operator-(DT toadd)
+{
+ mematrix<DT> temp(nrow, ncol);
+ temp.data = data.array() - 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;
+ temp.data = data - M.data;
+ temp.ncol=temp.data.cols();
+ temp.nrow=temp.data.rows();
+ temp.nelements=temp.nrow* temp.ncol;
+ return temp;
+}
+
+template<class DT>
+mematrix<DT> mematrix<DT>::operator*(DT multiplyer)
+{
+// MatrixXd add = MatrixXd::Constant(nrow, ncol, toadd);
+ mematrix<DT> temp(nrow, ncol);
+ temp.data = data * multiplyer;
+ 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);
+
+ }
+ mematrix<DT> temp;
+ temp.data = data * M.data;
+ temp.ncol=temp.data.cols();
+ temp.nrow=temp.data.rows();
+ temp.nelements=temp.nrow* temp.ncol;
+
+ 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;
+// std::cout << "[DEBUG] resizing" << std::endl;
+
+ data.resize(nr, nc);
+ data.setZero();
+
+}
+template<class DT>
+DT mematrix<DT>::get(int nr, int nc)
+{
+ if (nc < 0 || nc > ncol)
+ {
+ fprintf(stderr,
+ "mematrix::get: column out of range: %d not in (0,%d)\n", nc,
+ ncol);
+ exit(1);
+ }
+ if (nr < 0 || nr > nrow)
+ {
+ printf("mematrix::get: row out of range: %d not in (0,%d)\n", nr, nrow);
+ exit(1);
+ }
+ DT temp = data(nr, nc);
+ return temp;
+}
+template<class DT>
+void mematrix<DT>::put(DT value, int nr, int nc)
+{
+ if (nc < 0 || nc > ncol)
+ {
+ fprintf(stderr,
+ "mematrix::put: column out of range: %d not in (0,%d)\n", nc,
+ ncol);
+ exit(1);
+ }
+ if (nr < 0 || nr > nrow)
+ {
+ printf("mematrix::put: row out of range: %d not in (0,%d)\n", nr, nrow);
+ exit(1);
+ }
+ data(nr, 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);
+ }
+
+ return data.col(nc).mean();;
+}
+
+template<class DT>
+mematrix<DT> column_sum(mematrix<DT> &M)
+{
+ mematrix<DT> out;
+ out.reinit(1, M.ncol);
+ out.data = M.data.colwise().mean();
+ 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++)
+ cout << data.data()[i * ncol + j] << "\t";
+ cout << "\n";
+ }
+
+
+
+
+}
+
+//
+// other functions
+//
+
+template<class DT>
+mematrix<DT> transpose(mematrix<DT> &M)
+{
+// cout << "[DEBUG TRANSPOSE PRE]nrow=" << M.nrow << "; ncol=" << M.ncol << "; nelements=" << M.nelements;
+
+ mematrix<DT> temp = M;
+ temp.data.transposeInPlace();
+ int swap=temp.ncol;
+ temp.ncol=temp.nrow;
+ temp.nrow=swap;
+// cout << "[DEBUG TRANSPOSE post]nrow=" << temp.nrow << "; ncol=" << temp.ncol << "; nelements=" << temp.nelements;
+
+ return temp;
+}
+
+//template<class DT>
+//mematrix<DT> reorder(mematrix<DT> &M, mematrix<int> order)
+//{
+// if (M.nrow != order.nrow)
+// {
+// fprintf(stderr, "reorder: M & order have differet # 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<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> invert(mematrix<DT> &M)
+{
+ if (M.ncol != M.nrow)
+ {
+ fprintf(stderr, "invert: only square matrices possible\n");
+ exit(1);
+ }
+
+ mematrix<DT> temp = M;
+
+ temp.data = temp.data.inverse();
+ 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;
+// mematrix<DT> temp(M.nrow, M.ncol);
+
+ for (int i = 0; i < temp.nrow; i++)
+ {
+ temp.data.row(i) = M.data.row(i).array() * D.data.row(0).array();
+ }
+ return temp;
+}
+
+#endif
Property changes on: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h (rev 0)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -0,0 +1,64 @@
+#ifndef __EIGEN_MEMATRIX_H__
+#define __EIGEN_MEMATRIX_H__
+#include <iostream>
+#include <Eigen/Dense>
+
+using namespace Eigen;
+
+using namespace std;
+
+template<class DT> class mematrix
+{
+public:
+ int nrow;
+ int ncol;
+ int nelements;
+ Matrix<DT, Dynamic, Dynamic,RowMajor> data;
+
+ mematrix()
+ {
+ nrow = ncol = nelements = 0;
+// std::cout << "[DEBUG] creating object" << std::endl;
+
+ data.resize(1, 1);
+ }
+ 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);
+
+ 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);
+ DT column_sum(int nc);
+ void print(void);
+
+
+};
+
+//mematrix transpose(mematrix M);
+//mematrix invert(mematrix M);
+
+#endif
Property changes on: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -7,8 +7,13 @@
#include <string>
#include "gendata.h"
#include "fvlib/FileVector.h"
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#else
#include "mematrix.h"
#include "mematri1.h"
+#endif
#include "utilities.h"
void gendata::get_var(int var, float * data)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -9,7 +9,12 @@
#define GENDATA_H_
#include <string>
#include "fvlib/FileVector.h"
+
+#if EIGEN
+#include "eigen_mematrix.h"
+#else
#include "mematrix.h"
+#endif
class gendata
{
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -35,8 +35,13 @@
#include <stdio.h>
#include <vector>
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#else
#include "mematrix.h"
#include "mematri1.h"
+#endif
#include "data.h"
#include "reg1.h"
#include "comand_line_settings.h"
@@ -133,7 +138,7 @@
<< input_var.getSep() << "n" << input_var.getSep()
<< "Mean_predictor_allele";
if (input_var.getChrom() != "-1")
- (*outfile[i]) << input_var.getSep() << "input_var.getChrom()";
+ (*outfile[i]) << input_var.getSep() << "chrom";
if (input_var.getMapfilename() != NULL)
(*outfile[i]) << input_var.getSep() << "position";
}
@@ -321,7 +326,9 @@
#elif LINEAR
linear_reg nrd = linear_reg(nrgd);
-
+#if DEBUG
+ std::cout << "[DEBUG] linear_reg nrd = linear_reg(nrgd); DONE.";
+#endif
nrd.estimate(nrgd, 0, CHOLTOL, 0, input_var.getInteraction(),
input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif COXPH
@@ -707,20 +714,80 @@
rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
#elif LINEAR
//cout << (rgd.get_unmasked_data()).nids << " 1\n";
+#if DEBUG
+ rgd.X.print();
+ rgd.Y.print();
+#endif
linear_reg rd(rgd);
+#if DEBUG
+ rgd.X.print();
+ rgd.Y.print();
+#endif
//cout << (rgd.get_unmasked_data()).nids << " 2\n";
if (input_var.getScore())
+ {
+#if DEBUG
+ cout << "input_var.getScore/n";
+ nrd.residuals.print();
+ cout << CHOLTOL << " <-CHOLTOL\n";
+ cout << model << " <-model\n";
+ cout << input_var.getInteraction()
+ << " <-input_var.getInteraction()\n";
+ cout << input_var.getNgpreds()
+ << " <-input_var.getNgpreds()\n";
+ invvarmatrix.print();
+#endif
rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
input_var.getInteraction(), input_var.getNgpreds(),
invvarmatrix);
+#if DEBUG
+ rd.beta.print();
+ cout << rd.chi2_score << " <-chi2_scoren\n";
+ rd.covariance.print();
+ rd.residuals.print();
+ rd.sebeta.print();
+ cout << rd.loglik << " <-logliken\n";
+ cout << rd.sigma2 << " <-sigma2\n";
+#endif
+ }
else
{
// if(input_var.getInverseFilename()== NULL)
// {
//cout << (rgd.get_unmasked_data()).nids << " 3\n";
+#if DEBUG
+ cout << "rd.estimate\n";
+ cout << CHOLTOL << " <-CHOLTOL\n";
+ cout << model << " <-model\n";
+ cout << input_var.getInteraction()
+ << " <-input_var.getInteraction()\n";
+ cout << input_var.getNgpreds()
+ << " <-input_var.getNgpreds()\n";
+ cout << input_var.getRobust()
+ << " <-input_var.getRobust()\n";
+ cout << "start invarmatrix\n";
+ invvarmatrix.print();
+ cout << "end invarmatrix\n";
+ cout << rgd.is_interaction_excluded
+ << " <-rgd.is_interaction_excluded\n";
+#endif
rd.estimate(rgd, 0, CHOLTOL, model,
input_var.getInteraction(), input_var.getNgpreds(),
invvarmatrix, input_var.getRobust());
+
+#if DEBUG
+ cout << "rd.beta\n";
+ rd.beta.print();
+ cout << rd.chi2_score << " <-chi2_scoren\n";
+ cout << "rd.covariance\n";
+ rd.covariance.print();
+ cout << "rd.residuals\n";
+ rd.residuals.print();
+ cout << "rd.sebeta\n";
+ rd.sebeta.print();
+ cout << rd.loglik << " <-logliken\n";
+ cout << rd.sigma2 << " <-sigma2\n";
+#endif
//cout << (rgd.get_unmasked_data()).nids << " 4\n";
// }
// else
@@ -734,13 +801,22 @@
#endif
if (!input_var.getAllcov() && input_var.getInteraction() == 0)
+ {
start_pos = rd.beta.nrow - 1;
+ }
else if (!input_var.getAllcov()
&& input_var.getInteraction() != 0)
+ {
start_pos = rd.beta.nrow - 2;
+ }
else
+ {
start_pos = 0;
+ }
+#if DEBUG
+ cout << " start_pos;" << start_pos << "\n";
+#endif
for (int pos = start_pos; pos < rd.beta.nrow; pos++)
{
*beta_sebeta[0] << input_var.getSep() << rd.beta[pos]
@@ -764,6 +840,9 @@
//________________________________
if (input_var.getInverseFilename() == NULL)
{
+#if DEBUG
+ cout << " inverse_filename == NULL" << "\n";
+#endif
if (input_var.getScore() == 0)
{
*chi2[0] << rd.loglik; //2.*(rd.loglik-null_loglik);
@@ -831,6 +910,9 @@
else
{
*outfile[0] << beta_sebeta[0]->str() << "\n";
+#if DEBUG
+ cout << "Se beta" << beta_sebeta[0] << "\n";
+#endif
}
}
//clean chi2
Modified: branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -385,6 +385,7 @@
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");
Modified: branches/ProbABEL-refactoring/ProbABEL/src/phedata.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/phedata.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/phedata.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -8,7 +8,11 @@
#ifndef PHEDATA_H_
#define PHEDATA_H_
+#if EIGEN
+#include "eigen_mematrix.h"
+#else
#include "mematrix.h"
+#endif
class phedata
{
Modified: branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/reg1.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/reg1.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -40,7 +40,20 @@
// model 2 = dominant 1 df
// model 3 = recessive 1 df
// model 4 = over-dominant 1 df
+
{
+#if DEBUG
+
+ X.print();
+
+ std::cout << "apply_model():model" << model << std::endl;
+ std::cout << "apply_model():interaction" << interaction << std::endl;
+ std::cout << "apply_model():ngpreds" << ngpreds << std::endl;
+ std::cout << "apply_model():is_interaction_excluded"
+ << is_interaction_excluded << std::endl;
+ std::cout << "apply_model():iscox" << iscox << std::endl;
+ std::cout << "apply_model():nullmodel" << nullmodel << std::endl;
+#endif
if (model == 0)
{
if (interaction != 0 && !nullmodel)
@@ -235,6 +248,9 @@
return nX_without_interact_phe;
} //interaction_only, model!=0, ngpreds==2
//Oct 26, 2009
+#if DEBUG
+ nX.print();
+#endif
return nX;
}
@@ -378,7 +394,12 @@
//suda ineraction parameter
// model should come here
regdata rdata = rdatain.get_unmasked_data();
+ if (verbose)
+ {
+ cout << rdata.is_interaction_excluded
+ << " <-irdata.is_interaction_excluded\n";
+ }
mematrix<double> invvarmatrix;
if (invvarmatrixin.nrow != 0 && invvarmatrixin.ncol != 0)
{
@@ -398,11 +419,26 @@
}
}
+ if (verbose)
+ {
+ printf("invvarmatrix:\n");
+ invvarmatrix.print();
+ }
+ if (verbose)
+ {
+ printf("rdata.X:\n");
+ rdata.X.print();
+ }
//fprintf(stdout,"estimate: %i %i %i %i\n",rdata.nids,(rdata.X).nrow,(rdata.X).ncol,(rdata.Y).ncol);
mematrix<double> X = apply_model(rdata.X, model, interaction, ngpreds,
rdata.is_interaction_excluded, false, nullmodel);
+ if (verbose)
+ {
+ printf("X:\n");
+ X.print();
+ }
int length_beta = X.ncol;
beta.reinit(length_beta, 1);
sebeta.reinit(length_beta, 1);
@@ -421,7 +457,7 @@
}
//Oct 26, 2009
mematrix<double> tX = transpose(X);
-
+
if (invvarmatrix.nrow != 0 && invvarmatrix.ncol != 0)
{
tX = tX * invvarmatrix;
@@ -429,6 +465,11 @@
//tX = productXbySymM(tX,invvarmatrix);
// X = invvarmatrix*X; std::cout<<"new tX.nrow="<<X.nrow<<" tX.ncol="<<X.ncol<<"\n";
}
+ if (verbose)
+ {
+ printf("tX:\n");
+ tX.print();
+ }
mematrix<double> tXX = tX * X;
@@ -476,8 +517,15 @@
mematrix<double> ttX = transpose(tX);
mematrix<double> sigma2_matrix = rdata.Y;
mematrix<double> sigma2_matrix1 = ttX * beta;
+// printf("sigma2_matrix");
+// sigma2_matrix.print();
+//
+// printf("sigma2_matrix1");
+// sigma2_matrix1.print();
+
sigma2_matrix = sigma2_matrix - sigma2_matrix1;
-
+// printf("sigma2_matrix");
+// sigma2_matrix.print();
static double val;
// std::cout<<"sigma2_matrix.nrow="<<sigma2_matrix.nrow<<"sigma2_matrix.ncol"<<sigma2_matrix.ncol<<"\n";
@@ -485,7 +533,9 @@
for (int i = 0; i < sigma2_matrix.nrow; i++)
{
val = sigma2_matrix.get(i, 0);
+// printf("val = %f\n", val);
sigma2 += val * val;
+// printf("sigma2+= = %f\n", sigma2);
}
double sigma2_internal = sigma2 / (N - double(length_beta));
@@ -505,6 +555,7 @@
// replaced for ML
// sigma2_internal = sigma2/(N - double(length_beta) - 1);
+// printf("sigma2/=N = %f\n", sigma2);
sigma2 /= N;
// std::cout<<"N="<<N<<", length_beta="<<length_beta<<"\n";
@@ -573,6 +624,7 @@
{
if (robust)
{
+// cout << "estimate :robust\n";
double value = sqrt(robust_sigma2.get(i, i));
sebeta.put(value, i, 0);
//Han Chen
@@ -597,6 +649,7 @@
}
else
{
+// cout << "estimate :non-robust\n";
double value = sqrt(sigma2_internal * tXX_i.get(i, i));
sebeta.put(value, i, 0);
//Han Chen
@@ -1058,11 +1111,11 @@
double loglik_int[2];
int flag;
double sctest = 1.0;
-
- 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);
+//TODO(maarten): remove the following comment signs. This is done for testing purpose only EVIL!
+// 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);
for (int i = 0; i < X.nrow; i++)
sebeta[i] = sqrt(imat.get(i, i));
loglik = loglik_int[1];
Modified: branches/ProbABEL-refactoring/ProbABEL/src/regdata.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/regdata.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/regdata.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -28,6 +28,7 @@
ncov = obj.ncov;
ngpreds = obj.ngpreds;
noutcomes = obj.noutcomes;
+ is_interaction_excluded=obj.is_interaction_excluded;
X = obj.X;
Y = obj.Y;
masked_data = new unsigned short int[nids];
@@ -115,6 +116,7 @@
to.ncov = ncov;
to.ngpreds = ngpreds;
to.noutcomes = noutcomes;
+ to.is_interaction_excluded=is_interaction_excluded;
int dim2Y = Y.ncol;
int dim2X = X.ncol;
(to.X).reinit(to.nids, dim2X);
Modified: branches/ProbABEL-refactoring/ProbABEL/src/regdata.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/regdata.h 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/regdata.h 2012-05-14 19:16:26 UTC (rev 906)
@@ -8,8 +8,13 @@
#ifndef REGDATA_H_
#define REGDATA_H_
+#if EIGEN
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#else
#include "mematrix.h"
#include "mematri1.h"
+#endif
#include "gendata.h"
#include "phedata.h"
Modified: branches/ProbABEL-refactoring/ProbABEL/src/testchol.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/testchol.cpp 2012-05-02 06:48:53 UTC (rev 905)
+++ branches/ProbABEL-refactoring/ProbABEL/src/testchol.cpp 2012-05-14 19:16:26 UTC (rev 906)
@@ -5,8 +5,13 @@
#include <iomanip>
#include <stdio.h>
#include <getopt.h>
+#if EIGEN
#include "mematrix.h"
#include "mematri1.h"
+#else
+#include "eigen_mematrix.h"
+#include "eigen_mematri1.h"
+#endif
#include "data.h"
#include "reg1.h"
#include "usage.h"
More information about the Genabel-commits
mailing list