fixed ard bug

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;
     }
 

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

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);

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);