Skip to content
Snippets Groups Projects
Commit 16cccdfb authored by Christian Beckmann's avatar Christian Beckmann
Browse files

updated to ln(12) beta: percent explained variance etc

parent f7ba695d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 95
− 1
View file @ 16cccdfb
......@@ -129,7 +129,8 @@ namespace Melodic{
}
}
stdDevi = pow(stdev(Data - tmpDeWhite*WS),-1);
stdDev = stdev(Data - tmpDeWhite*WS);
stdDevi = pow(stdDev,-1);
Data = Data + meanC*ones(1,Data.Ncols());
}
......@@ -253,6 +254,14 @@ namespace Melodic{
logger.appendDir(opts.outputfname.value() + "_FTmix") <<endl);
}
//Output ICstats
if(ICstats.Storage()>0){
write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_ICstats"),
ICstats);
message(" "<<
logger.appendDir(opts.outputfname.value() + "_ICstats") <<endl);
}
//Output unmixMatrix
if(opts.output_unmix.value() && unmixMatrix.Storage()>0){
write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_unmix"),unmixMatrix);
......@@ -457,6 +466,91 @@ namespace Melodic{
}
} //void create_mask()
void MelodicData::sort()
{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
numTime = mixMatrix.Nrows(), i,j;
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
if(IC.Row(ctr_i).Sum()>0){
flipres(ctr_i); };}
// cerr << "HERE2" << endl << endl;
// re-order wrt standard deviation of IC maps
message(" sorting IC maps" << endl << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
tmpscales = stdev(IC,2);
ICstats = tmpscales;
//cerr << "SCLAES2 " << tmpscales << endl;
double max_val, min_val = tmpscales.Minimum()-1;
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
max_val = tmpscales.Maximum2(i,j);
ICstats(ctr_i,1)=max_val;
//cerr << endl << "scales: " << tmpscales << " ctr_i " << ctr_i <<
// " i " << i << " j " << j << " max " << max_val << endl << endl;
tmpICrow = IC.Row(ctr_i);
tmpMIXcol = mixMatrix.Column(ctr_i);
IC.SubMatrix(ctr_i,ctr_i,1,numVox) = IC.SubMatrix(i,i,1,numVox);
mixMatrix.SubMatrix(1,numTime,ctr_i,ctr_i) =
mixMatrix.SubMatrix(1,numTime,i,i);
IC.SubMatrix(i,i,1,numVox) = tmpICrow.SubMatrix(1,1,1,numVox);
mixMatrix.SubMatrix(1,numTime,i,i) = tmpMIXcol.SubMatrix(1,numTime,1,1);
tmpscales(i,1)=tmpscales(ctr_i,1);
tmpscales(ctr_i,1)=min_val;
}
//cerr << " ICstats " << ICstats << endl << endl;
ICstats /= ICstats.Column(1).Sum();
ICstats *= 100;
//cerr << " ICstats " << ICstats << endl << endl;
if(EVP.Storage()>0){
tmpscales = ICstats.Column(1).AsMatrix(ICstats.Nrows(),1) * EVP(1,numComp);
ICstats |= tmpscales;
}
if(DataVN.Storage()>0&&stdDev.Storage()>0){
//cerr << " ICstats " << ICstats << endl << endl;
Matrix copeP(tmpscales), copeN(tmpscales);
Matrix max_ICs(tmpscales), min_ICs(tmpscales);
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
int i,j;
max_ICs(ctr_i,1) = IC.Row(ctr_i).Maximum2(i,j);
//cerr << " ICstats " << ICstats << endl << endl;
//cerr << endl <<(pinv(mixMatrix)*DataVN.Column(j)) << endl;
copeP(ctr_i,1) = std::abs((pinv(mixMatrix)*DataVN.Column(j)).Row(ctr_i).AsScalar()*stdDev(1,j)*100*(mixMatrix.Column(ctr_i).Maximum()-mixMatrix.Column(ctr_i).Minimum())/meanR(1,j));
min_ICs(ctr_i,1) = IC.Row(ctr_i).Minimum2(i,j);
copeN(ctr_i,1) = -1.0*std::abs((pinv(mixMatrix)*DataVN.Column(j)).Row(ctr_i).AsScalar()*stdDev(1,j)*100*(mixMatrix.Column(ctr_i).Maximum()-mixMatrix.Column(ctr_i).Minimum())/meanR(1,j));
}
ICstats |= copeP;
ICstats |= copeN;
}
mixFFT=calc_FFT(mixMatrix);
unmixMatrix = pinv(mixMatrix);
//if(ICstats.Storage()>0){cout << "ICstats: " << ICstats.Nrows() <<"x" << ICstats.Ncols() << endl;}else{cout << "ICstats empty " <<endl;}
}
void MelodicData::status(const string &txt)
{
......
meldata.h
+ 14
− 0
View file @ 16cccdfb
......@@ -101,6 +101,9 @@ namespace Melodic{
inline Matrix& get_RXweight() {return RXweight;}
inline void set_RXweight(Matrix& Arg) {RXweight = Arg;}
inline Matrix& get_ICstats() {return ICstats;}
inline void set_ICstats(Matrix& Arg) {ICstats = Arg;}
inline Matrix& get_EVP() {return EVP;}
inline void set_EVP(Matrix& Arg) {EVP = Arg;}
......@@ -128,6 +131,13 @@ 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);
ICstats(num,3) = -1.0*ICstats(num,4);
ICstats(num,4) = -1.0*tmp;
}
}
......@@ -135,6 +145,8 @@ namespace Melodic{
Data = SP(Data,ones(Data.Nrows(),1)*stdDevi);
}
void sort();
volumeinfo tempInfo;
Matrix calc_FFT(const Matrix& Mat);
......@@ -150,12 +162,14 @@ namespace Melodic{
Matrix dewhiteMatrix;
Matrix meanC;
Matrix meanR;
Matrix stdDev;
Matrix stdDevi;
Matrix RXweight;
Matrix mixMatrix;
Matrix unmixMatrix;
Matrix mixFFT;
Matrix IC;
Matrix ICstats;
Matrix EVP;
Matrix EV;
Matrix stdNoisei;
......
melica.cc
+ 16
− 21
View file @ 16cccdfb
......@@ -437,26 +437,6 @@ write_ascii_matrix("Xim",Data.Row(1));
return Res;
}
void MelodicICA::sort()
{
int numComp = melodat.get_mix().Ncols();
Matrix tmpIC = melodat.get_IC();
Matrix tmpA = melodat.get_mix();
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
if(tmpIC.Row(ctr_i).Sum()>0){
tmpIC.Row(ctr_i) = -tmpIC.Row(ctr_i);
tmpA.Column(ctr_i) = -tmpA.Column(ctr_i);
};
}
melodat.set_mix(tmpA);
melodat.set_IC(tmpIC);
Matrix tmpW = pinv(tmpA);
melodat.set_unmix(tmpW);
}
void MelodicICA::perf_ica(const Matrix &Data)
{
message("Starting ICA estimation using " << opts.approach.value()
......@@ -482,8 +462,23 @@ write_ascii_matrix("Xim",Data.Row(1));
// Add the mean time course again
temp += melodat.get_unmix()*melodat.get_meanC()*ones(1,temp.Ncols());
//re-normalise the decomposition to std(mix)=1
Matrix scales;
scales = stdev(melodat.get_mix());
// cerr << " SCALES 1 " << scales << endl;
Matrix tmp;
tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
temp = SP(temp,scales.t()*ones(1,temp.Ncols()));
scales = scales.t();
melodat.set_mix(tmp);
melodat.set_IC(temp);
sort();
melodat.set_ICstats(scales);
melodat.sort();
}
}
}
......
melodic.cc
+ 2
− 10
View file @ 16cccdfb
......@@ -135,15 +135,6 @@ int main(int argc, char *argv[])
if(!no_conv){
//re-normalise the decomposition to std(mix)=1
Matrix scales;
scales = stdev(melodat.get_mix());
Matrix tmp;
tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
melodat.set_mix(tmp);
tmp = SP(melodat.get_IC(),scales.t()*ones(1,melodat.get_IC().Ncols()));
melodat.set_IC(tmp);
//save raw IC results
melodat.save();
......@@ -269,6 +260,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
melodat.set_fmix(fmixMatrix);
fmixMatrix = pinv(mixMatrix);
melodat.set_unmix(fmixMatrix);
melodat.sort();
// write_ascii_matrix("ICs",ICs);
......@@ -414,7 +406,7 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
if( bool(opts.genreport.value()) ){
message(" creating report page ... ");
report.IC_rep(mixmod,ctr,melodat.get_IC().Nrows());
report.IC_rep(mixmod,ctr,melodat.get_IC().Nrows(),melodat.get_ICstats());
message(" done" << endl);
}
}
......
melodic.h
+ 1
− 1
View file @ 16cccdfb
......@@ -32,7 +32,7 @@
namespace Melodic{
const string version = "ln(11)";
const string version = "ln(12) beta";
}
......
melreport.cc
+ 13
− 2
View file @ 16cccdfb
......@@ -15,7 +15,7 @@
namespace Melodic{
void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim)
void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats)
{
if( bool(opts.genreport.value()) ){
addlink(mmodel.get_prefix()+".html",num2str(cnum));
......@@ -44,7 +44,18 @@ namespace Melodic{
IChtml << "<hr><p>" << endl;
}
{//output IC stats
if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
IChtml << ICstats(cnum,1) << " % of explained variance";
if(ICstats.Ncols()>1)
IChtml << "; &nbsp; &nbsp; " << ICstats(cnum,2) << " % of total variance";
if(ICstats.Ncols()>2){
IChtml << "<p>" <<endl;
IChtml << " &nbsp; &nbsp; " << ICstats(cnum,3) << " % signal change (pos peak voxel); &nbsp; &nbsp;" << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;}
IChtml << "<hr><p>" << endl;
}
}
volume4D<float> tempVol;
if(mmodel.get_threshmaps().Storage()>0&&
......
melreport.h
+ 1
− 1
View file @ 16cccdfb
......@@ -202,7 +202,7 @@ namespace Melodic{
return report.getDir();
}
void IC_rep(MelGMix &mmodel, int cnum, int dim);
void IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats);
void IC_simplerep(string prefix, int cnum, int dim);
void PPCA_rep();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment