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

Melodic 3.03

parent b0075278
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 1
− 1
View file @ f5c6cf39
......@@ -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;
......
melhlprfns.cc
+ 4
− 8
View file @ f5c6cf39
......@@ -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(),
......
melodic.cc
+ 1
− 0
View file @ f5c6cf39
......@@ -10,6 +10,7 @@
/* CCOPYRIGHT */
#include <iostream>
#include <iomanip>
#include "newmatap.h"
#include "newmatio.h"
#include "newimage/newimageall.h"
......
melodic.h
+ 1
− 1
View file @ f5c6cf39
......@@ -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
......
meloptions.cc
+ 1
− 0
View file @ f5c6cf39
......@@ -11,6 +11,7 @@
/* CCOPYRIGHT */
#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
......
meloptions.h
+ 1
− 0
View file @ f5c6cf39
......@@ -17,6 +17,7 @@
#include <string>
#include <strstream>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
......
melreport.cc
+ 103
− 68
View file @ f5c6cf39
......@@ -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 << "; &nbsp; &nbsp; " << 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 << " &beta;(" <<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 &beta;'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"),
......
melreport.h
+ 2
− 2
View file @ f5c6cf39
......@@ -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;
}
......
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