From 70978cde28e60291eeeff3259c07ab5a48c9fed9 Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 30 Jul 2007 16:42:49 +0000
Subject: [PATCH] MELODIC 3.02 with local FDR
---
fsl_regfilt.cc | 12 ++++----
meldata.cc | 83 ++++++++++++++++++++++++++++++++++++--------------
meldata.h | 7 ++---
melgmix.cc | 60 +++++++++++++++++++++++-------------
melhlprfns.cc | 16 +++++++---
melhlprfns.h | 1 +
melodic.cc | 2 ++
melodic.h | 19 ++++++++++--
meloptions.cc | 7 ++---
meloptions.h | 12 ++------
melreport.cc | 30 +++++++++++++-----
melreport.h | 1 +
test.cc | 8 +----
13 files changed, 168 insertions(+), 90 deletions(-)
diff --git a/fsl_regfilt.cc b/fsl_regfilt.cc
index e8f4c1c..bb16d13 100644
--- a/fsl_regfilt.cc
+++ b/fsl_regfilt.cc
@@ -23,17 +23,17 @@ using namespace std;
// printed out as the help or usage message
string title=string("fsl_regfilt (Version 1.0)")+
- string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
- string(" Data filtering by regressing out part of a design matrix\n \n")+
- string(" using simple OLS regression on 4D images\n\n");
- string examples="fsl_regfilt -i <input> -d <design> -f -o <out> [options]";
+ string("\n\n Copyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
+ string(" Data de-noising by regressing out part of a design matrix\n")+
+ string(" using simple OLS regression on 4D images");
+ string examples="fsl_regfilt -i <input> -d <design> -f <components> -o <out> [options]";
//Command line Options {
Option<string> fnin(string("-i,--in"), string(""),
string(" input file name (4D image)"),
true, requires_argument);
Option<string> fnout(string("-o,--out"), string(""),
- string(" output file name for the filtered data"),
+ string("output file name for the filtered data"),
true, requires_argument);
Option<string> fndesign(string("-d,--design"), string(""),
string("file name of the GLM design matrix (time courses)"),
@@ -42,7 +42,7 @@ using namespace std;
string("mask image file name"),
false, requires_argument);
Option<string> filter(string("-f,--filter"),string(""),
- string("filter out part of the regression model"),
+ string("filter out part of the regression model, e.g. -f \"1,2,3\""),
true, requires_argument);
Option<bool> perfvn(string("--vn"),FALSE,
string(" perfrom variance-normalisation on data"),
diff --git a/meldata.cc b/meldata.cc
index 4cbb2ba..1cb150b 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -62,6 +62,9 @@ namespace Melodic{
message(" done" << endl);
}
+ meanC = mean(tmpData,2);
+ tmpData = remmean(tmpData,2);
+
//convert to power spectra
if(opts.pspec.value()){
message(" Converting data to powerspectra ...");
@@ -69,8 +72,6 @@ namespace Melodic{
message(" done" << endl);
}
- meanC = mean(tmpData,2);
-
//switch dimension in case temporal ICA is required
if(opts.temporal.value()){
message(string(" Switching dimensions for temporal ICA") << endl);
@@ -85,7 +86,7 @@ namespace Melodic{
//variance - normalisation
if(opts.varnorm.value()){
message(" Normalising by voxel-wise variance ...");
- stdDev = varnorm(tmpData,tmpData.Nrows(),3.1);
+ stdDev = varnorm(tmpData,tmpData.Nrows(),2.3);
stdDevi = pow(stdDev,-1);
message(" done" << endl);
}
@@ -171,38 +172,52 @@ namespace Melodic{
opts.approach.set_T("tica");
Matrix alldat, tmpData;
+// bool tmpvarnorm = opts.varnorm.value();
+// opts.varnorm.set_T(false);
alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
for(int ctr = 1; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
alldat += tmpData / numfiles;
}
+ // opts.varnorm.set_T(tmpvarnorm);
+ //variance - normalisation
+ // if(opts.varnorm.value()){
+ // message(" Normalising by voxel-wise variance ...");
+ // stdDev = varnorm(alldat,alldat.Nrows(),3.1);
+ // stdDevi = pow(stdDev,-1);
+ // message(" done" << endl);
+ // }
+
//update mask
if(opts.update_mask.value()){
- message(" Excluding voxels with constant value ...");
- update_mask(Mask, alldat);
+ message("Excluding voxels with constant value ...");
+ // update_mask(Mask, alldat);
message(" done" << endl);
}
- message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
+ if(numfiles>1)
+ message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
//estimate model order
- ColumnVector PPCA;
+ ColumnVector tmpPPCA;
RowVector AdjEV, PercEV;
Matrix Corr, tmpE;
int order;
+ order = ppca_dim(alldat, RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
+
if(opts.pca_dim.value() == 0){
- order = ppca_dim(alldat, RXweight, PPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
- calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
opts.pca_dim.set_T(order);
- }
- else{
- order = opts.pca_dim.value();
- std_pca(alldat, RXweight, Corr, pcaE, pcaD);
- calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
- }
+ PPCA=tmpPPCA;
+ }
+ order = opts.pca_dim.value();
+ calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
+ EV = AdjEV;
+ EVP = PercEV;
+
+
if(numfiles == 1){
Data = alldat;
Matrix tmp = Identity(Data.Nrows());
@@ -212,16 +227,19 @@ namespace Melodic{
else {
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
-
+// if(opts.varnorm.value())
+// varnorm(tmpData,stdDev);
+
// whiten (separate / joint)
if(!opts.joined_whiten.value()){
+ message(" Individual whitening in a " << order << " dimensional subspace " << endl);
std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
}
tmpData = whiteMatrix * tmpData;
DWM.push_back(dewhiteMatrix);
WM.push_back(whiteMatrix);
-
+
//concatenate Data
if(Data.Storage() == 0)
Data = tmpData;
@@ -450,10 +468,21 @@ namespace Melodic{
}while(p);
message(endl);
Matrix newData;
+
+ outMsize("DATA",Data);
+ outMsize("IC",IC);
+ outMsize("noiseIC",noiseIC);
+ outMsize("noiseMix",noiseMix);
+ outMsize("meanR",meanR);
+ outMsize("meanC",meanC);
+
newData = Data - noiseMix * noiseIC.t();
- newData = newData + meanC*ones(1,newData.Ncols());
- newData = newData + ones(newData.Nrows(),1)*meanR;
+ if(meanC.Storage()>0)
+ newData = newData + meanC*ones(1,newData.Ncols());
+ if(meanR.Storage()>0)
+ newData = newData + ones(newData.Nrows(),1)*meanR;
+ cerr << " HERE REMOVE " << endl;
volume4D<float> tmp;
read_volume4D(tmp,opts.inputfname.value().at(0));
tmp.setmatrix(newData,Mask);
@@ -699,14 +728,23 @@ namespace Melodic{
//flip IC maps to be positive (on average)
//flip Subject/Session modes to be positive (on average)
//have time courses accordingly
+
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++)
- if(IC.Row(ctr_i).Sum()<0)
+ if(IC.Row(ctr_i).MaximumAbsoluteValue()>IC.Row(ctr_i).Maximum()){
flipres(ctr_i);
-
- // re-order wrt standard deviation of IC maps
+ }
message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
+ if(numfiles > 1){
+ set_TSmode();
+ Matrix allmodes = Smodes.at(0);
+ for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
+ allmodes |= Smodes.at(ctr);
+ tmpscales = mean(allmodes).t();
+ } else {
+ // re-order wrt standard deviation of IC maps
tmpscales = stdev(IC,2);
+ }
ICstats = tmpscales;
double max_val, min_val = tmpscales.Minimum()-1;
@@ -762,7 +800,6 @@ namespace Melodic{
mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
unmixMatrix = pinv(mixMatrix);
-
}
void MelodicData::status(const string &txt)
diff --git a/meldata.h b/meldata.h
index 623db92..809389e 100644
--- a/meldata.h
+++ b/meldata.h
@@ -158,8 +158,8 @@ namespace Melodic{
inline void set_EV(Matrix& Arg) {if(EV.Storage()==0)
EV = Arg;}
- inline Matrix& get_PPCA() {return PPCA;}
- inline void set_PPCA(Matrix& Arg) {if(PPCA.Storage()==0)
+ inline ColumnVector& get_PPCA() {return PPCA;}
+ inline void set_PPCA(ColumnVector& Arg) {if(PPCA.Storage()==0)
PPCA = Arg;}
inline Matrix& get_stdNoisei() {return stdNoisei;}
@@ -180,7 +180,6 @@ namespace Melodic{
mixMatrix.Column(num) = -mixMatrix.Column(num);
mixFFT=calc_FFT(mixMatrix);
unmixMatrix = pinv(mixMatrix);
-
if(ICstats.Storage()>0&&ICstats.Ncols()>3){
double tmp;
tmp = ICstats(num,3);
@@ -225,7 +224,7 @@ namespace Melodic{
volume<float> background;
Matrix Data;
- Matrix PPCA;
+ ColumnVector PPCA;
Matrix jointCC;
bool after_mm;
diff --git a/melgmix.cc b/melgmix.cc
index 89ed1e5..bbe13dc 100644
--- a/melgmix.cc
+++ b/melgmix.cc
@@ -103,7 +103,7 @@ namespace Melodic{
char *p;
char t[1024];
- const char *discard = ", [];{(})abcdeghijklmopqstuvwxyzABCEGHIJKLMNOQSTUVWXYZ~!@#$%^&*_-=+|\':></?";
+ const char *discard = ", [];{(})abceghijklmoqstuvwxyzABCEGHIJKLMNOQSTUVWXYZ~!@#$%^&*_-=+|\':></?";
strcpy(t, tmpstr.c_str());
p=strtok(t,discard);
@@ -146,7 +146,7 @@ namespace Melodic{
if(fdr.Storage()>0){levls = levls | fdr;}
if(nht.Storage()>0){levls = levls | nht;}
- //cerr << " levles : " << levls << endl << endl;
+ // cerr << " levles : " << levls << endl << endl;
Res = threshold(data, levls);
set_threshmaps(Res);
@@ -186,10 +186,10 @@ namespace Melodic{
add_infstr(" alternative hypothesis test at p > "+float2str(tests(1,4+next),0,2,false));
Matrix tmpNull;
tmpNull = dat;
-
+/*
float cutoffpos, cutoffneg;
- cutoffpos = means(1)+6*std::sqrt(vars(1)+0.0000001);
- cutoffneg = means(1)-6*std::sqrt(vars(1)+0.0000001);
+ cutoffpos = means(1)+100*std::sqrt(vars(1)+0.0000001);
+ cutoffneg = means(1)-100*std::sqrt(vars(1)+0.0000001);
for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
if( probmap(ctr) > tests(1,4+next) ){
@@ -202,7 +202,12 @@ namespace Melodic{
for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
if( (dat(ctr) > cutoffneg) && (dat(ctr) < cutoffpos) )
tmpNull(1,ctr)=0.0;
-
+*/
+ for(int ctr=1; ctr<=tmpNull.Ncols(); ctr++)
+ if( probmap(ctr) < tests(1,4+next) ){
+ tmpNull(1,ctr)=0.0;
+ }
+
Res.Row(next) << tmpNull;
}
next++;
@@ -225,23 +230,36 @@ namespace Melodic{
for(int ctr1=1;ctr1<=tests(1,3);ctr1++){
if(4+next <=tests.Ncols()){
- cerr << " false discovery rate " << tests(1,4+next)<< endl;
- Matrix newcdf(Nprobs);
- for(int ctr=1;ctr<=Nprobs.Nrows();ctr++){
- newcdf.Row(ctr) << 1-props(ctr)*normcdf(dat,means(ctr),vars(ctr));
- }
- Matrix tmpNull;
- tmpNull = newcdf.Row(1);
- Matrix tmpAlt;
- tmpAlt = sum(newcdf.Rows(2,Nprobs.Nrows()),1);
- Matrix tmp;
- tmp = dat;
- for(int ctr=1; ctr<=tmp.Ncols(); ctr++)
- if(tmpNull(1,ctr) > (1-tests(1,4+next)) * tmpAlt(1,ctr))
- tmp(1,ctr)=0.0;
+ message(" Local False Discovery Rate control at p < " << tests(1,4+next) << endl);
+ add_infstr(" Local False Discovery Rate control at p < "+float2str(tests(1,4+next),0,2,false));
+ RowVector tmp=dat;
+ SortAscending(tmp);
+ RowVector newcdf(tmp);
+ newcdf << normcdf(tmp,means(1),vars(1));
+ write_ascii_matrix(dat,"dat");
+ write_ascii_matrix(newcdf,"newcdf");
+ write_ascii_matrix(tmp,"tmpmat");
+ float thrp = tmp(tmp.Ncols())+0.01;
+ float thrn = tmp(1)-0.01;
+ int ctr=tmp.Ncols();
+ do{
+ thrp = tmp(ctr);
+ ctr-=1;
+ }while(ctr>0 && ( (1.0-newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*(tmp.Ncols()-ctr+1)) ));
+
+ ctr=1;
+ do{
+ thrn = tmp(ctr);
+ ctr+=1;
+ }while(ctr<=tmp.Ncols() && ( (newcdf(ctr))*tmp.Ncols() < (tests(1,4+next)*ctr)));
+
+ tmp = dat;
+ for(ctr=1; ctr<=tmp.Ncols();ctr++)
+ if((tmp(ctr) < thrp)&&(tmp(ctr) > thrn))
+ tmp(ctr) = 0.0;
Res.Row(next) << tmp;
- next++;
}
+ next++;
}
for(int ctr1=1;ctr1<=tests(1,4);ctr1++){
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 05a73d4..eaa1d0e 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -121,6 +121,12 @@ namespace Melodic{
return out;
} //RowVector varnorm
+ void varnorm(Matrix& in, RowVector& vars)
+ {
+ Matrix tmp;
+ in = SP(in,pow(ones(in.Nrows(),1) * vars,-1));
+ }
+
RowVector varnorm(Matrix& in, Matrix& Corr, int dim, float level)
{
Matrix tmpE, white, dewhite;
@@ -132,15 +138,15 @@ namespace Melodic{
Matrix ws = white * in;
for(int ctr1 = 1; ctr1<=ws.Ncols(); ctr1++)
for(int ctr2 = 1; ctr2<=ws.Nrows(); ctr2++)
- if(std::abs(ws(ctr2,ctr1)) < level)
- ws(ctr2,ctr1)=0;
+ if(std::abs(ws(ctr2,ctr1)) < level)
+ ws(ctr2,ctr1)=0;
tmpD = stdev(in - dewhite*ws);
for(int ctr1 = 1; ctr1<=tmpD.Ncols(); ctr1++)
if(tmpD(ctr1) < 0.01){
- tmpD(ctr1) = 1.0;
- in.Column(ctr1) = 0.0*in.Column(ctr1);
+ tmpD(ctr1) = 1.0;
+ in.Column(ctr1) = 0.0*in.Column(ctr1);
}
- in = SP(in,pow(ones(in.Nrows(),1) * tmpD,-1));
+ varnorm(in,tmpD);
return tmpD;
} //RowVector varnorm
diff --git a/melhlprfns.h b/melhlprfns.h
index 6f8d5a4..eb1327e 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -29,6 +29,7 @@ namespace Melodic{
Matrix convert_to_pbsc(Matrix& Mat);
RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6);
+ void varnorm(Matrix& in, RowVector& vars);
RowVector varnorm(Matrix& in, Matrix& Corr, int dim = 30, float level = 1.6);
Matrix SP2(const Matrix& in, const Matrix& weights, bool econ = 0);
diff --git a/melodic.cc b/melodic.cc
index 70808e5..4837a28 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -159,6 +159,8 @@ int main(int argc, char *argv[]){
if( bool(opts.genreport.value()) ){
report.analysistxt();
+ if(melodat.get_numfiles()>1)
+ report.Smode_rep();
report.PPCA_rep();
}
diff --git a/melodic.h b/melodic.h
index 1d47215..949a344 100644
--- a/melodic.h
+++ b/melodic.h
@@ -34,8 +34,23 @@
namespace Melodic{
-const string version = "3.0 beta";
-
+const string version = "3.02";
+
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+const string title=string("MELODIC (Version ")+version+")"+
+ string("\n Multivariate Exploratory Linear Optimised Decomposition into Independent Components\n")+
+ string(" Copyright(c) 2007, University of Oxford (Christian F. Beckmann)");
+
+const string usageexmpl=string(" melodic -i <filename> <options>")+
+ string("\n \t \t to run melodic")+
+ // string("\n melodic -i <filename> --mix=melodic_mix")+
+ // string(" --filter=\"string of component numbers\"")+
+ // string("\n \t \t to remove estimated ICs from input")+
+ string("\n melodic -i <filename> --ICs=melodic_IC")+
+ string(" --mix=melodic_mix <options>")+
+ string("\n \t \t to run Mixture Model based inference on estimated ICs")+
+ string("\n melodic --help ");
}
#endif
diff --git a/meloptions.cc b/meloptions.cc
index 5742c17..758f163 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -200,7 +200,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
void MelodicOptions::print_usage(int argc, char *argv[]){
- print_copyright();
+ // print_copyright();
options.usage();
/* cout << "Usage: " << argv[0] << " ";
Have own usage output here
@@ -208,13 +208,12 @@ void MelodicOptions::print_usage(int argc, char *argv[]){
}
void MelodicOptions::print_version(){
- cout << endl <<"MELODIC Version " << version;
+ cout << endl <<"MELODIC Version " << version << endl;
cout.flush();
}
void MelodicOptions::print_copyright(){
- print_version();
- cout << endl << "Copyright (C) 1999-2007 University of Oxford " << endl;
+ cout << endl << title << endl;
cout.flush();
}
diff --git a/meloptions.h b/meloptions.h
index 3ca7a3f..1a6913b 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -133,6 +133,7 @@ class MelodicOptions {
void print_version();
void print_copyright();
void status();
+
};
inline MelodicOptions& MelodicOptions::getInstance(){
@@ -333,16 +334,7 @@ class MelodicOptions {
string("perform temporal ICA"),
false, no_argument, false),
retrystep(3),
- options(string(""),
- string(" melodic -i <filename> <options>")+
- string("\n \t \t to run melodic")+
-// string("\n melodic -i <filename> --mix=melodic_mix")+
-// string(" --filter=\"string of component numbers\"")+
-// string("\n \t \t to remove estimated ICs from input")+
- string("\n melodic -i <filename> --ICs=melodic_IC")+
- string(" --mix=melodic_mix <options>")+
- string("\n \t \t to run Mixture Model based inference on estimated ICs")+
- string("\n melodic --help "))
+ options(title, usageexmpl)
{
try {
options.add(logdir);
diff --git a/melreport.cc b/melreport.cc
index a61db40..fc819f0 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -73,11 +73,11 @@ namespace Melodic{
miscpic newpic;
float map1min = std::max((map1 + binarise(tempVol[0],tempVol[0].min(),
- float(0.0)) * map1.max()).min(),float(0.01));
- float map1max = std::max(map1.max(),float(0.01));
- float map2min = std::min(map2.min(),float(-0.01));
+ float(0.0)) * map1.max()).min(),float(0.001));
+ float map1max = std::max(map1.max(),float(0.001));
+ float map2min = std::min(map2.min(),float(-0.001));
float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
- tempVol[0].max()) * map2.min()).max(),float(-0.01));
+ tempVol[0].max()) * map2.min()).max(),float(-0.001));
newpic.overlay(newvol, melodat.get_bg(), map1, map2,
melodat.get_bg().percentile(0.02),
@@ -96,12 +96,11 @@ namespace Melodic{
newpic.set_title(string("Component No. "+num2str(cnum)+
" - thresholded IC map ") + mmodel.get_infstr(1));
newpic.set_cbar(string("ysb"));
- //cerr << instr << endl;
- if(map1.max()-map1.min()>0.01)
+ if((std::abs(map1.max()-map1.min())>0.01) || (std::abs(map2.max()-map2.min())>0.01))
newpic.slicer(newvol, instr, &melodat.tempInfo);
else
- newpic.slicer(map1, instr, &melodat.tempInfo);
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
@@ -489,7 +488,7 @@ namespace Melodic{
report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
string(mmodel.get_prefix() +
" GGM("+num2str(mmodel.mixtures())+") fit"),
- true, float(0.0), float(2.0));
+ true, float(0.0), float(0.0));
}
else{
@@ -693,6 +692,7 @@ namespace Melodic{
newplot.add_label("dim. estimate");
}
+ newplot.grid_swapdefault();
newplot.timeseries(what,
report.appendDir("EVplot.png"),
string("Eigenspectrum Analysis"),
@@ -701,4 +701,18 @@ namespace Melodic{
report << "<img BORDER=0 SRC=\"EVplot.png\"><p>" << endl;
}//time series plot
}
+
+ void MelodicReport::Smode_rep(){
+ report << "<p> <H3>TICA Subject/Session modes </H3> <p>" << endl;
+ miscplot newplot;
+ Matrix allmodes = melodat.get_Smodes().at(0);
+ for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr)
+ allmodes |= melodat.get_Smodes().at(ctr);
+
+ newplot.set_xysize(100+30*allmodes.Ncols(),300);
+ newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")),
+ string("Subject/Session modes"));
+
+ report << "<img BORDER=0 SRC=\"bp_Smodes.png\"><p>" << endl;
+ }
}
diff --git a/melreport.h b/melreport.h
index 68662e9..dff72ce 100644
--- a/melreport.h
+++ b/melreport.h
@@ -186,6 +186,7 @@ namespace Melodic{
void IC_simplerep(string prefix, int cnum, int dim);
void PPCA_rep();
+ void Smode_rep();
private:
MelodicData &melodat;
diff --git a/test.cc b/test.cc
index d0e0165..33ff5a6 100644
--- a/test.cc
+++ b/test.cc
@@ -80,14 +80,8 @@ int do_work(int argc, char* argv[]) {
cerr << weights << endl <<endl << Rows << endl;
- cerr << SP(weights,pow(Rows,-1)) << endl;
+ cerr << remmean(SP(weights,pow(Rows,-1))) << endl;
- vector<int> vi;
- for(int ctr=1;ctr<=5;ctr++)
- vi.push_back(ctr);
-
- for(int ctr=0;ctr<5;ctr++)
- cerr << " vi["<<ctr<<"] = " <<vi.at(ctr)<<endl;
return 0;
}
--
GitLab