diff --git a/fibre.h b/fibre.h
index bf0d284cbb897cc3f3c9b3fc86047a7e697db8e0..c23f7b9cebd7c84058677e1bc3c1fd7046a0229c 100644
--- a/fibre.h
+++ b/fibre.h
@@ -23,32 +23,6 @@ const float minlogfloat=-23;
namespace FIBRE{
- //Returns the natural log of the 0th order modified Bessel function of first kind for an argument x
- //Follows the exponential implementation of the Bessel function in Numerical Recipes, Ch. 6
- float logIo(const float& x){
- float y,b;
-
- b=std::fabs(x);
- if (b<3.75){
- float a=x/3.75;
- a*=a;
- //Bessel function evaluation
- y=1.0+a*(3.5156229+a*(3.0899424+a*(1.2067492+a*(0.2659732+a*(0.0360768+a*0.0045813)))));
- y=std::log(y);
- }
- else{
- float a=3.75/b;
- //Bessel function evaluation
- //y=(exp(b)/sqrt(b))*(0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377))))))));
- //Logarithm of Bessel function
- y=b+std::log((0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377))))))))/std::sqrt(b));
- }
-
- return y;
- }
-
-
-
/////////////////////////////////////////////////////////////////////////////////////////////
////Class that represents one anisotropic compartment of the PVM model in an MCMC framework
////////////////////////////////////////////////////////////////////////////////////////////
@@ -163,6 +137,8 @@ namespace FIBRE{
m_f_prior=0;
compute_f_prior();
+ //cc OUT(m_f_prior);
+ //cc OUT(m_ardfudge);
m_lam_prior=0;
compute_lam_prior();
@@ -302,22 +278,27 @@ namespace FIBRE{
return true;
else{
//m_f_prior=-(log(m_lam) + (m_lam-1)*log(1-m_f));
- if(!can_use_ard)
+ if(!can_use_ard ){
m_f_prior=0;
+ }
else{
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;
}
}
@@ -456,6 +437,7 @@ namespace FIBRE{
inline bool propose_f( bool can_use_ard=true){
+ //cout<<"f"<<endl;
m_f_old=m_f;
m_f+=normrnd().AsScalar()*m_f_prop;
bool rejflag=compute_f_prior(can_use_ard);
@@ -478,6 +460,7 @@ namespace FIBRE{
inline bool propose_lam(){
+ //cout<<"lam"<<endl;
if(m_lam_jump){
m_lam_old=m_lam;
m_lam+=normrnd().AsScalar()*m_lam_prop;
@@ -566,7 +549,6 @@ namespace FIBRE{
float m_d_old_prior;
float m_d_acc;
float m_d_rej;
-
float m_d_std; //second d param in multi-shell model - std in gamma dist of diffusivities - to get non-monoexponential decay.
float m_d_std_old;
float m_d_std_prop;
@@ -574,7 +556,6 @@ namespace FIBRE{
float m_d_std_old_prior;
float m_d_std_acc;
float m_d_std_rej;
-
float m_S0;
float m_S0_old;
float m_S0_prop;
@@ -582,23 +563,6 @@ namespace FIBRE{
float m_S0_old_prior;
float m_S0_acc;
float m_S0_rej;
-
- float m_tau; //Precision parameter (1/sigma^2) for Rician noise
- float m_tau_old;
- float m_tau_prop;
- float m_tau_prior;
- float m_tau_old_prior;
- float m_tau_acc;
- float m_tau_rej;
-
- float m_f0; //Volume fraction for the unattenuated signal compartment
- float m_f0_old;
- float m_f0_prop;
- float m_f0_prior;
- float m_f0_old_prior;
- float m_f0_acc;
- float m_f0_rej;
-
float m_prior_en; //Joint Prior
float m_old_prior_en;
float m_likelihood_en; //Likelihood
@@ -613,29 +577,20 @@ namespace FIBRE{
const ColumnVector& m_alpha; //Theta angles of bvecs
const ColumnVector& m_beta; //Phi angles of bvecs
const Matrix& m_bvals; //b values
- const int m_modelnum; //1 for single-shell, 2 for multi-shell model
- const bool& m_rician; //If true, use Rician noise model
- const bool& m_includef0; //If true, include an unattenuated signal compartment in the model with fraction f0
- const bool& m_ardf0; //If true, use ard on the f0 compartment
+ const int& m_modelnum; //1 for single-shell, 2 for multi-shell model
+
public:
- //Constructor
+ // empty constructor
+
Multifibre(const ColumnVector& data,const ColumnVector& alpha,
- const ColumnVector& beta, const Matrix& b, int N ,float ardfudge=1, int modelnum=1, bool rician=false, bool inclf0=false, bool ardf0=false):
- m_ardfudge(ardfudge),m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b), m_modelnum(modelnum), m_rician(rician), m_includef0(inclf0), m_ardf0(ardf0){
+ const ColumnVector& beta, const Matrix& b, int N ,float ardfudge=1, int modelnum=1):
+ m_ardfudge(ardfudge),m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b), m_modelnum(modelnum){
+ // initialise(Amat,N);
m_iso_Signal.ReSize(alpha.Nrows());
m_iso_Signal=0;
m_iso_Signal_old=m_iso_Signal; //Initialize vectors that keep the signal from the isotropic compartment
-
- /* m_d_acc=0; m_d_rej=0;
- m_d_std_acc=0; m_d_std_rej=0;
- m_S0_acc=0; m_S0_rej=0;
- m_tau_acc=0; m_tau_rej=0;
- m_f0_acc=0; m_f0_rej=0;
- m_d_prior=0; m_d_std_prior=0; m_S0_prior=0; m_f0_prior=0; m_tau_prior=0; */
-
- m_f0=0.0; //Set this parameter to 0, it will be initialized later only if m_includef0 is true
- }
+ }
~Multifibre(){}
@@ -656,16 +611,11 @@ namespace FIBRE{
Fibre fib(m_alpha,m_beta,m_bvals,m_d,m_ardfudge,m_modelnum,m_d_std);
m_fibres.push_back(fib);
}
-
+
void initialise_energies(){
compute_d_prior();
if(m_modelnum==2)
compute_d_std_prior();
- if (m_rician){
- compute_tau_prior();
- }
- if (m_includef0)
- compute_f0_prior();
compute_S0_prior();
m_prior_en=0;
compute_prior();
@@ -678,11 +628,9 @@ namespace FIBRE{
//Initialize standard deviations for the proposal distributions
void initialise_props(){
- m_S0_prop=m_S0/10.0; //must have set initial values before this is called;
+ m_S0_prop=m_S0/10.0; //must have set inital values before this is called;
m_d_prop=m_d/10.0;
m_d_std_prop=m_d_std/10.0;
- m_tau_prop=m_tau/2.0;
- m_f0_prop=0.2;
}
@@ -701,12 +649,6 @@ namespace FIBRE{
inline float get_S0() const{ return m_S0;}
inline void set_S0(const float S0){ m_S0=S0; }
- inline float get_tau() const { return m_tau; }
- inline void set_tau(const float tau){ m_tau=tau;}
-
- inline float get_f0() const{ return m_f0;}
- inline void set_f0(const float f0){ m_f0=f0; }
-
inline void report() const{
OUT(m_d);
OUT(m_d_old);
@@ -749,20 +691,6 @@ namespace FIBRE{
}
}
-
- //More restrictive prior for the model with f0
- //particularly useful for the CSF voxels
- inline bool compute_d_prior_f0(){
- m_d_old_prior=m_d_prior;
- if(m_d<0 || m_d>0.008)
- return true;
- else{
- m_d_prior=0;
- return false;
- }
- }
-
-
inline bool compute_d_std_prior(){
m_d_std_old_prior=m_d_std_prior;
if(m_d_std<0)
@@ -784,34 +712,9 @@ namespace FIBRE{
}
- inline bool compute_tau_prior(){
- m_tau_old_prior=m_tau_prior;
- if(m_tau<=0)
- return true;
- else{
- m_tau_prior=0;
- return false;
- }
- }
-
-
- inline bool compute_f0_prior(){
- m_f0_old_prior=m_f0_prior;
- if (m_f0<=0 || m_f0>=1)
- return true;
- else{
- if(!m_ardf0) //Without ARD
- m_f0_prior=0;
- else //With ARD
- m_f0_prior=std::log(m_f0);
- return false;
- }
- }
-
-
//Check if sum of volume fractions is >1
bool reject_f_sum(){
- float fsum=m_f0;
+ float fsum=0;
for(unsigned int f=0;f<m_fibres.size();f++){
fsum+=m_fibres[f].get_f();
}
@@ -825,10 +728,6 @@ namespace FIBRE{
m_prior_en=m_d_prior+m_S0_prior;
if(m_modelnum==2)
m_prior_en=m_prior_en+m_d_std_prior;
- if(m_rician)
- m_prior_en=m_prior_en+m_tau_prior;
- if (m_includef0)
- m_prior_en=m_prior_en+m_f0_prior;
for(unsigned int f=0;f<m_fibres.size(); f++){
m_prior_en=m_prior_en+m_fibres[f].get_prior();
}
@@ -853,37 +752,29 @@ namespace FIBRE{
void compute_likelihood(){
m_old_likelihood_en=m_likelihood_en;
- ColumnVector pred(m_alpha.Nrows()),pred2;
- pred=0; pred2=pred;
- float fsum=m_f0;
+ ColumnVector pred(m_alpha.Nrows());
+ pred=0;
+ float fsum=0;
for(unsigned int f=0;f<m_fibres.size();f++){ //Signal from the anisotropic compartments
pred=pred+m_fibres[f].get_f()*m_fibres[f].getSignal();
fsum+=m_fibres[f].get_f(); //Total anisotropic volume fraction
}
for(int i=1;i<=pred.Nrows();i++){
- pred(i)=m_S0*(pred(i)+(1-fsum)*m_iso_Signal(i)+m_f0); //Add the signal from the isotropic compartment
- } //and multiply by S0 to get the total signal
+ pred(i)=m_S0*(pred(i)+(1-fsum)*m_iso_Signal(i)); //Add the signal from the isotropic compartment
+ } //and multiply by S0 to get the total signal
- if (!m_rician){ //Likelihood energy for Gaussian noise
- float sumsquares=(m_data-pred).SumSquare(); //Sum of squared residuals
- m_likelihood_en=(m_data.Nrows()/2.0)*log(sumsquares/2.0);
- }
- else{ //Likelihood energy for Rician noise
- for(int i=1;i<=pred.Nrows();i++)
- pred2(i)=std::log(m_data(i))-0.5*m_tau*(m_data(i)*m_data(i)+pred(i)*pred(i))+logIo(m_tau*pred(i)*m_data(i));
- m_likelihood_en=-m_data.Nrows()*std::log(m_tau)-pred2.Sum();
- }
- }
+ float sumsquares=(m_data-pred).SumSquare(); //Sum of squared residuals
+ m_likelihood_en=(m_data.Nrows()/2.0)*log(sumsquares/2.0);
+ }
- inline void compute_energy(){
+ inline void compute_energy(){
m_old_energy=m_energy;
m_energy=m_prior_en+m_likelihood_en;
- //cout<<m_prior_en<<" "<<m_likelihood_en<<endl;
}
- /*
+
void evaluate_energy(){
m_old_energy=m_energy;
for(unsigned int f=0;f<m_fibres.size();f++){
@@ -897,12 +788,13 @@ namespace FIBRE{
compute_prior();
compute_likelihood();
m_energy=m_prior_en+m_likelihood_en;
- } */
+ }
//Test whether the candidate state will be accepted
bool test_energy(){
float tmp=exp(m_old_energy-m_energy);
+ //cout<<m_old_energy<< " " << m_energy <<endl;
return (tmp>unifrnd().AsScalar());
}
@@ -922,13 +814,9 @@ namespace FIBRE{
inline bool propose_d(){
//cout<<"d"<<endl;
- bool rejflag;
m_d_old=m_d;
m_d+=normrnd().AsScalar()*m_d_prop;
- if (m_includef0)
- rejflag=compute_d_prior_f0();
- else
- rejflag=compute_d_prior();
+ bool rejflag=compute_d_prior();//inside this it stores the old prior
for(unsigned int f=0;f<m_fibres.size();f++)
m_fibres[f].compute_signal();
compute_iso_signal();
@@ -953,6 +841,7 @@ namespace FIBRE{
inline bool propose_d_std(){
+ //cout<<"d_std"<<endl;
m_d_std_old=m_d_std;
m_d_std+=normrnd().AsScalar()*m_d_std_prop;
bool rejflag=compute_d_std_prior();//inside this it stores the old prior
@@ -980,6 +869,7 @@ namespace FIBRE{
inline bool propose_S0(){
+ //cout<<"S0"<<endl;
m_S0_old=m_S0;
m_S0+=normrnd().AsScalar()*m_S0_prop;
bool rejflag=compute_S0_prior();//inside this it stores the old prior
@@ -1000,47 +890,6 @@ namespace FIBRE{
}
- inline bool propose_tau(){
- m_tau_old=m_tau;
- m_tau+=normrnd().AsScalar()*m_tau_prop;
- bool rejflag=compute_tau_prior();//inside this it stores the old prior
- return rejflag;
- };
-
-
- inline void accept_tau(){
- m_tau_acc++;
- }
-
- inline void reject_tau(){
- m_tau=m_tau_old;
- m_tau_prior=m_tau_old_prior;
- m_prior_en=m_old_prior_en;
- m_tau_rej++;
- }
-
-
- inline bool propose_f0(){
- m_f0_old=m_f0;
- m_f0+=normrnd().AsScalar()*m_f0_prop;
- bool rejflag=compute_f0_prior();
- compute_prior();
- return rejflag;
- };
-
-
- inline void accept_f0(){
- m_f0_acc++;
- }
-
- inline void reject_f0(){
- m_f0=m_f0_old;
- m_f0_prior=m_f0_old_prior;
- m_prior_en=m_old_prior_en;
- m_f0_rej++;
- }
-
-
//Adapt standard deviation of proposal distributions during MCMC execution
//to avoid over-rejection/over-acceptance of MCMC samples
inline void update_proposals(){
@@ -1052,18 +901,6 @@ namespace FIBRE{
m_d_std_acc=0;
m_d_std_rej=0;
}
- if (m_rician){
- m_tau_prop*=sqrt(float(m_tau_acc+1)/float(m_tau_rej+1));
- m_tau_prop=min(m_tau_prop,maxfloat);
- m_tau_acc=0;
- m_tau_rej=0;
- }
- if (m_includef0){
- m_f0_prop*=sqrt(float(m_f0_acc+1)/float(m_f0_rej+1));
- m_f0_prop=min(m_f0_prop,maxfloat);
- m_f0_acc=0;
- m_f0_rej=0;
- }
m_S0_prop*=sqrt(float(m_S0_acc+1)/float(m_S0_rej+1));
m_S0_prop=min(m_S0_prop,maxfloat);
m_d_acc=0;
@@ -1079,43 +916,9 @@ namespace FIBRE{
//Function that performs a single MCMC iteration
void jump( bool can_use_ard=true ){
- if (m_includef0){ //Try f0
- if(!propose_f0()){
- if(!reject_f_sum()){
- compute_prior();
- compute_likelihood();
- compute_energy();
- if(test_energy())
- accept_f0();
- else{
- reject_f0();
- restore_energy();
- }
- }
- else //else for rejecting fsum>1
- reject_f0();
- }
- else //else for rejecting rejflag returned from propose_f()
- reject_f0();
- }
-
-
- if(m_rician){ //Try tau
- if(!propose_tau()){
- compute_prior();
- compute_likelihood();
- compute_energy();
- if(test_energy())
- accept_tau();
- else{
- reject_tau();
- restore_energy();
- }
- }
- else
- reject_tau();
- }
-
+ // cout<<"modelnum "<<m_modelnum<<endl;
+ // cout<<"fmodelnum "<<m_fibres[0].get_modelnum()<<endl;
+ // cout<<"fmodelnum2 "<<m_fibres[1].get_modelnum()<<endl;
if(!propose_d()){ //Try d
compute_prior();
compute_likelihood();
@@ -1128,9 +931,9 @@ namespace FIBRE{
restore_energy();
}
}
- else
+ else{
reject_d();
-
+ }
if(m_modelnum==2){ //Try d_std
if(!propose_d_std()){
@@ -1145,11 +948,10 @@ namespace FIBRE{
restore_energy();
}
}
- else
+ else{
reject_d_std();
+ }
}
-
-
if(!propose_S0()){ //Try S0
compute_prior();
compute_likelihood();
@@ -1161,29 +963,44 @@ namespace FIBRE{
restore_energy();
}
}
- else
+ else{
reject_S0();
+ }
-
-
-
-
for(unsigned int f=0;f<m_fibres.size();f++){ //For each fibre
+ //cout<<"pre th "<<m_fibres[f].get_th()<<endl;
+ //cout << m_fibres[f].getSignal()<<endl;
+ //cout<<"prior"<<m_prior_en<<endl;
+ //cout<<"lik"<<m_likelihood_en<<endl;
+ //cout<<"energy"<<m_energy<<endl;
+ //cout <<endl;
if(!m_fibres[f].propose_th()){ //Try theta
+
compute_prior();
compute_likelihood();
compute_energy();
- if(test_energy())
+ //cout<<"post th "<<m_fibres[f].get_th()<<endl;
+ //cout<< m_d<< " " <<m_d_std<<endl;
+ //cout << m_fibres[f].getSignal()<<endl;
+ //cout<<"prior"<<m_prior_en<<endl;
+ //cout<<"lik"<<m_likelihood_en<<endl;
+ //cout<<"energy"<<m_energy<<endl;
+
+ if(test_energy()){
m_fibres[f].accept_th();
+ //cout<<"Accepted"<<endl;}
+ }
else{
m_fibres[f].reject_th();
restore_energy();
+ //cout<<"rejected"<<endl;
}
+ //cout <<endl<<endl;
}
- else
+ else {
m_fibres[f].reject_th();
-
+ }
if(!m_fibres[f].propose_ph()){ //Try phi
compute_prior();
@@ -1196,28 +1013,30 @@ namespace FIBRE{
restore_energy();
}
}
- else
- m_fibres[f].reject_ph();
-
+ else{
+ m_fibres[f].reject_ph();
+ }
if(!m_fibres[f].propose_f( can_use_ard )){ //Try f
if(!reject_f_sum()){
compute_prior();
compute_likelihood();
compute_energy();
- if(test_energy())
+ if(test_energy()){
m_fibres[f].accept_f();
+ }
else{
m_fibres[f].reject_f();
restore_energy();
}
}
- else //else for rejectin fsum>1
+ else{//else for rejectin fsum>1
m_fibres[f].reject_f();
- }
- else //else for rejecting rejflag returned from propose_f()
- m_fibres[f].reject_f();
-
+ }
+ }
+ else{//else for rejecting rejflag returned from propose_f()
+ m_fibres[f].reject_f();
+ }
// if(!m_fibres[f].propose_lam()){
// compute_prior();
@@ -1236,7 +1055,7 @@ namespace FIBRE{
}
}
-
+
Multifibre& operator=(const Multifibre& rhs){
m_fibres=rhs.m_fibres;
m_d=rhs.m_d;
@@ -1260,20 +1079,6 @@ namespace FIBRE{
m_S0_old_prior=rhs.m_S0_old_prior;
m_S0_acc=rhs.m_S0_acc;
m_S0_rej=rhs.m_S0_rej;
- m_f0=rhs.m_f0;
- m_f0_old=rhs.m_f0_old;
- m_f0_prop=rhs.m_f0_prop;
- m_f0_prior=rhs.m_f0_prior;
- m_f0_old_prior=rhs.m_f0_old_prior;
- m_f0_acc=rhs.m_f0_acc;
- m_f0_rej=rhs.m_f0_rej;
- m_tau=rhs.m_tau;
- m_tau_old=rhs.m_tau_old;
- m_tau_prop=rhs.m_tau_prop;
- m_tau_prior=rhs.m_tau_prior;
- m_tau_old_prior=rhs.m_tau_old_prior;
- m_tau_acc=rhs.m_tau_acc;
- m_tau_rej=rhs.m_tau_rej;
m_prior_en=rhs.m_prior_en;
m_old_prior_en=rhs.m_old_prior_en;
m_likelihood_en=rhs.m_likelihood_en;
@@ -1284,7 +1089,7 @@ namespace FIBRE{
m_iso_Signal=rhs.m_iso_Signal;
m_iso_Signal_old=rhs.m_iso_Signal_old;
return *this;
- }
+ }
};