[Genabel-commits] r871 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 29 02:07:08 CEST 2012
Author: maartenk
Date: 2012-03-29 02:07:07 +0200 (Thu, 29 Mar 2012)
New Revision: 871
Added:
branches/ProbABEL-refactoring/ProbABEL/src/utilities.cpp
branches/ProbABEL-refactoring/ProbABEL/src/utilities.h
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
branches/ProbABEL-refactoring/ProbABEL/src/mematrix.h
branches/ProbABEL-refactoring/ProbABEL/src/phedata.cpp
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
Log:
merged with upstream commits(made last 2 weeks
)
Modified: branches/ProbABEL-refactoring/ProbABEL/src/data.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-03-28 22:46:52 UTC (rev 870)
+++ branches/ProbABEL-refactoring/ProbABEL/src/data.cpp 2012-03-29 00:07:07 UTC (rev 871)
@@ -2,8 +2,6 @@
#include <sstream>
#include <fstream>
-
-
#include "fvlib/AbstractMatrix.h"
#include "fvlib/CastUtils.h"
#include "fvlib/const.h"
@@ -24,37 +22,38 @@
extern bool is_interaction_excluded;
-
unsigned int Nmeasured(char * fname, int nphenocols, int npeople) {
int ncov = nphenocols - 2;
int nids_all = npeople;
- FILE * infile;
// first pass -- find unmeasured people
- if ((infile = fopen(fname, "r")) == NULL) {
- fprintf(stderr, "Nmeasured: cannot open file %s\n", fname);
+ std::ifstream infile(fname);
+ if (!infile) {
+ std::cerr << "Nmeasured: cannot open file " << fname << endl;
}
char tmp[100];
for (int i = 0; i < nphenocols; i++) {
- fscanf(infile, "%s", tmp);
- // printf("%s ",tmp);
- } //printf("\n");
+ infile >> tmp;
+ }
unsigned short int * allmeasured = new unsigned short int[npeople];
int nids = 0;
for (int i = 0; i < npeople; i++) {
allmeasured[i] = 1;
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
for (int j = 1; j < nphenocols; j++) {
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
if (tmp[0] == 'N' || tmp[0] == 'n')
allmeasured[i] = 0;
}
if (allmeasured[i] == 1)
nids++;
}
- fclose(infile);
+ infile.close();
+
+ delete[] allmeasured;
+
return (nids);
}
@@ -171,14 +170,12 @@
int cmpfun(const void *a, const void *b) {
double el1 = *(double*) a;
double el2 = *(double*) b;
- int returnvalue;
if (el1 > el2)
- returnvalue = 1;
+ return 1;
if (el1 < el2)
- returnvalue = -1;
+ return -1;
if (el1 == el2)
- returnvalue = 0;
- return returnvalue;
+ return 0;
}
coxph_data::coxph_data(const coxph_data &obj) {
@@ -225,7 +222,9 @@
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 ...)\n", phed.noutcomes);
+ fprintf(stderr,
+ "coxph_data: status not 0/1 (right order: id, fuptime, status ...)\n",
+ phed.noutcomes);
exit(1);
}
}
@@ -369,23 +368,28 @@
}
mlinfo::mlinfo(char * filename, char * mapname) {
- FILE * infile = fopen(filename, "r");
- if (infile == NULL) {
- fprintf(stderr, "mlinfo: cannot open file %s", filename);
- exit(1);
- }
+
char tmp[100];
unsigned int nlin = 0;
- while (fscanf(infile, "%s", tmp) != EOF) {
- nlin++;
+ std::ifstream infile(filename);
+ if (infile.is_open()) {
+ while (infile.good()) {
+ infile >> tmp;
+ nlin++;
+ }
+ nlin--; // Subtract one, the previous loop added 1 too much
+ } else {
+ std::cerr << "mlinfo: cannot open file " << filename << endl;
+ exit(1);
}
- fclose(infile);
+ infile.close();
+
if (nlin % 7) {
- fprintf(stderr, "mlinfo: number of columns != 7 in %s", filename);
+ std::cerr << "mlinfo: number of columns != 7 in " << filename << endl;
exit(1);
}
nsnps = int(nlin / 7) - 1;
- printf("Number of SNPs = %d\n", nsnps);
+ std::cout << "Number of SNPs = " << nsnps << endl;
name = new std::string[nsnps];
A1 = new std::string[nsnps];
A2 = new std::string[nsnps];
@@ -394,36 +398,41 @@
Quality = new double[nsnps];
Rsq = new double[nsnps];
map = new std::string[nsnps];
- if ((infile = fopen(filename, "r")) == NULL) {
- fprintf(stderr, "mlinfo: cannot open file %s", filename);
+
+ infile.open(filename);
+ if (!infile) { // file couldn't be opened
+ std::cerr << "mlinfo: cannot open file " << filename << endl;
exit(1);
}
+ /* Read the header and discard it */
for (int i = 0; i < 7; i++)
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
+
for (int i = 0; i < nsnps; i++) {
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
name[i] = tmp;
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
A1[i] = tmp;
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
A2[i] = tmp;
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
Freq1[i] = atof(tmp);
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
MAF[i] = atof(tmp);
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
Quality[i] = atof(tmp);
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
Rsq[i] = atof(tmp);
map[i] = "-999";
}
- fclose(infile);
+ infile.close();
+
if (mapname != NULL) {
std::ifstream instr(mapname);
int BFS = 1000;
char line[BFS], tmp[BFS];
if (!instr.is_open()) {
- fprintf(stderr, "mlinfo: cannot open file %s", mapname);
+ std::cerr << "mlinfo: cannot open file " << mapname << endl;
exit(1);
}
instr.getline(line, BFS);
@@ -445,65 +454,66 @@
delete[] mlinfo::Rsq;
delete[] mlinfo::map;
}
+};
//_________________________________________Maksim_start
InvSigma::InvSigma(const char * filename_, phedata * phe) {
- filename = filename_;
- npeople = phe->nids;
- std::ifstream myfile(filename_);
- char * line = new char[MAXIMUM_PEOPLE_AMOUNT];
- double val;
- std::string id;
- unsigned row = 0, col = 0;
+filename = filename_;
+npeople = phe->nids;
+std::ifstream myfile(filename_);
+char * line = new char[MAXIMUM_PEOPLE_AMOUNT];
+double val;
+std::string id;
+unsigned row = 0, col = 0;
- matrix.reinit(npeople, npeople);
+matrix.reinit(npeople, npeople);
- //idnames[k], if (allmeasured[i]==1)
+//idnames[k], if (allmeasured[i]==1)
- if (myfile.is_open()) {
- while (myfile.getline(line, MAXIMUM_PEOPLE_AMOUNT)) {
+if (myfile.is_open()) {
+ while (myfile.getline(line, MAXIMUM_PEOPLE_AMOUNT)) {
- std::stringstream line_stream(line);
- line_stream >> id;
+ std::stringstream line_stream(line);
+ line_stream >> id;
- if (phe->idnames[row] != id) {
- std::cerr << "error:in row " << row << " id="
- << phe->idnames[row]
- << " in inverse variance matrix but id=" << id
- << " must be there. Wrong inverse variance matrix (only measured id must be there)\n";
- exit(1);
- }
+ if (phe->idnames[row] != id) {
+ std::cerr << "error:in row " << row << " id=" << phe->idnames[row]
+ << " in inverse variance matrix but id=" << id
+ << " must be there. Wrong inverse variance matrix (only measured id must be there)\n";
+ exit(1);
+ }
- while (line_stream >> val) {
- matrix.put(val, row, col);
- col++;
- }
+ while (line_stream >> val) {
+ matrix.put(val, row, col);
+ col++;
+ }
- if (col != npeople) {
- fprintf(stderr,
- "error: inv file: Number of columns in row %d equals to %d but people amount is %d\n",
- row, col, npeople);
- myfile.close();
- exit(1);
- }
- col = 0;
- row++;
+ if (col != npeople) {
+ fprintf(stderr,
+ "error: inv file: Number of columns in row %d equals to %d but people amount is %d\n",
+ row, col, npeople);
+ myfile.close();
+ exit(1);
}
- myfile.close();
- } else {
- fprintf(stderr, "error: inv file: cannot open file '%s'\n", filename_);
+ col = 0;
+ row++;
}
+ myfile.close();
+} else {
+ fprintf(stderr, "error: inv file: cannot open file '%s'\n", filename_);
+}
+delete[] line;
}
;
InvSigma::~InvSigma() {
- //af af
+//af af
}
mematrix<double> & InvSigma::get_matrix(void) {
- return matrix;
+return matrix;
}
//________________________________________Maksim_end
Modified: branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-03-28 22:46:52 UTC (rev 870)
+++ branches/ProbABEL-refactoring/ProbABEL/src/gendata.cpp 2012-03-29 00:07:07 UTC (rev 871)
@@ -11,8 +11,6 @@
#include "mematri1.h"
#include "utilities.h"
-
-
void gendata::get_var(int var, float * data) {
if (DAG == NULL)
for (int i = 0; i < G.nrow; i++)
@@ -84,10 +82,11 @@
G.reinit(nids, (nsnps * ngpreds));
- FILE * infile;
+ std::ifstream infile;
- if ((infile = fopen(fname, "r")) == NULL) {
- fprintf(stderr, "gendata: cannot open file %s\n", fname);
+ infile.open(fname);
+ if (!infile) {
+ std::cerr << "gendata: cannot open file " << fname << endl;
}
char tmp[100], tmpn[100];
@@ -99,7 +98,7 @@
if (skipd > 0) {
// int ttt;
char ttt[100];
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
// sscanf(tmp,"%d->%s",&ttt, tmpn);
// these changes are thanks to BMM & BP :)
// sscanf(tmp,"%s->%s",&ttt, tmpn);
@@ -117,19 +116,20 @@
"phenofile and dosefile did not match at line %d ",
i + 2);
cerr << "(" << tmpid << " != " << idnames[k] << ")\n";
- fclose(infile);
+ infile.close();
exit(1);
}
}
for (int j = 1; j < skipd; j++) {
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
}
for (int j = 0; j < (nsnps * ngpreds); j++) {
- int a = fscanf(infile, "%s", tmp);
- if (!a || a == EOF) {
- fprintf(stderr,
- "cannot read dose-file: check skipd and ngpreds parameters\n");
- fclose(infile);
+ if (infile.good()) {
+ infile >> tmp;
+ } else {
+ std::cerr
+ << "cannot read dose-file: check skipd and ngpreds parameters\n";
+ infile.close();
exit(1);
}
G.put(atof(tmp), k, j);
@@ -137,11 +137,11 @@
k++;
} else {
for (int j = 0; j < skipd; j++)
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
for (int j = 0; j < (nsnps * ngpreds); j++)
- fscanf(infile, "%s", tmp);
+ infile >> tmp;
}
- fclose(infile);
+ infile.close();
}
// HERE NEED A NEW CONSTRUCTOR BASED ON DATABELBASECPP OBJECT
gendata::~gendata() {
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-03-28 22:46:52 UTC (rev 870)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-03-29 00:07:07 UTC (rev 871)
@@ -11,7 +11,7 @@
// Author: Yurii S. Aulchenko (cox, log, lin regressions)
// Modified by: L.C. Karssen,
// Maksim V. Struchalin
-//
+//
// modified 14-May-2009 by MVS: interaction with SNP, interaction with SNP with exclusion of interacted covariates,
// mmscore implemented (poor me)
// modified 20-Jul-2009 by YSA: small changes, bug fix in mmscore option
@@ -264,7 +264,9 @@
mlinfo mli(mlinfofilename, mapfilename);
int nsnps = mli.nsnps;
- phedata phd(phefilename, noutcomes, npeople, interaction, iscox);
+ phedata phd;
+ phd.set_is_interaction_excluded(is_interaction_excluded);
+ phd.setphedata(phefilename, noutcomes, npeople, interaction, iscox);
int interaction_cox = interaction;
#if COXPH
@@ -349,21 +351,24 @@
// estimate null model
double null_loglik = 0.;
#if COXPH
- coxph_data nrgd(phd,gtd,-1);
-#else
- regdata nrgd(phd, gtd, -1);
+ coxph_data nrgd=coxph_data(phd,gtd,-1);
+#else
+ regdata nrgd = regdata(phd, gtd, -1);
#endif
std::cout << " loaded null data ...";
#if LOGISTIC
- logistic_reg nrd(nrgd);
+ logistic_reg nrd=logistic_reg(nrgd);
nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, interaction, ngpreds, invvarmatrix, robust, 1);
#elif LINEAR
- linear_reg nrd(nrgd);
+
+ linear_reg nrd=linear_reg(nrgd);
+
nrd.estimate(nrgd,0,CHOLTOL,0, interaction, ngpreds, invvarmatrix, robust, 1);
#elif COXPH
coxph_reg nrd(nrgd);
+
nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, interaction, ngpreds, 1);
#endif
null_loglik = nrd.loglik;
@@ -672,10 +677,8 @@
}
*/
// exit(1);
-
//________________________________________________________________
//Maksim, 9 Jan, 2009
-
int maxmod = 5;
int start_pos, end_pos;
Modified: branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h 2012-03-28 22:46:52 UTC (rev 870)
+++ branches/ProbABEL-refactoring/ProbABEL/src/mematri1.h 2012-03-29 00:07:07 UTC (rev 871)
@@ -10,142 +10,141 @@
// 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);
- }
- // fprintf(stderr,"mematrix(nr,nc): can allocate memory (%d,%d)\n",nrow,ncol);
+ 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);
+ }
+ // fprintf(stderr,"mematrix(nr,nc): can allocate memory (%d,%d)\n",nrow,ncol);
}
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);
- }
- // fprintf(stderr,"mematrix const(mematrix): can allocate memory (%d,%d)\n",M.nrow,M.ncol);
- for (int i = 0; i < M.ncol * M.nrow; i++)
- data[i] = M.data[i];
+ 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);
+ }
+ // fprintf(stderr,"mematrix const(mematrix): can allocate memory (%d,%d)\n",M.nrow,M.ncol);
+ 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];
- // fprintf(stderr,"mematrix=: can allocate memory (%d,%d)\n",M.nrow,M.ncol);
- }
- return *this;
+ 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];
+ // fprintf(stderr,"mematrix=: can allocate memory (%d,%d)\n",M.nrow,M.ncol);
+ }
+ 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];
+ 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;
+ 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;
+ 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;
+ 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;
+ 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;
+ // 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;
+ 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;
}
//
@@ -153,132 +152,132 @@
//
template<class DT>
void mematrix<DT>::reinit(int nr, int nc) {
- if (nelements > 0)
- delete[] data;
- 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);
- }
+ if (nelements > 0)
+ delete[] data;
+ 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>
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 * ncol + nc];
- return temp;
+ 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 * ncol + 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 * ncol + nc] = value;
+ 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 * 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;
+ 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++)
- cout << data[i * ncol + j] << "\t";
- cout << "\n";
- }
+ 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[i * ncol + j] << "\t";
+ cout << "\n";
+ }
}
template<class DT>
void mematrix<DT>::delete_column(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];
- }
+ 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(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];
- }
+ 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];
+ }
}
@@ -287,210 +286,209 @@
//
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;
+ 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;
+ 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;
+ 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) {
- 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;
+ 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<DT> productMatrDiag(mematrix<DT> &M, mematrix<DT> &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;
+ 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;
+ 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);
- }
+ 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;
+ 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);
- }
- }
+ 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;
+ 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];
- }
- /*
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 871
More information about the Genabel-commits
mailing list