[Genabel-commits] r1018 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 12 22:54:03 CET 2012
Author: maartenk
Date: 2012-11-12 22:54:03 +0100 (Mon, 12 Nov 2012)
New Revision: 1018
Added:
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.cpp
branches/ProbABEL-refactoring/ProbABEL/src/reg1.cpp
Removed:
branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp
branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/Makefile.am
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/coxph_data.h
branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.h
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/maskedmatrix.cpp
branches/ProbABEL-refactoring/ProbABEL/src/maskedmatrix.h
branches/ProbABEL-refactoring/ProbABEL/src/phedata.h
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
branches/ProbABEL-refactoring/ProbABEL/src/regdata.h
branches/ProbABEL-refactoring/ProbABEL/src/testchol.cpp
Log:
-split reg1.h to reg1.cpp and reg1.h
-merge coxph_reg.* to coxph_reg.* (code is interacting heavely with each other: I see no reason not to merge these files)
- renamed eigen_mematri1.h to eigen_mematrix.cpp
-- ./configure --enable-pacoxph compiles, but example do not run yet
Modified: branches/ProbABEL-refactoring/ProbABEL/src/Makefile.am
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/Makefile.am 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/Makefile.am 2012-11-12 21:54:03 UTC (rev 1018)
@@ -6,9 +6,9 @@
REGFILES = data.h data.cpp gendata.h gendata.cpp mematrix.h mematri1.h \
command_line_settings.h command_line_settings.cpp reg1.h usage.h \
usage.cpp main.cpp utilities.h utilities.cpp phedata.h phedata.cpp \
- cholesky.h cholesky.cpp regdata.h regdata.cpp coxph_data.h \
- coxph_data.cpp maskedmatrix.cpp maskedmatrix.h eigen_mematrix.h \
- eigen_mematri1.h
+ cholesky.h cholesky.cpp regdata.h regdata.cpp \
+ maskedmatrix.cpp maskedmatrix.h eigen_mematrix.h \
+ eigen_mematrix.cpp reg1.cpp
## Filevector files:
FVSRC = fvlib/AbstractMatrix.cpp fvlib/CastUtils.cpp \
@@ -47,7 +47,7 @@
palogist_CPPFLAGS = $(AM_CPPFLAGS)
pacoxph_SOURCES = $(COXSRC) $(REGFILES) $(FVSRC) $(FVHEADERS) \
- $(RHEADERS) survS.h survproto.h
+ $(RHEADERS) survS.h survproto.h coxph_data.h coxph_data.cpp
pacoxph_CXXFLAGS = -DCOXPH -I $(top_srcdir)/src/include $(AM_CXXFLAGS)
pacoxph_CPPFLAGS = $(AM_CPPFLAGS)
pacoxph_CFLAGS = -DCOXPH -I $(top_srcdir)/src/include $(AM_CFLAGS)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.cpp 2012-11-12 21:54:03 UTC (rev 1018)
@@ -6,7 +6,7 @@
*/
#if EIGEN
#include "eigen_mematrix.h"
-#include "eigen_mematri1.h"
+#include "eigen_mematrix.cpp"
#else
#include "mematrix.h"
#include "mematri1.h"
@@ -73,8 +73,7 @@
matrix[i * n + i] = 0;
if (pivot < -8 * eps)
nonneg = -1;
- }
- else
+ } else
{
rank++;
for (j = (i + 1); j < n; j++)
@@ -133,8 +132,7 @@
matrix[j * n + i] = 0;
for (j = i; j < n; j++)
matrix[i * n + j] = 0;
- }
- else
+ } else
{
for (j = (i + 1); j < n; j++)
{
Modified: branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/cholesky.h 2012-11-12 21:54:03 UTC (rev 1018)
@@ -10,6 +10,7 @@
#if EIGEN
#include "eigen_mematrix.h"
+#include "eigen_mematrix.cpp"
#else
#include "mematrix.h"
#endif
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.cpp 2012-11-12 21:54:03 UTC (rev 1018)
@@ -4,8 +4,13 @@
* Created on: Mar 31, 2012
* Author: mkooyman
*/
+#include "coxph_data.h"
+#include <cmath>
+extern "C" {
+#include "survproto.h"
+}
-#include "coxph_data.h"
+// #include "reg1.h"
#include "fvlib/AbstractMatrix.h"
#include "fvlib/CastUtils.h"
#include "fvlib/const.h"
@@ -23,22 +28,22 @@
double el2 = *(double*) b;
if (el1 > el2)
{
- return 1;
+ return 1;
}
if (el1 < el2)
{
- return -1;
+ return -1;
}
if (el1 == el2)
{
- return 0;
+ return 0;
}
// You should never come here...
return -9;
}
-coxph_data::coxph_data( coxph_data &obj)
+coxph_data::coxph_data(const coxph_data &obj)
{
nids = obj.nids;
ncov = obj.ncov;
@@ -52,25 +57,25 @@
order = obj.order;
masked_data = new unsigned short int[nids];
for (int i = 0; i < nids; i++)
- masked_data[i] = 0;
+ masked_data[i] = 0;
}
-coxph_data::coxph_data( phedata &phed, gendata &gend, int snpnum)
+coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum)
{
nids = gend.nids;
masked_data = new unsigned short int[nids];
for (int i = 0; i < nids; i++)
- masked_data[i] = 0;
+ masked_data[i] = 0;
ngpreds = gend.ngpreds;
if (snpnum >= 0)
- ncov = phed.ncov + ngpreds;
+ ncov = phed.ncov + ngpreds;
else
- ncov = phed.ncov;
+ ncov = phed.ncov;
if (phed.noutcomes != 2)
{
- fprintf(stderr,
- "coxph_data: number of outcomes should be 2 (now: %d)\n",
- phed.noutcomes);
- exit(1);
+ fprintf(stderr,
+ "coxph_data: number of outcomes should be 2 (now: %d)\n",
+ phed.noutcomes);
+ exit(1);
}
// X.reinit(nids,(ncov+1));
X.reinit(nids, ncov);
@@ -82,29 +87,29 @@
order.reinit(nids, 1);
for (int i = 0; i < nids; i++)
{
- // 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)
- {
- fprintf(stderr,
- "coxph_data: status not 0/1 (right order: id, fuptime, status ...) %d \n",
- phed.noutcomes);
- exit(1);
- }
+ // 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)
+ {
+ fprintf(stderr,
+ "coxph_data: status not 0/1 (right order: id, fuptime, status ...) %d \n",
+ phed.noutcomes);
+ exit(1);
+ }
}
for (int j = 0; j < phed.ncov; j++)
- for (int i = 0; i < nids; i++)
- X.put((phed.X).get(i, j), i, j);
+ for (int i = 0; i < nids; i++)
+ X.put((phed.X).get(i, j), i, j);
if (snpnum > 0)
- for (int j = 0; j < ngpreds; j++)
- {
- float snpdata[nids];
- gend.get_var(snpnum * ngpreds + j, snpdata);
- for (int i = 0; i < nids; i++)
- X.put(snpdata[i], i, (ncov - ngpreds + j));
- }
+ for (int j = 0; j < ngpreds; j++)
+ {
+ float snpdata[nids];
+ gend.get_var(snpnum * ngpreds + j, snpdata);
+ for (int i = 0; i < nids; i++)
+ X.put(snpdata[i], i, (ncov - ngpreds + j));
+ }
// for (int i=0;i<nids;i++)
// for (int j=0;j<ngpreds;j++)
@@ -112,36 +117,36 @@
for (int i = 0; i < nids; i++)
{
- weights[i] = 1.0;
- offset[i] = 0.0;
- strata[i] = 0;
+ weights[i] = 1.0;
+ offset[i] = 0.0;
+ strata[i] = 0;
}
// sort by time
double tmptime[nids];
int passed_sorted[nids];
for (int i = 0; i < nids; i++)
{
- tmptime[i] = stime[i];
- passed_sorted[i] = 0;
+ tmptime[i] = stime[i];
+ passed_sorted[i] = 0;
}
qsort(tmptime, nids, sizeof(double), cmpfun);
for (int i = 0; i < nids; i++)
{
- int passed = 0;
- for (int j = 0; j < nids; j++)
- if (tmptime[j] == stime[i])
- if (!passed_sorted[j])
- {
- order[i] = j;
- passed_sorted[j] = 1;
- passed = 1;
- break;
- }
- if (passed != 1)
- {
- fprintf(stderr, "cannot recover element %d\n", i);
- exit(1);
- }
+ int passed = 0;
+ for (int j = 0; j < nids; j++)
+ if (tmptime[j] == stime[i])
+ if (!passed_sorted[j])
+ {
+ order[i] = j;
+ passed_sorted[j] = 1;
+ passed = 1;
+ break;
+ }
+ if (passed != 1)
+ {
+ fprintf(stderr, "cannot recover element %d\n", i);
+ exit(1);
+ }
}
stime = reorder(stime, order);
sstat = reorder(sstat, order);
@@ -156,22 +161,22 @@
// stime.print();
// sstat.print();
}
-void coxph_data::update_snp( gendata &gend, int snpnum)
+void coxph_data::update_snp(gendata &gend, int snpnum)
{
// note this sorts by "order"!!!
for (int j = 0; j < ngpreds; j++)
{
- float snpdata[nids];
- for (int i = 0; i < nids; i++)
- masked_data[i] = 0;
- gend.get_var(snpnum * ngpreds + j, snpdata);
- for (int i = 0; i < nids; i++)
- {
- X.put(snpdata[i], (ncov - ngpreds + j), order[i]);
- if (isnan(snpdata[i]))
- masked_data[order[i]] = 1;
- }
+ float snpdata[nids];
+ for (int i = 0; i < nids; i++)
+ masked_data[i] = 0;
+ gend.get_var(snpnum * ngpreds + j, snpdata);
+ for (int i = 0; i < nids; i++)
+ {
+ X.put(snpdata[i], (ncov - ngpreds + j), order[i]);
+ if (isnan(snpdata[i]))
+ masked_data[order[i]] = 1;
+ }
}
// for (int i=0;i<nids;i++)
// for (int j=0;j<ngpreds;j++)
@@ -197,8 +202,8 @@
int nmeasured = 0;
for (int i = 0; i < nids; i++)
- if (masked_data[i] == 0)
- nmeasured++;
+ if (masked_data[i] == 0)
+ nmeasured++;
to.nids = nmeasured;
// std::cout << "nmeasured=" << nmeasured << "\n";
to.ncov = ncov;
@@ -222,25 +227,75 @@
for (int i = 0; i < nids; i++)
{
// std::cout << nids << " " << i << " " << masked_data[i] << "\n";
- if (masked_data[i] == 0)
- {
- (to.weights).put(weights.get(i, 1), j, 1);
- (to.stime).put(stime.get(i, 1), j, 1);
- (to.sstat).put(sstat.get(i, 1), j, 1);
- (to.offset).put(offset.get(i, 1), j, 1);
- (to.strata).put(strata.get(i, 1), j, 1);
- (to.order).put(order.get(i, 1), j, 1);
- for (int nc = 0; nc < dim1X; nc++)
- (to.X).put(X.get(nc, i), nc, j);
- j++;
- }
+ if (masked_data[i] == 0)
+ {
+ (to.weights).put(weights.get(i, 1), j, 1);
+ (to.stime).put(stime.get(i, 1), j, 1);
+ (to.sstat).put(sstat.get(i, 1), j, 1);
+ (to.offset).put(offset.get(i, 1), j, 1);
+ (to.strata).put(strata.get(i, 1), j, 1);
+ (to.order).put(order.get(i, 1), j, 1);
+ for (int nc = 0; nc < dim1X; nc++)
+ (to.X).put(X.get(nc, i), nc, j);
+ j++;
+ }
}
// std::cout << " !!! just after cycle !!! ";
//delete [] to.masked_data;
to.masked_data = new unsigned short int[to.nids];
for (int i = 0; i < to.nids; i++)
- to.masked_data[i] = 0;
+ to.masked_data[i] = 0;
//fprintf(stdout,"get_unmasked: %i %i %i\n",to.nids,dim2X,dim2Y);
return (to);
}
+
+coxph_reg::coxph_reg(coxph_data &cdatain)
+{
+ coxph_data cdata = cdatain.get_unmasked_data();
+ beta.reinit(cdata.X.nrow, 1);
+ sebeta.reinit(cdata.X.nrow, 1);
+ loglik = -9.999e+32;
+ sigma2 = -1.;
+ chi2_score = -1.;
+ niter = 0;
+}
+
+void coxph_reg::estimate(coxph_data &cdatain, int verbose, int maxiter,
+ double eps, double tol_chol, int model, int interaction, int ngpreds,
+ bool iscox, int nullmodel)
+{
+ coxph_data cdata = cdatain.get_unmasked_data();
+ mematrix<double> X = t_apply_model(cdata.X, model, interaction, ngpreds,
+ iscox, nullmodel);
+ // X.print();
+ int length_beta = X.nrow;
+ beta.reinit(length_beta, 1);
+ sebeta.reinit(length_beta, 1);
+ mematrix<double> newoffset = cdata.offset;
+ newoffset = cdata.offset - (cdata.offset).column_mean(0);
+ mematrix<double> means(X.nrow, 1);
+ for (int i = 0; i < X.nrow; i++)
+ beta[i] = 0.;
+ mematrix<double> u(X.nrow, 1);
+ mematrix<double> imat(X.nrow, X.nrow);
+ double work[X.ncol * 2 + 2 * (X.nrow) * (X.nrow) + 3 * (X.nrow)];
+ double loglik_int[2];
+ int flag;
+ double sctest = 1.0;
+
+ //TODO(maarten): this function works only in combination with eigen remove .data() from arguments to make is available under old matrix
+ 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);
+ for (int i = 0; i < X.nrow; i++)
+ {
+ sebeta[i] = sqrt(imat.get(i, i));
+ }
+ loglik = loglik_int[1];
+ niter = maxiter;
+}
+
Modified: branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_data.h 2012-11-12 21:54:03 UTC (rev 1018)
@@ -10,18 +10,32 @@
#if EIGEN
#include "eigen_mematrix.h"
-#include "eigen_mematri1.h"
+#include "eigen_mematrix.cpp"
#else
#include "mematrix.h"
#include "mematri1.h"
#endif
+#include "reg1.h"
#include "gendata.h"
#include "phedata.h"
-class coxph_data
-{
+class coxph_data {
public:
+
+ coxph_data get_unmasked_data();
+ coxph_data()
+ {
+ nids = 0;
+ ncov = 0;
+ ngpreds = 0;
+ masked_data = NULL;
+ }
+ coxph_data(const coxph_data &obj);
+ coxph_data(phedata &phed, gendata &gend, int snpnum);
+ void update_snp(gendata &gend, int snpnum);
+ ~coxph_data();
+
int nids;
int ncov;
int ngpreds;
@@ -33,15 +47,21 @@
mematrix<double> X;
mematrix<int> order;
unsigned short int * masked_data;
- coxph_data get_unmasked_data();
- coxph_data()
- {
- }
- coxph_data( coxph_data &obj);
- coxph_data( phedata &phed, gendata &gend, int snpnum);
- void update_snp( gendata &gend, int snpnum);
- ~coxph_data();
+};
+class coxph_reg {
+public:
+ mematrix<double> beta;
+ mematrix<double> sebeta;
+ mematrix<double> residuals;
+ double sigma2;
+ double loglik;
+ double chi2_score;
+ int niter;
+ coxph_reg(coxph_data &cdatain);
+ void estimate(coxph_data &cdatain, int verbose, int maxiter, double eps,
+ double tol_chol, int model, int interaction, int ngpreds,
+ bool iscox, int nullmodel = 0);
};
#endif /* COXPH_DATA_H_ */
Deleted: branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.cpp 2012-11-12 21:54:03 UTC (rev 1018)
@@ -1,55 +0,0 @@
-
-#include "coxph_reg.h"
-#include "coxph_data.h"
-extern "C"
-{
-#include "survproto.h"
-}
-
-coxph_reg::coxph_reg(const coxph_data &cdatain) {
- coxph_data cdata = cdatain.get_unmasked_data();
- beta.reinit(cdata.X.nrow, 1);
- sebeta.reinit(cdata.X.nrow, 1);
- loglik = -9.999e+32;
- sigma2 = -1.;
- chi2_score = -1.;
- niter = 0;
-}
-coxph_reg::~coxph_reg() {
- // delete beta;
- // delete sebeta;
-}
-void coxph_reg::estimate( coxph_data &cdatain, int verbose, int maxiter,
- double eps, double tol_chol, int model, int interaction, int ngpreds,
- bool iscox, int nullmodel = 0) {
- // cout << "model = " << model << "\n";
- // cdata.X.print();
- coxph_data cdata = cdatain.get_unmasked_data();
- mematrix<double> X = t_apply_model(cdata.X, model, interaction, ngpreds,
- iscox, nullmodel);
- // X.print();
- int length_beta = X.nrow;
- beta.reinit(length_beta, 1);
- sebeta.reinit(length_beta, 1);
- mematrix<double> newoffset = cdata.offset;
- newoffset = cdata.offset - (cdata.offset).column_mean(0);
- mematrix<double> means(X.nrow, 1);
- for (int i = 0; i < X.nrow; i++)
- beta[i] = 0.;
- mematrix<double> u(X.nrow, 1);
- mematrix<double> imat(X.nrow, X.nrow);
- double work[X.ncol * 2 + 2 * (X.nrow) * (X.nrow) + 3 * (X.nrow)];
- double loglik_int[2];
- int flag;
- double sctest = 1.0;
-//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];
- niter = maxiter;
-}
-
Deleted: branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/coxph_reg.h 2012-11-12 21:54:03 UTC (rev 1018)
@@ -1,35 +0,0 @@
-/*
- * coxph_reg.h
- *
- * Created on: Oct 6, 2012
- * Author: mkooyman
- */
-
-#ifndef COXPH_REG_H_
-#define COXPH_REG_H_
-#include "coxph_data.h"
-
-#if EIGEN
-#include "eigen_mematrix.h"
-#else
-#include "mematrix.h"
-#endif
-
-class coxph_reg {
-public:
- mematrix<double> beta;
- mematrix<double> sebeta;
- mematrix<double> residuals;
- double sigma2;
- double loglik;
- double chi2_score;
- int niter;
-
- coxph_reg(const coxph_data &cdatain);
- virtual ~coxph_reg();
- void coxph_reg::estimate( coxph_data &cdatain, int verbose, int maxiter,
- double eps, double tol_chol, int model, int interaction, int ngpreds,
- bool iscox, int nullmodel = 0);
-};
-
-#endif /* COXPH_REG_H_ */
Modified: branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-11-12 21:54:03 UTC (rev 1018)
@@ -17,7 +17,7 @@
#if EIGEN
#include "eigen_mematrix.h"
-#include "eigen_mematri1.h"
+#include "eigen_mematrix.cpp"
#else
#include "mematrix.h"
#include "mematri1.h"
@@ -80,8 +80,7 @@
nlin++;
}
nlin--; // Subtract one, the previous loop added 1 too much
- }
- else
+ } else
{
std::cerr << "mlinfo: cannot open file " << filename << endl;
exit(1);
@@ -217,8 +216,7 @@
row++;
}
myfile.close();
- }
- else
+ } else
{
fprintf(stderr, "error: inv file: cannot open file '%s'\n", filename_);
}
Deleted: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h 2012-11-08 21:36:15 UTC (rev 1017)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h 2012-11-12 21:54:03 UTC (rev 1018)
@@ -1,349 +0,0 @@
-#ifndef EIGEN_MEMATRI1_H
-#define EIGEN_MEMATRI1_H
-#include "eigen_mematrix.h"
-#include <Eigen/Dense>
-#include <Eigen/LU>
-#include <string>
-#include <cstdarg>
-#include <cstdio>
-#include <cstdlib>
-
-
-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);
- }
- this->nrow = nr;
- this->ncol = nc;
- this->nelements = nr * nc;
- this->data.resize(nr,nc);
-}
-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;
- }
- 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+(const 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 = data.cols();
- temp.nrow = data.rows();
- temp.nelements = temp.nrow * temp.ncol;
-
- return temp;
-}
-
-template<class DT>
-mematrix<DT> mematrix<DT>::operator-(const DT toadd)
-{
- mematrix<DT> temp(nrow, ncol);
- temp.data = data.array() - toadd;
- return temp;
-}
-template<class DT>
-mematrix<DT> mematrix<DT>::operator-(const 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*(const 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;
-}
-
-template<class DT>
-mematrix<DT> mematrix<DT>::operator*(const 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)
-{
-// printf("put val:%f nr=%i nc=%i \n", value, nr, nc);
-// printf("mat nr=%i nc=%i \n", data.rows(),data.cols());
-// printf("mat nr=%i nc=%i \n", nrow,ncol);
- 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(const 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(const mematrix<DT> &M)
-{
-//cout << "[DEBUG TRANSPOSE PRE]nrow=" << M.nrow << "; ncol=" << M.ncol << "; nelements=" << M.nelements;
-
- mematrix<DT> temp;
- temp.data = M.data.transpose();
- temp.ncol = M.nrow;
- temp.nrow = M.ncol;
- temp.nelements = M.nelements;
-//cout << "[DEBUG TRANSPOSE post]nrow=" << temp.nrow << "; ncol=" << temp.ncol << "; nelements=" << temp.nelements;
-
- return temp;
-}
-
-template<class DT>
-mematrix<DT> reorder(const mematrix<DT> &M, const 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);
-//FIXME(maarten): commented out to get compilation running
-// 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(const 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( const mematrix<DT> &M, const 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;
- //make a array of the first row of D in the same way orientation as M.data.row(i).array()
- Array<DT,Dynamic,Dynamic> row=D.data.block(0,0,M.ncol,1).transpose();
-
- for (int i = 0; i < temp.nrow; i++)
- {
- temp.data.row(i) = M.data.row(i).array() * row;
- }
- return temp;
-}
-
-#endif
Copied: branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.cpp (from rev 1008, branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematri1.h)
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.cpp (rev 0)
+++ branches/ProbABEL-refactoring/ProbABEL/src/eigen_mematrix.cpp 2012-11-12 21:54:03 UTC (rev 1018)
@@ -0,0 +1,355 @@
+#ifndef EIGEN_MEMATRI1_H
+#define EIGEN_MEMATRI1_H
+#include "eigen_mematrix.h"
+#include <Eigen/Dense>
+#include <Eigen/LU>
+#include <string>
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+
+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);
+ }
+ this->nrow = nr;
+ this->ncol = nc;
+ this->nelements = nr * nc;
+ this->data.resize(nr, nc);
+}
+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;
+ }
+ 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+(const 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 = data.cols();
+ temp.nrow = data.rows();
+ temp.nelements = temp.nrow * temp.ncol;
+
+ return temp;
+}
+
+template<class DT>
+mematrix<DT> mematrix<DT>::operator-(const DT toadd)
+{
+ mematrix<DT> temp(nrow, ncol);
+ temp.data = data.array() - toadd;
+ return temp;
+}
+template<class DT>
+mematrix<DT> mematrix<DT>::operator-(const 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*(const 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;
+// fprintf(stderr, "mematrix*: (%d,%d) and (%d,%d):result%d\n", nrow,
+// ncol, M.nrow, M.ncol,temp.nrow * temp.ncol);
+// std::cout.flush();
+
+ return temp;
+}
+
+template<class DT>
+mematrix<DT> mematrix<DT>::operator*(const 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;
+// fprintf(stderr, "mematrix*: (%d,%d) and (%d,%d):result%d\n", nrow,
+// ncol, M->nrow, M->ncol,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 !NDEBUG
+ 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);
+ }
+#endif
+ DT temp = data(nr, nc);
+ return temp;
+}
+template<class DT>
+void mematrix<DT>::put(DT value, int nr, int nc)
+{
+#if !NDEBUG
+ 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);
+ }
+#endif
+ 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(const 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>
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1018
More information about the Genabel-commits
mailing list