From 9d2966e8f3c58f375c04a80b99bcd5d31c80c736 Mon Sep 17 00:00:00 2001
From: Duncan Mortimer <duncan.mortimer@ndcn.ox.ac.uk>
Date: Thu, 19 Jul 2007 08:36:24 +0000
Subject: [PATCH] Reverted to 1.5 code and changed exit code for unsuccessful
running to 1
---
ccops.cc | 168 ++++---------------------------------------------------
1 file changed, 10 insertions(+), 158 deletions(-)
diff --git a/ccops.cc b/ccops.cc
index 349f89f..810f9f0 100644
--- a/ccops.cc
+++ b/ccops.cc
@@ -9,13 +9,11 @@
#include "ccopsOptions.h"
#include <vector>
#include <algorithm>
-#include "dpm_gibbs.h"
using namespace std;
using namespace NEWIMAGE;
using namespace NEWMAT;
using namespace CCOPS;
-
void spect_reord(SymmetricMatrix& A,ColumnVector& r,ColumnVector& y){
SymmetricMatrix Q=-A;
DiagonalMatrix t(Q.Nrows());
@@ -72,117 +70,6 @@ using namespace CCOPS;
}
-// calculate pca decomposition using the covariance method
-// X is the data dxn matrix with n observations (data points) and d variables (dimensions)
-ReturnMatrix pcacov(Matrix& X,const float& perc){
- int n=X.Nrows();
- int d=X.Ncols();
- // de-mean
- cout<<"de-mean"<<endl;
- ColumnVector Xmean(d);
- for(int j=1;j<=d;j++)
- Xmean(j) = X.Column(j).Sum()/n;
- for(int j=1;j<=d;j++){
- for(int i=1;i<=n;i++)
- X(i,j) -= Xmean(j);
- }
-
- // calculate covariance
- cout<<"covariance"<<endl;
- SymmetricMatrix C(d);
-
- if(d<n)
- C << (X.t()*X)/n;
- else
- C << X*X.t()/n;
-
- // eigenvalues
- cout<<"eigenvalues"<<endl;
- Matrix V;
- DiagonalMatrix D;
- EigenValues(C,D,V);
-
- // select subset
- cout<<"subset"<<endl;
- float cumsum=0,total=D.Trace();
- int dim=0;
- for(int i=D.Nrows();i>=1;i--){
- cumsum += D(i);
- if(cumsum/total < perc){dim++;}
- else{break;}
- }
- if(dim<=2)dim=2;
-
- dim=5;
-
- Matrix v(V.Nrows(),dim);
- ColumnVector lam(dim);
- for(int j=1;j<=dim;j++){
- v.Column(j) = V.Column(V.Ncols()-j+1);
- lam(j)=D(V.Ncols()-j+1);
- }
-
- //cout<<"zscores"<<endl;
- // convert to z-scores
- //for(int i=1;i<=d;i++)
- //X.Row(i) /= sqrt(C(i,i));
- // reconstruct data
- cout<<"data"<<endl;
- Matrix data(dim,n);
-
- if(!(d<n)){
- v = (X.t() * v);
- for(int j=1;j<=dim;j++)
- v.Column(j) /= sqrt(n*lam(j));
- }
- data = X*v;
-
-
- data.Release();
- return data;
-}
-
-void dpm_reord(Matrix& A,ColumnVector& r,ColumnVector& y){
- cout << "start dpm reordering" <<endl;
- // pca preprocess
- cout << "pca preprocessing" <<endl;
- float perc=.90;
- Matrix data;
- OUT(A.Nrows());
- OUT(A.Ncols());
- write_ascii_matrix(A,"nonpreprocessed_data");
- OUT(A.Nrows());
- OUT(A.Ncols());
- data = pcacov(A,perc);
- write_ascii_matrix(data,"preprocessed_data");
-
- // dpm
- cout << "dpm clustering" <<endl;
- int numiter=2000;
- int burnin=1000;
- int sampleevery=1;
- DPM_GibbsSampler gs(data,numiter,burnin,sampleevery);
- gs.init();
- gs.run();
- // save
- int n = data.Nrows();
- vector< pair<float,int> > myvec;
- ColumnVector z(n);
- z = gs.get_dataindex();
- //z = gs.get_mldataindex();
- for(int i=1;i<=n;i++){
- pair<float,int> mypair;
- mypair.first = z(i);
- mypair.second = i;
- myvec.push_back(mypair);
- }
- sort(myvec.begin(),myvec.end());
- r.ReSize(n);y.ReSize(n);
- for(int i=1;i<=n;i++){
- y(i)=myvec[i-1].first;
- r(i)=myvec[i-1].second;
- }
-}
void rem_zrowcol(const Matrix& myOMmat,const Matrix& coordmat,const Matrix& tractcoordmat,const bool coordbool,const bool tractcoordbool,Matrix& newOMmat,Matrix& newcoordmat, Matrix& newtractcoordmat)
{
@@ -413,12 +300,10 @@ void add_connexity(SymmetricMatrix& CtCt,const Matrix& coord,const float p=.5){
int main ( int argc, char **argv ){
ccopsOptions& opts = ccopsOptions::getInstance();
int success=opts.parse_command_line(argc,argv);
- if(!success) return 0;
+ if(!success) return 1;
string ip=opts.inmatrix.value();
make_basename(ip);
- srand(time(NULL));
-
ColumnVector y1,r1,y2,r2;
volume<int> myOM;
volume<int> coordvol;
@@ -525,9 +410,9 @@ int main ( int argc, char **argv ){
rem_zrowcol(myOMmat,mycoordmat,mytractcoordmat,coordbool,tractcoordbool,newOMmat,newcoordmat,newtractcoordmat);
- //cerr<<"NOW"<<endl;
- //cerr<<myOMmat.MaximumAbsoluteValue()<<endl;
- //cerr<<newOMmat.MaximumAbsoluteValue()<<endl;
+ // cerr<<"NOW"<<endl;
+ // cerr<<myOMmat.MaximumAbsoluteValue()<<endl;
+ // cerr<<newOMmat.MaximumAbsoluteValue()<<endl;
//write_ascii_matrix("ncm",newcoordmat);
// write_ascii_matrix("nctm",newtractcoordmat);
@@ -542,12 +427,12 @@ int main ( int argc, char **argv ){
cerr<<"Computing correlation"<<endl;
- SymmetricMatrix CtCt(newOMmat.Nrows());
- //CtCt << corrcoef(newOMmat.t());
- //CtCt << CtCt+1;
+ SymmetricMatrix CtCt;
+ CtCt << corrcoef(newOMmat.t());
+ CtCt << CtCt+1;
// adding connexity constraint
- if(opts.connexity.value()!=0 && !coordbool){
+ if(!coordbool){
cerr<<"WARNING !! No coordinates provided. I cannot apply any connexity constraint."<<endl;
}
else{
@@ -582,29 +467,7 @@ int main ( int argc, char **argv ){
}
else{
cerr<<"Starting First Reordering"<<endl;
- if(opts.scheme.value()=="dpm"){
- OUT(myOMmat.Ncols());
- OUT(myOMmat.Nrows());
- OUT(newOMmat.Ncols());
- OUT(newOMmat.Nrows());
- dpm_reord(newOMmat,r1,y1);
- //Matrix A;
- //A=CtCt;
- //dpm_reord(A,r1,y1);
- }
- else if(opts.scheme.value()=="kmeans"){
- int kk=opts.kmeans.value();
- if(kk<=1){
- cerr << "Error using kmeans. Number of clusters must be >=2" << endl;
- return -1;
- }
- cout << "... using kmeans" << endl;
- //kmeans_reord(CtCt,r1,y1,kk);
- }
- else
- spect_reord(CtCt,r1,y1);
-
-
+ spect_reord(CtCt,r1,y1);
cerr<<"Permuting seed CC matrix"<<endl;
@@ -628,18 +491,7 @@ int main ( int argc, char **argv ){
write_ascii_matrix(y1,base+"y1");
save_volume(outCCvol,"reord_CC_"+base);
save_volume(outcoords,"coords_for_reord_"+base);
- if(opts.scheme.value()=="dpm" || opts.scheme.value()=="kmeans"){
- volume<int> mask;
- read_volume(mask,opts.mask.value());
- mask = 0;
- for(int i=0;i<outcoords.xsize();i++){
- mask(outcoords(i,0,0),
- outcoords(i,1,0),
- outcoords(i,2,0)) = (int)y1(i+1) + 1;
- }
- save_volume(mask,"reord_mask_"+base);
- }
-
+
}
if(opts.reord2.value()){
--
GitLab