diff --git a/meldata.cc b/meldata.cc
index 582434954a3ae25f3f7c6fc819ce2f93f7a43660..fb227a159ffbac56b613707b6141e6d690616f3b 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -341,7 +341,6 @@ namespace Melodic{
SconF = read_ascii_matrix(opts.fn_SconF.value());
Tdes = remmean(Tdes,1);
- Sdes = remmean(Sdes,1);
}
void MelodicData::save()
@@ -781,6 +780,7 @@ namespace Melodic{
// re-order wrt standard deviation of IC maps
tmpscales = stdev(IC,2);
}
+
ICstats = tmpscales;
double max_val, min_val = tmpscales.Minimum()-1;
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 95666b57a75c397d1981b62a0f40a0d773b40b06..e353681668569dc8d8b7623b582dbedf32c382e1 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -805,23 +805,19 @@ namespace Melodic{
dof = (int)-1; cbeta = -1.0*ones(1);
if(data.Nrows()==design.Nrows()){
- Matrix dat = remmean(data,1);
+ Matrix dat = data;
Matrix tmp = design.t()*design;
Matrix pinvdes = tmp.i()*design.t();
beta = pinvdes * dat;
residu = dat - design*beta;
-// Matrix R = Identity(design.Nrows()) - design * pinvdes;
-// Matrix R2 = R*R;
-// float tR = R.Trace();
-// sigsq = sum(SP(residu,residu))/tR;
-
-// dof = (int)(tR*tR/R2.Trace() - DOFadjust);
+
dof = design.Nrows() - design.Ncols();
sigsq = sum(SP(residu,residu))/dof;
float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
- f_fmf = SP(var(design*beta),pow(var(residu),-1))* fact;
+ f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
+
pf_fmf = f_fmf.Row(1);
for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
diff --git a/melodic.cc b/melodic.cc
index 63d3132dff8fb9260708536c958903bd4981e0c1..670f4bb51da421e5c31f092fca970a3307f9c2eb 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -10,6 +10,7 @@
/* CCOPYRIGHT */
#include <iostream>
+#include <iomanip>
#include "newmatap.h"
#include "newmatio.h"
#include "newimage/newimageall.h"
diff --git a/melodic.h b/melodic.h
index 16d8516d4eb01b07d7abd9b53f16043e37d60f4e..0eefaffbfeb6e9cdb4d1d6d3a1d5549527763469 100644
--- a/melodic.h
+++ b/melodic.h
@@ -40,7 +40,7 @@
namespace Melodic{
-const string version = "3.02";
+const string version = "3.03";
// The two strings below specify the title and example usage that is
// printed out as the help or usage message
diff --git a/meloptions.cc b/meloptions.cc
index 758f163263a66d107b79e4e7d9e05fefc4f857b9..6d50d52b70524d8bbe542a551a02acb8487914e1 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -11,6 +11,7 @@
/* CCOPYRIGHT */
#include <iostream>
+#include <iomanip>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
diff --git a/meloptions.h b/meloptions.h
index a6256aed0b081b3c83eac3725311e05cb11f380b..6ffd12ebc28826af237d76cad44ccab6f9efbfd5 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -17,6 +17,7 @@
#include <string>
#include <strstream>
#include <iostream>
+#include <iomanip>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
diff --git a/melreport.cc b/melreport.cc
index 272598e5ce088275c133439c44fa3dac237e80ab..01cf3aa02f5254180053c502e3647c78725a267c 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -45,7 +45,7 @@ namespace Melodic{
}
{//output IC stats
if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
- IChtml << ICstats(cnum,1) << " % of explained variance";
+ IChtml << fixed << setprecision(2) << ICstats(cnum,1) << " % of explained variance";
if(ICstats.Ncols()>1)
IChtml << "; " << ICstats(cnum,2) << " % of total variance";
if(ICstats.Ncols()>2&&melodat.get_numfiles()==1){
@@ -103,6 +103,7 @@ namespace Melodic{
IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
}
+
{//plot time course
IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
miscplot newplot;
@@ -114,25 +115,32 @@ namespace Melodic{
newplot.add_label(string("IC ")+num2str(cnum)+" time course");
newplot.add_label("full model fit");
}
-
- if(opts.tr.value()>0.0)
- newplot.timeseries(tmptc,
- report.appendDir(string("t")+num2str(cnum)+".png"),
- string("Timecourse (in seconds); TR = ")+
- float2str(opts.tr.value(),0,2,0)+" s",
- opts.tr.value(),150,4,1);
- else
- newplot.timeseries(tmptc,
- report.appendDir(string("t")+num2str(cnum)+".png"),
- string("Timecourse (in TRs)"));
- write_ascii_matrix(report.appendDir(string("t")
- +num2str(cnum)+".txt"),tmptc.t());
- IChtml << "<A HREF=\"" << string("t")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+
+ if(opts.tr.value()>0.0)
+ newplot.add_xlabel(string("Time (seconds); TR = ")+
+ float2str(opts.tr.value(),0,2,0)+" s");
+ else
+ newplot.add_xlabel(string("Time (TRs)"));
+
+ newplot.add_ylabel("Normalised Response");
+ newplot.set_yrange(tmptc.Minimum()-0.05*(tmptc.Maximum()-tmptc.Minimum()),
+ tmptc.Maximum()+0.05*(tmptc.Maximum()-tmptc.Minimum()));
+ newplot.grid_swapdefault();
+ newplot.timeseries(tmptc,
+ report.appendDir(string("t")+num2str(cnum)+".png"),
+ string("Timecourse No. ")+num2str(cnum),
+ opts.tr.value(),150,12,1,false);
+
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),tmptc.t());
+ IChtml << "<A HREF=\"" << string("t")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
}//time series plot
- {//plot frequency
+
+ if(!opts.pspec.value())
+ {//plot frequency
miscplot newplot;
RowVector empty(1);
empty = 0.0;
@@ -144,6 +152,9 @@ namespace Melodic{
Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1), opts.logPower.value());
+ newplot.set_Ylabel_fmt("%.0f");
+ newplot.set_yrange(0.0,1.02*fmixtc.Maximum());
+ newplot.grid_swapdefault();
if(opts.tr.value()>0.0){
newplot.add_xlabel(string("Frequency (in Hz / ")+num2str(fact)+ " )");
newplot.timeseries(empty | fmixtc.t(),
@@ -188,17 +199,17 @@ namespace Melodic{
" <BR> dof1 = " << melodat.Tdes.Ncols() << "; dof2 = "
<< melodat.glmT.get_dof() << "<BR>" <<endl;
if(melodat.glmT.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
- IChtml << "<b> p < " << melodat.glmT.get_pf_fmf().Column(cnum) <<
+ IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmT.get_pf_fmf().Column(cnum) <<
"<BR> (uncorrected for #comp.)<b></TD>" << endl;
else
- IChtml << " p < " <<
+ IChtml << fixed << setprecision(5) << " p < " <<
melodat.glmT.get_pf_fmf().Column(cnum) <<
"<BR> (uncorrected for #comp.)</TD>" << endl;
}
if(melodat.Tcon.Storage() > 0){
- IChtml << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
+ IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
- IChtml << "con(" << melodat.Tcon.Row(ctr) << "): <br>" << endl;
+ IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
IChtml << "<td align=right>" << endl;
for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
IChtml <<" z = <BR>" <<endl;
@@ -208,10 +219,10 @@ namespace Melodic{
IChtml << "<td align=right>" << endl;
for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
if(melodat.glmT.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
- IChtml << "<b> p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
"</b><BR>" << endl;
else
- IChtml << " p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ IChtml << fixed << setprecision(5) <<" p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
"<BR>" << endl;
IChtml << "</TABLE></td></tr>" << endl;
}
@@ -227,18 +238,26 @@ namespace Melodic{
miscplot newplot;
//add GLM OLS fit
+ //newplot.setscatter(smode,2);
+
if(melodat.Sdes.Storage() > 0){
smode |= melodat.Sdes * melodat.glmS.get_beta().Column(cnum);
newplot.add_label(string("IC ")+num2str(cnum)+" subject/session-mode");
newplot.add_label("full model fit");
}
- newplot.setscatter(smode,5);
+ newplot.grid_swapdefault();
+ newplot.set_Ylabel_fmt("%.2f");
+ newplot.add_xlabel(" Subject Number");
+// newplot.set_xysize(smode.Nrows()*80,150);
newplot.timeseries(smode.t(),
report.appendDir(string("s")+num2str(cnum)+".png"),
string("Subject/Session mode"));
+ newplot.clear_xlabel();
+ newplot.clear_labels();
newplot.set_xysize(120,200);
newplot.set_minmaxscale(1.1);
- newplot.boxplot(smode,
+
+ newplot.boxplot((Matrix)smode.Column(1),
report.appendDir(string("b")+num2str(cnum)+".png"),
string("Subject/Session mode"));
write_ascii_matrix(report.appendDir(string("s")
@@ -253,39 +272,52 @@ namespace Melodic{
+string("b")+num2str(cnum)+".png\"></A><p>" << endl;
}
{//add S-mode GLM F-stats for full model fit & contrasts
- if(melodat.Sdes.Storage() > 0){
- IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
- "<CAPTION><EM> <b>GLM (OLS) on subject/session-mode </b></EM></CAPTION>" << endl
- << "<TR valign=middle><TH colspan=2>Betas <TH> <EM> F-test on <br> full model fit </em>";
- if(melodat.Scon.Storage() > 0)
- IChtml << "<TH colspan=3><EM>Contrasts</EM>"<<endl;
- IChtml << "<TR><TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
- IChtml << " β(" <<num2str(ctr)+"): <br>" << endl;
- IChtml << "<TD align=right>" << endl;
- for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
- IChtml << melodat.glmS.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
- IChtml <<
- " <TD align=center> F = "<< melodat.glmS.get_f_fmf().Column(cnum) <<
- " <BR> dof1 = " << melodat.Sdes.Ncols() << "; dof2 = "
- << melodat.glmS.get_dof() << "<BR> p < " <<
- melodat.glmS.get_pf_fmf().Column(cnum) << "</TD>" << endl;
- }
- if(melodat.Scon.Storage() > 0){
- IChtml << "<td>" <<endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml << "con(" << melodat.Scon.Row(ctr) << ") <br>" << endl;
- IChtml << "<td align=center>" << endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml << " z = " << melodat.glmS.get_z().Column(cnum).Row(ctr) <<
- "<BR>" <<endl;
- IChtml << "<td align=right>" << endl;
- for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
- IChtml << " p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
- "<BR>" << endl;
- IChtml << "</td></tr>" << endl;
+ if(melodat.Sdes.Storage() > 0){
+ IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
+ "<CAPTION><EM> <b>GLM (OLS) on subject/session-mode </b></EM></CAPTION>" << endl
+ << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
+ if(melodat.Scon.Storage() > 0)
+ IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
+ IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
+ IChtml << "<TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << melodat.glmS.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
+ IChtml << "</TABLE>" <<
+ " <TD align=center> F = "<< melodat.glmS.get_f_fmf().Column(cnum) <<
+ " <BR> dof1 = " << melodat.Sdes.Ncols() << "; dof2 = "
+ << melodat.glmS.get_dof() << "<BR>" <<endl;
+ if(melodat.glmS.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) <<"<b> p < " << melodat.glmS.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)<b></TD>" << endl;
+ else
+ IChtml << fixed << setprecision(5) << " p < " <<
+ melodat.glmS.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)</TD>" << endl;
+ }
+ if(melodat.Scon.Storage() > 0){
+ IChtml << fixed << setprecision(2) << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << "COPE(" << num2str(ctr) << "): <br>" << endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml <<" z = <BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << melodat.glmS.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ if(melodat.glmS.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
+ IChtml << fixed << setprecision(5) << "<b> p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
+ "</b><BR>" << endl;
+ else
+ IChtml << fixed << setprecision(5) <<" p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
+ "<BR>" << endl;
+ IChtml << "</TABLE></td></tr>" << endl;
}
- IChtml << "</TABLE><p>" << endl;
+ IChtml << "</TABLE><p>" << endl;
+
}
}//subject mode plot
@@ -348,8 +380,8 @@ namespace Melodic{
{ //finish IC page
IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
} //finish IC page
@@ -486,7 +518,7 @@ namespace Melodic{
mmodel.get_pi(),
report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
string(mmodel.get_prefix() +
- " GGM("+num2str(mmodel.mixtures())+") fit"),
+ " Gaussian/Gamma Mixture Model("+num2str(mmodel.mixtures())+") fit"),
true, float(0.0), float(0.0));
}
@@ -500,7 +532,7 @@ namespace Melodic{
mmodel.get_pi(),
report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
string(mmodel.get_prefix() +
- " GMM("+num2str(mmodel.mixtures())+") fit"),
+ " Gaussian Mixture Model("+num2str(mmodel.mixtures())+") fit"),
false, float(0.0), float(2.0));
}
@@ -520,8 +552,8 @@ namespace Melodic{
{ //finish IC2 page
IChtml2<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
} //finish IC2 page
@@ -663,8 +695,8 @@ namespace Melodic{
{ //finish IC page
IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
} //finish IC page
@@ -690,7 +722,10 @@ namespace Melodic{
what &= melodat.get_PPCA().t();
newplot.add_label("dim. estimate");
}
-
+
+ newplot.set_Ylabel_fmt("%.2f");
+ newplot.add_xlabel("Number of included components");
+ newplot.set_yrange(0.0,1.02);
newplot.grid_swapdefault();
newplot.timeseries(what,
report.appendDir("EVplot.png"),
diff --git a/melreport.h b/melreport.h
index dff72ce4e3a43a1557907a8f4d59cb2a0d494d5f..93cd08e0432a96bdb549b028356cdca58a6e622f 100644
--- a/melreport.h
+++ b/melreport.h
@@ -62,8 +62,8 @@ namespace Melodic{
~MelodicReport(){
if( bool(opts.genreport.value()) ){
report << "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC</A> Version "
+ << version << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
}