From f4af04864c524eae190b5d5e44c6d619b8ad1b58 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Wed, 26 Jul 2006 12:43:42 +0000
Subject: [PATCH] fixed ard bug
---
fibre.h | 27 +++++++++++++++++++++------
old_bedpostX_proc | 4 ++--
xfibres.cc | 10 ++++++++--
xfibresoptions.h | 4 ++++
4 files changed, 35 insertions(+), 10 deletions(-)
diff --git a/fibre.h b/fibre.h
index 38005af..6a32207 100644
--- a/fibre.h
+++ b/fibre.h
@@ -90,10 +90,12 @@ namespace FIBRE{
compute_ph_prior();
m_f_prior=0;
- compute_f_prior();
-
- m_lam_prior=0;
- compute_lam_prior();
+ compute_f_prior();
+ //cc OUT("COCK1");
+ //cc OUT(m_f_prior);
+ //cc OUT(m_ardfudge);
+ m_lam_prior=0;
+ compute_lam_prior();
m_Signal.ReSize(alpha.Nrows());
m_Signal=0;
@@ -132,7 +134,9 @@ namespace FIBRE{
m_f_prior=0;
compute_f_prior();
-
+ //cc OUT("COCK2");
+ //cc OUT(m_f_prior);
+ //cc OUT(m_ardfudge);
m_lam_prior=0;
compute_lam_prior();
@@ -269,13 +273,23 @@ namespace FIBRE{
m_f_prior=0;
}
else{
- if(m_lam_jump)
+ if(m_lam_jump){
// m_f_prior=log(1-m_f)+2*log(fabs(log(1-m_f))); //marginalised with uniform prior on lambda
m_f_prior=std::log(m_f);
+
+ //cc OUT(m_f);
+ //cc OUT(m_ardfudge);
+ //cc float mmk=m_ardfudge*m_f_prior;
+ //cc OUT(mmk);
+ //cc OUT(m_f_old_prior);
+
+ }
else
m_f_prior=0;
+
}
m_f_prior=m_ardfudge*m_f_prior;
+
return false;
}
}
@@ -454,6 +468,7 @@ namespace FIBRE{
m_lam_jump=rhs.m_lam_jump;
m_Signal=rhs.m_Signal;
m_Signal_old=rhs.m_Signal_old;
+ m_ardfudge=rhs.m_ardfudge;
return *this;
}
diff --git a/old_bedpostX_proc b/old_bedpostX_proc
index c9868ef..4cdbd5f 100644
--- a/old_bedpostX_proc
+++ b/old_bedpostX_proc
@@ -28,8 +28,8 @@ while [ $slice -lt $nslices ];do
echo `hostname`_${$} > ${subjdir}.bedpostX/logs/.$slicezp
sleep 10
if [ `hostname`_${$} = `cat ${subjdir}.bedpostX/logs/.$slicezp | sed -n '1p'` ] ; then
- nice ${FSLDIR}/bin/xfibres --data=$subjdir/data_slice_$slicezp --mask=$subjdir/nodif_brain_mask_slice_$slicezp -b $subjdir/bvals -r $subjdir/bvecs --forcedir --logdir=$subjdir.bedpostX/diff_slices/data_slice_$slicezp --nj=1000 --bi=2000 --bn=0 --se=20 --nfibres=2 --fudge=1.5> $subjdir.bedpostX/logs/log$slicezp
-
+ nice ${FSLDIR}/bin/xfibres --data=$subjdir/data_slice_$slicezp --mask=$subjdir/nodif_brain_mask_slice_$slicezp -b $subjdir/bvals -r $subjdir/bvecs --forcedir --logdir=$subjdir.bedpostX/diff_slices/data_slice_$slicezp --nj=1000 --bi=600 --bn=0 --se=20 --nfibres=2 --fudge=1.5> $subjdir.bedpostX/logs/log$slicezp
+#bi can be changed to 2000 if convergence looks dodgy
touch ${subjdir}.bedpostX/logs/.${slicezp}_finished
fi
fi
diff --git a/xfibres.cc b/xfibres.cc
index a8a8573..b4da466 100644
--- a/xfibres.cc
+++ b/xfibres.cc
@@ -300,7 +300,6 @@ public:
}
}
ret=(maxfib>1); //
- cout<<ret<<" ret"<<endl;
if(ret){
mfibre.set_d(m_mean_dsamples(voxbest));
mfibre.set_S0(m_mean_S0samples(voxbest));
@@ -441,9 +440,13 @@ class xfibresVoxelManager{
void initialise(const Matrix& Amat){
+ if(!opts.localinit.value()){
if(!m_samples.neighbour_initialise(m_voxelnumber,m_multifibre)){
initialise_tensor(Amat);
}
+ }else{
+ initialise_tensor(Amat);
+ }
m_multifibre.initialise_energies();
m_multifibre.initialise_props();
}
@@ -465,14 +468,16 @@ class xfibresVoxelManager{
logS(i)=0;
}
}
+
Dvec = -pinv(Amat)*logS;
+
if( Dvec(7) > -maxlogfloat ){
S0=exp(-Dvec(7));
}
else{
S0=m_data.MaximumAbsoluteValue();
}
-
+
for ( int i = 1; i <= logS.Nrows(); i++)
{
if(S0<m_data.Sum()/m_data.Nrows()){ S0=m_data.MaximumAbsoluteValue(); }
@@ -481,6 +486,7 @@ class xfibresVoxelManager{
Dvec = -pinv(Amat)*logS;
S0=exp(-Dvec(7));
+
if(S0<m_data.Sum()/m_data.Nrows()){ S0=m_data.Sum()/m_data.Nrows(); }
tens = vec2tens(Dvec);
EigenValues(tens,Dd,Vd);
diff --git a/xfibresoptions.h b/xfibresoptions.h
index f70c3a6..f74b8e0 100644
--- a/xfibresoptions.h
+++ b/xfibresoptions.h
@@ -44,6 +44,7 @@ class xfibresOptions {
Option<int> updateproposalevery;
Option<int> seed;
Option<bool> no_ard;
+ Option<bool> localinit;
void parse_command_line(int argc, char** argv, Log& logger);
private:
@@ -112,6 +113,8 @@ class xfibresOptions {
false,requires_argument),
no_ard(string("--noard"),false,string("Turn ARD off on all fibres"),
false,no_argument),
+ localinit(string("--nospat"),false,string("Initialise with tensor, not spatially"),
+ false,no_argument),
options("xfibres", "xfibres -k <filename>\n xfibres --verbose\n")
{
@@ -130,6 +133,7 @@ class xfibresOptions {
options.add(njumps);
options.add(nburn);
options.add(nburn_noard);
+ options.add(localinit);
options.add(sampleevery);
options.add(updateproposalevery);
options.add(seed);
--
GitLab