Skip to content
Snippets Groups Projects
Commit c08f6e83 authored by Stamatios Sotiropoulos's avatar Stamatios Sotiropoulos
Browse files

Added function to compute signal from isotropic compartment. Trigonmetric... 

Added function to compute signal from isotropic compartment. Trigonmetric identities used to speed up likelihood calc
parent 9cc09df6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 150 additions and 75 deletions
fibre.h
+ 150
− 75
View file @ c08f6e83
...@@ -23,20 +23,23 @@ const float minlogfloat=-23; ...@@ -23,20 +23,23 @@ const float minlogfloat=-23;
namespace FIBRE{ namespace FIBRE{
class Fibre{ /////////////////////////////////////////////////////////////////////////////////////////////
float m_th; ////Class that represents one anisotropic compartment of the PVM model in an MCMC framework
////////////////////////////////////////////////////////////////////////////////////////////
class Fibre{
float m_th; //Current/candidate MCMC state
float m_ph; float m_ph;
float m_f; float m_f;
float m_lam; float m_lam;
float m_th_prop; float m_th_prop; //Standard deviation for Gaussian proposal distributions of parameters
float m_ph_prop; float m_ph_prop;
float m_f_prop; float m_f_prop;
float m_lam_prop; float m_lam_prop;
float m_th_old; float m_th_old; //Last accepted value. If a sample is rejected this value is restored
float m_ph_old; float m_ph_old;
float m_f_old; float m_f_old;
float m_lam_old; float m_lam_old;
float m_th_prior; float m_th_prior; //Priors for the model parameters
float m_ph_prior; float m_ph_prior;
float m_f_prior; float m_f_prior;
float m_lam_prior; float m_lam_prior;
...@@ -44,10 +47,10 @@ namespace FIBRE{ ...@@ -44,10 +47,10 @@ namespace FIBRE{
float m_ph_old_prior; float m_ph_old_prior;
float m_f_old_prior; float m_f_old_prior;
float m_lam_old_prior; float m_lam_old_prior;
float m_prior_en; float m_prior_en; //Joint Prior
float m_old_prior_en; float m_old_prior_en;
float m_ardfudge; float m_ardfudge;
int m_th_acc; int m_th_acc;
int m_th_rej; int m_th_rej;
int m_ph_acc; int m_ph_acc;
int m_ph_rej; int m_ph_rej;
...@@ -56,14 +59,14 @@ namespace FIBRE{ ...@@ -56,14 +59,14 @@ namespace FIBRE{
int m_lam_acc; int m_lam_acc;
int m_lam_rej; int m_lam_rej;
bool m_lam_jump; bool m_lam_jump;
ColumnVector m_Signal; ColumnVector m_Signal; //Vector that stores the signal from the anisotropic compartment during the candidate/current MCMC state
ColumnVector m_Signal_old; ColumnVector m_Signal_old;
const float& m_d; const float& m_d;
const float& m_d_std; //second d param in multi-shell model - std in gamma dist of diffusivities const float& m_d_std; //second d param in multi-shell model - std in gamma dist of diffusivities
const ColumnVector& m_alpha; const ColumnVector& m_alpha; //Theta angles of bvecs
const ColumnVector& m_beta; const ColumnVector& m_beta; //Phi angles of bvecs
const Matrix& m_bvals; const Matrix& m_bvals; //b values
const int& m_modelnum; const int& m_modelnum; //1 for single-shell, 2 for multi-shell model
public: public:
//constructors:: //constructors::
...@@ -228,15 +231,17 @@ namespace FIBRE{ ...@@ -228,15 +231,17 @@ namespace FIBRE{
inline float get_modelnum() const{ return m_modelnum;} inline float get_modelnum() const{ return m_modelnum;}
inline void update_proposals(){ //Adapt the standard deviation of the proposal distributions during MCMC execution
//to avoid over-rejection/over-acceptance of samples
inline void update_proposals(){
m_th_prop*=sqrt(float(m_th_acc+1)/float(m_th_rej+1)); m_th_prop*=sqrt(float(m_th_acc+1)/float(m_th_rej+1));
m_th_prop=min(m_th_prop,maxfloat); m_th_prop=min(m_th_prop,maxfloat);
m_ph_prop*=sqrt(float(m_ph_acc+1)/float(m_ph_rej+1)); m_ph_prop*=sqrt(float(m_ph_acc+1)/float(m_ph_rej+1));
m_ph_prop=min(m_ph_prop,maxfloat); m_ph_prop=min(m_ph_prop,maxfloat);
m_f_prop*=sqrt(float(m_f_acc+1)/float(m_f_rej+1)); m_f_prop*=sqrt(float(m_f_acc+1)/float(m_f_rej+1));
m_f_prop=min(m_f_prop,maxfloat); m_f_prop=min(m_f_prop,maxfloat);
m_lam_prop*=sqrt(float(m_lam_acc+1)/float(m_lam_rej+1)); //m_lam_prop*=sqrt(float(m_lam_acc+1)/float(m_lam_rej+1)); //STAM: REmoved Update for m_lam_prior
m_lam_prop=min(m_lam_prop,maxfloat); //m_lam_prop=min(m_lam_prop,maxfloat);
m_th_acc=0; m_th_acc=0;
m_th_rej=0; m_th_rej=0;
m_ph_acc=0; m_ph_acc=0;
...@@ -256,11 +261,15 @@ namespace FIBRE{ ...@@ -256,11 +261,15 @@ namespace FIBRE{
} }
return false; //instant rejection flag return false; //instant rejection flag
} }
inline bool compute_ph_prior(){ inline bool compute_ph_prior(){
m_ph_old_prior=m_ph_prior; m_ph_old_prior=m_ph_prior;
m_ph_prior=0; m_ph_prior=0;
return false; return false;
} }
inline bool compute_f_prior(bool can_use_ard=true){ inline bool compute_f_prior(bool can_use_ard=true){
//note(gamma(lam+1)/(gamma(1)*gamma(lam)) = lam //note(gamma(lam+1)/(gamma(1)*gamma(lam)) = lam
// the following is a beta distribution with alpha=0 // the following is a beta distribution with alpha=0
...@@ -273,7 +282,7 @@ namespace FIBRE{ ...@@ -273,7 +282,7 @@ namespace FIBRE{
m_f_prior=0; m_f_prior=0;
} }
else{ else{
if(m_lam_jump){ //if(m_lam_jump){ //STAM: Removed the if-else statement (always true since m_lam_jump=true)
// m_f_prior=log(1-m_f)+2*log(fabs(log(1-m_f))); //marginalised with uniform prior on lambda // 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); m_f_prior=std::log(m_f);
...@@ -283,9 +292,9 @@ namespace FIBRE{ ...@@ -283,9 +292,9 @@ namespace FIBRE{
//cc OUT(mmk); //cc OUT(mmk);
//cc OUT(m_f_old_prior); //cc OUT(m_f_old_prior);
} //}
else //else
m_f_prior=0; //m_f_prior=0;
} }
m_f_prior=m_ardfudge*m_f_prior; m_f_prior=m_ardfudge*m_f_prior;
...@@ -293,7 +302,8 @@ namespace FIBRE{ ...@@ -293,7 +302,8 @@ namespace FIBRE{
return false; return false;
} }
} }
inline bool compute_lam_prior(){ inline bool compute_lam_prior(){
m_lam_old_prior=m_lam_prior; m_lam_old_prior=m_lam_prior;
if(m_lam <0 | m_lam > 1e16 ) if(m_lam <0 | m_lam > 1e16 )
...@@ -304,23 +314,31 @@ namespace FIBRE{ ...@@ -304,23 +314,31 @@ namespace FIBRE{
} }
} }
inline void compute_prior(){ inline void compute_prior(){
m_old_prior_en=m_prior_en; m_old_prior_en=m_prior_en;
m_prior_en=m_th_prior+m_ph_prior+m_f_prior+m_lam_prior; //m_prior_en=m_th_prior+m_ph_prior+m_f_prior+m_lam_prior; //STAM:Removed m_lam_prior
m_prior_en=m_th_prior+m_ph_prior+m_f_prior;
} }
void compute_signal(){
//Compute model predicted signal, only due to the anisotropic compartment
void compute_signal(){
m_Signal_old=m_Signal; m_Signal_old=m_Signal;
for (int i = 1; i <= m_alpha.Nrows(); i++){ for (int i = 1; i <= m_alpha.Nrows(); i++){
float angtmp=cos(m_ph-m_beta(i))*sin(m_alpha(i))*sin(m_th) + cos(m_alpha(i))*cos(m_th); float cos_alpha_minus_theta=cos(m_alpha(i)-m_th); //STAM: Used trigonometric identity to reduce the number of sinusoids evaluated
float cos_alpha_plus_theta=cos(m_alpha(i)+m_th);
//float angtmp=cos(m_ph-m_beta(i))*sin(m_alpha(i))*sin(m_th) + cos(m_alpha(i))*cos(m_th);
float angtmp=cos(m_ph-m_beta(i))*(cos_alpha_minus_theta-cos_alpha_plus_theta)/2 + (cos_alpha_minus_theta+cos_alpha_plus_theta)/2;
angtmp=angtmp*angtmp; angtmp=angtmp*angtmp;
if(m_modelnum==1){ if(m_modelnum==1){
m_Signal(i)=exp(-m_d*m_bvals(1,i)*angtmp); m_Signal(i)=exp(-m_d*m_bvals(1,i)*angtmp);
// cout<<"here"<<endl; // cout<<"here"<<endl;
} }
else if(m_modelnum==2){ else if(m_modelnum==2){
float dalpha=m_d*m_d/(m_d_std*m_d_std);
float dbeta=m_d/(m_d_std*m_d_std); float dbeta=m_d/(m_d_std*m_d_std);
float dalpha=m_d*dbeta; //STAM: Removed one division
m_Signal(i)=exp( log(dbeta/(dbeta + m_bvals(1,i)*angtmp)) *dalpha ); //gamma distribution of diffusivities - params alpha (m_d) and beta (m_d_beta): Expected signal is (beta./(beta+b*g(th,ph))).^alpha; m_Signal(i)=exp( log(dbeta/(dbeta + m_bvals(1,i)*angtmp)) *dalpha ); //gamma distribution of diffusivities - params alpha (m_d) and beta (m_d_beta): Expected signal is (beta./(beta+b*g(th,ph))).^alpha;
//cout<<i<<" "<<angtmp<<" "<<m_d_std<<" "<<m_d<<endl; //cout<<i<<" "<<angtmp<<" "<<m_d_std<<" "<<m_d<<endl;
//cout<<"balh "<<m_Signal(i)<<endl; //cout<<"balh "<<m_Signal(i)<<endl;
...@@ -341,10 +359,12 @@ namespace FIBRE{ ...@@ -341,10 +359,12 @@ namespace FIBRE{
return rejflag; return rejflag;
}; };
inline void accept_th(){ inline void accept_th(){
m_th_acc++; m_th_acc++;
} }
inline void reject_th(){ inline void reject_th(){
m_th=m_th_old; m_th=m_th_old;
m_th_prior=m_th_old_prior; m_th_prior=m_th_old_prior;
...@@ -353,6 +373,7 @@ namespace FIBRE{ ...@@ -353,6 +373,7 @@ namespace FIBRE{
m_th_rej++; m_th_rej++;
} }
inline bool propose_ph(){ inline bool propose_ph(){
//cout<<"ph"<<endl; //cout<<"ph"<<endl;
m_ph_old=m_ph; m_ph_old=m_ph;
...@@ -363,10 +384,12 @@ namespace FIBRE{ ...@@ -363,10 +384,12 @@ namespace FIBRE{
return rejflag; return rejflag;
}; };
inline void accept_ph(){ inline void accept_ph(){
m_ph_acc++; m_ph_acc++;
} }
inline void reject_ph(){ inline void reject_ph(){
m_ph=m_ph_old; m_ph=m_ph_old;
m_ph_prior=m_ph_old_prior; m_ph_prior=m_ph_old_prior;
...@@ -388,10 +411,14 @@ namespace FIBRE{ ...@@ -388,10 +411,14 @@ namespace FIBRE{
return rejflag_th | rejflag_ph; return rejflag_th | rejflag_ph;
} }
void accept_th_ph(){ void accept_th_ph(){
m_th_acc++; m_th_acc++;
m_ph_acc++; m_ph_acc++;
} }
void reject_th_ph(){ void reject_th_ph(){
m_ph=m_ph_old; m_ph=m_ph_old;
m_ph_prior=m_ph_old_prior; m_ph_prior=m_ph_old_prior;
...@@ -406,6 +433,7 @@ namespace FIBRE{ ...@@ -406,6 +433,7 @@ namespace FIBRE{
m_ph_rej++; m_ph_rej++;
} }
inline bool propose_f( bool can_use_ard=true){ inline bool propose_f( bool can_use_ard=true){
//cout<<"f"<<endl; //cout<<"f"<<endl;
m_f_old=m_f; m_f_old=m_f;
...@@ -415,9 +443,11 @@ namespace FIBRE{ ...@@ -415,9 +443,11 @@ namespace FIBRE{
return rejflag; return rejflag;
}; };
inline void accept_f(){ inline void accept_f(){
m_f_acc++; m_f_acc++;
} }
inline void reject_f(){ inline void reject_f(){
m_f=m_f_old; m_f=m_f_old;
...@@ -425,6 +455,7 @@ namespace FIBRE{ ...@@ -425,6 +455,7 @@ namespace FIBRE{
m_prior_en=m_old_prior_en; m_prior_en=m_old_prior_en;
m_f_rej++; m_f_rej++;
} }
inline bool propose_lam(){ inline bool propose_lam(){
//cout<<"lam"<<endl; //cout<<"lam"<<endl;
...@@ -438,10 +469,12 @@ namespace FIBRE{ ...@@ -438,10 +469,12 @@ namespace FIBRE{
} }
else {return true;} else {return true;}
}; };
inline void accept_lam(){ inline void accept_lam(){
m_lam_acc++; m_lam_acc++;
} }
inline void reject_lam(){ inline void reject_lam(){
m_lam=m_lam_old; m_lam=m_lam_old;
...@@ -492,6 +525,7 @@ namespace FIBRE{ ...@@ -492,6 +525,7 @@ namespace FIBRE{
friend ostream& operator<<(ostream& ostr,const Fibre& p); friend ostream& operator<<(ostream& ostr,const Fibre& p);
}; };
//overload << //overload <<
inline ostream& operator<<(ostream& ostr,const Fibre& p){ inline ostream& operator<<(ostream& ostr,const Fibre& p){
ostr<<p.m_th<<" "<<p.m_ph<<" "<<p.m_f<<endl; ostr<<p.m_th<<" "<<p.m_ph<<" "<<p.m_f<<endl;
...@@ -499,6 +533,11 @@ namespace FIBRE{ ...@@ -499,6 +533,11 @@ namespace FIBRE{
} }
/////////////////////////////////////////////////////////////////////////////////////////////////////
//Class that represents the PVM model in an MCMC framework. An isotropic compartment and M>=1
//anisotropic compartments (Fibre objects) are considered. Functions are defined for performing
//a single MCMC iteration
/////////////////////////////////////////////////////////////////////////////////////////////////////
class Multifibre{ class Multifibre{
vector<Fibre> m_fibres; vector<Fibre> m_fibres;
float m_d; float m_d;
...@@ -522,19 +561,21 @@ namespace FIBRE{ ...@@ -522,19 +561,21 @@ namespace FIBRE{
float m_S0_old_prior; float m_S0_old_prior;
float m_S0_acc; float m_S0_acc;
float m_S0_rej; float m_S0_rej;
float m_prior_en; float m_prior_en; //Joint Prior
float m_old_prior_en; float m_old_prior_en;
float m_likelihood_en; float m_likelihood_en; //Likelihood
float m_old_likelihood_en; float m_old_likelihood_en;
float m_energy; float m_energy; //Posterior
float m_old_energy; float m_old_energy;
float m_ardfudge; float m_ardfudge;
ColumnVector m_iso_Signal; //STAM:Vector that stores the signal from the isotropic compartment during the candidate/current state
ColumnVector m_iso_Signal_old;
const ColumnVector& m_data; const ColumnVector& m_data; //Data vector
const ColumnVector& m_alpha; const ColumnVector& m_alpha; //Theta angles of bvecs
const ColumnVector& m_beta; const ColumnVector& m_beta; //Phi angles of bvecs
const Matrix& m_bvals; const Matrix& m_bvals; //b values
const int& m_modelnum; const int& m_modelnum; //1 for single-shell, 2 for multi-shell model
public: public:
...@@ -544,7 +585,9 @@ namespace FIBRE{ ...@@ -544,7 +585,9 @@ namespace FIBRE{
const ColumnVector& beta, const Matrix& b, int N ,float ardfudge=1, int modelnum=1): 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){ m_ardfudge(ardfudge),m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b), m_modelnum(modelnum){
// initialise(Amat,N); // initialise(Amat,N);
m_iso_Signal.ReSize(alpha.Nrows());
m_iso_Signal=0;
m_iso_Signal_old=m_iso_Signal; //STAM: Initialize vectors that keep the signal from the isotropic compartment
} }
~Multifibre(){} ~Multifibre(){}
...@@ -552,13 +595,14 @@ namespace FIBRE{ ...@@ -552,13 +595,14 @@ namespace FIBRE{
const vector<Fibre>& fibres() const{ const vector<Fibre>& fibres() const{
return m_fibres; return m_fibres;
} }
vector<Fibre>& get_fibres(){ vector<Fibre>& get_fibres(){
return m_fibres; return m_fibres;
} }
void addfibre(const float th, const float ph, const float f,const float lam, const bool lam_jump=true){ void addfibre(const float th, const float ph, const float f,const float lam, const bool lam_jump=true){
Fibre fib(m_alpha,m_beta,m_bvals,m_d,m_ardfudge,th,ph,f,lam,lam_jump,m_modelnum,m_d_std); Fibre fib(m_alpha,m_beta,m_bvals,m_d,m_ardfudge,th,ph,f,lam,lam_jump,m_modelnum,m_d_std);
m_fibres.push_back(fib); m_fibres.push_back(fib);
} }
void addfibre(){ void addfibre(){
...@@ -574,16 +618,20 @@ namespace FIBRE{ ...@@ -574,16 +618,20 @@ namespace FIBRE{
m_prior_en=0; m_prior_en=0;
compute_prior(); compute_prior();
m_likelihood_en=0; m_likelihood_en=0;
compute_iso_signal(); //STAM:Initialize the isotropic compartment signal before getting the likelihood for the first time
compute_likelihood(); compute_likelihood();
compute_energy(); compute_energy();
} }
void initialise_props(){
m_S0_prop=m_S0/10; //must have set inital values before this is called; //Initialize standard deviations for the proposal distributions
m_d_prop=m_d/10; void initialise_props(){
m_d_std_prop=m_d_std/10; 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;
} }
inline float get_d() const{ return m_d;} inline float get_d() const{ return m_d;}
inline void set_d(const float d){ m_d=d; } inline void set_d(const float d){ m_d=d; }
...@@ -661,6 +709,8 @@ namespace FIBRE{ ...@@ -661,6 +709,8 @@ namespace FIBRE{
} }
} }
//Check if sum of volume fractions is >1
bool reject_f_sum(){ bool reject_f_sum(){
float fsum=0; float fsum=0;
for(unsigned int f=0;f<m_fibres.size();f++){ for(unsigned int f=0;f<m_fibres.size();f++){
...@@ -669,6 +719,8 @@ namespace FIBRE{ ...@@ -669,6 +719,8 @@ namespace FIBRE{
return fsum>1; //true if sum(f) > 1 and therefore, we should reject f return fsum>1; //true if sum(f) > 1 and therefore, we should reject f
} }
//Compute Joint Prior
inline void compute_prior(){ inline void compute_prior(){
m_old_prior_en=m_prior_en; m_old_prior_en=m_prior_en;
m_prior_en=m_d_prior+m_S0_prior; m_prior_en=m_d_prior+m_S0_prior;
...@@ -679,7 +731,23 @@ namespace FIBRE{ ...@@ -679,7 +731,23 @@ namespace FIBRE{
} }
} }
void compute_likelihood(){
void compute_iso_signal(){ //STAM: Function that computes the signal from the isotropic compartment
m_iso_Signal_old=m_iso_Signal;
for(int i=1;i<=m_alpha.Nrows();i++)
{
if(m_modelnum==1)
m_iso_Signal(i)=exp(-m_d*m_bvals(1,i));
else if(m_modelnum==2)
{
float dbeta=m_d/(m_d_std*m_d_std);
float dalpha=m_d*dbeta;
m_iso_Signal(i)=exp(log(dbeta/(dbeta+m_bvals(1,i)))*dalpha);
}
}
}
void compute_likelihood(){
m_old_likelihood_en=m_likelihood_en; m_old_likelihood_en=m_likelihood_en;
ColumnVector pred(m_alpha.Nrows()); ColumnVector pred(m_alpha.Nrows());
...@@ -690,63 +758,59 @@ namespace FIBRE{ ...@@ -690,63 +758,59 @@ namespace FIBRE{
fsum+=m_fibres[f].get_f(); fsum+=m_fibres[f].get_f();
} }
for(int i=1;i<=pred.Nrows();i++){ for(int i=1;i<=pred.Nrows();i++)
if(m_modelnum==1) pred(i)=pred(i)+(1-fsum)*m_iso_Signal(i); //STAM: Use the precalculated isotropic compartment
pred(i)=pred(i)+(1-fsum)*exp(-m_d*m_bvals(1,i));
else if(m_modelnum==2){
float dalpha=m_d*m_d/(m_d_std*m_d_std);
float dbeta=m_d/(m_d_std*m_d_std);
pred(i)=pred(i)+(1-fsum)*exp(log(dbeta/(dbeta+m_bvals(1,i)))*dalpha);
}
}
pred=pred*m_S0; pred=pred*m_S0;
float sumsquares=(m_data-pred).SumSquare(); float sumsquares=(m_data-pred).SumSquare();
m_likelihood_en=(m_data.Nrows()/2)*log(sumsquares/2); 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_old_energy=m_energy;
m_energy=m_prior_en+m_likelihood_en; m_energy=m_prior_en+m_likelihood_en;
} }
void evaluate_energy(){ void evaluate_energy(){
m_old_energy=m_energy; m_old_energy=m_energy;
for(unsigned int f=0;f<m_fibres.size();f++){ for(unsigned int f=0;f<m_fibres.size();f++){
m_fibres[f].compute_th_prior(); m_fibres[f].compute_th_prior();
m_fibres[f].compute_ph_prior(); m_fibres[f].compute_ph_prior();
m_fibres[f].compute_f_prior(); m_fibres[f].compute_f_prior();
m_fibres[f].compute_lam_prior(); //m_fibres[f].compute_lam_prior(); //STAM: Removed the lam_prior calulation
m_fibres[f].compute_signal(); m_fibres[f].compute_signal();
m_fibres[f].compute_prior(); m_fibres[f].compute_prior();
} }
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
m_energy=m_prior_en+m_likelihood_en; m_energy=m_prior_en+m_likelihood_en;
} }
//Test whether the candidate state will be accepted
bool test_energy(){ bool test_energy(){
float tmp=exp(m_old_energy-m_energy); float tmp=exp(m_old_energy-m_energy);
//cout<<m_old_energy<< " " << m_energy <<endl; //cout<<m_old_energy<< " " << m_energy <<endl;
return (tmp>unifrnd().AsScalar()); return (tmp>unifrnd().AsScalar());
} }
inline void restore_energy(){ inline void restore_energy(){
m_prior_en=m_old_prior_en; m_prior_en=m_old_prior_en;
m_likelihood_en=m_old_likelihood_en; m_likelihood_en=m_old_likelihood_en;
m_energy=m_old_energy; m_energy=m_old_energy;
} }
inline void restore_energy_no_lik(){ inline void restore_energy_no_lik(){
m_prior_en=m_old_prior_en; m_prior_en=m_old_prior_en;
m_energy=m_old_energy; m_energy=m_old_energy;
} }
inline bool propose_d(){ inline bool propose_d(){
//cout<<"d"<<endl; //cout<<"d"<<endl;
m_d_old=m_d; m_d_old=m_d;
...@@ -754,12 +818,15 @@ namespace FIBRE{ ...@@ -754,12 +818,15 @@ namespace FIBRE{
bool rejflag=compute_d_prior();//inside this it stores the old prior bool rejflag=compute_d_prior();//inside this it stores the old prior
for(unsigned int f=0;f<m_fibres.size();f++) for(unsigned int f=0;f<m_fibres.size();f++)
m_fibres[f].compute_signal(); m_fibres[f].compute_signal();
compute_iso_signal(); //STAM: Compute the signal from the isotropic compartment
return rejflag; return rejflag;
}; };
inline void accept_d(){ inline void accept_d(){
m_d_acc++; m_d_acc++;
} }
inline void reject_d(){ inline void reject_d(){
m_d=m_d_old; m_d=m_d_old;
...@@ -767,8 +834,11 @@ namespace FIBRE{ ...@@ -767,8 +834,11 @@ namespace FIBRE{
m_prior_en=m_old_prior_en; m_prior_en=m_old_prior_en;
for(unsigned int f=0;f<m_fibres.size();f++) for(unsigned int f=0;f<m_fibres.size();f++)
m_fibres[f].restoreSignal(); m_fibres[f].restoreSignal();
m_iso_Signal=m_iso_Signal_old; //STAM:Restore the signal from the isotropic compartment
m_d_rej++; m_d_rej++;
} }
inline bool propose_d_std(){ inline bool propose_d_std(){
//cout<<"d_std"<<endl; //cout<<"d_std"<<endl;
m_d_std_old=m_d_std; m_d_std_old=m_d_std;
...@@ -776,22 +846,27 @@ namespace FIBRE{ ...@@ -776,22 +846,27 @@ namespace FIBRE{
bool rejflag=compute_d_std_prior();//inside this it stores the old prior bool rejflag=compute_d_std_prior();//inside this it stores the old prior
for(unsigned int f=0;f<m_fibres.size();f++) for(unsigned int f=0;f<m_fibres.size();f++)
m_fibres[f].compute_signal(); m_fibres[f].compute_signal();
compute_iso_signal(); //STAM: Compute the signal from the isotropic compartment
return rejflag; return rejflag;
}; };
inline void accept_d_std(){ inline void accept_d_std(){
m_d_std_acc++; m_d_std_acc++;
} }
inline void reject_d_std(){ inline void reject_d_std(){
m_d_std=m_d_std_old; m_d_std=m_d_std_old;
m_d_std_prior=m_d_std_old_prior; m_d_std_prior=m_d_std_old_prior;
m_prior_en=m_old_prior_en; m_prior_en=m_old_prior_en;
for(unsigned int f=0;f<m_fibres.size();f++) for(unsigned int f=0;f<m_fibres.size();f++)
m_fibres[f].restoreSignal(); m_fibres[f].restoreSignal();
m_iso_Signal=m_iso_Signal_old; //STAM:Restore the signal from the isotropic compartment
m_d_std_rej++; m_d_std_rej++;
} }
inline bool propose_S0(){ inline bool propose_S0(){
//cout<<"S0"<<endl; //cout<<"S0"<<endl;
m_S0_old=m_S0; m_S0_old=m_S0;
...@@ -800,10 +875,12 @@ namespace FIBRE{ ...@@ -800,10 +875,12 @@ namespace FIBRE{
return rejflag; return rejflag;
}; };
inline void accept_S0(){ inline void accept_S0(){
m_S0_acc++; m_S0_acc++;
} }
inline void reject_S0(){ inline void reject_S0(){
m_S0=m_S0_old; m_S0=m_S0_old;
m_S0_prior=m_S0_old_prior; m_S0_prior=m_S0_old_prior;
...@@ -811,6 +888,9 @@ namespace FIBRE{ ...@@ -811,6 +888,9 @@ namespace FIBRE{
m_S0_rej++; m_S0_rej++;
} }
//Adapt standard deviation of proposal distributions during MCMC execution
//to avoid over-rejection/over-acceptance of MCMC samples
inline void update_proposals(){ inline void update_proposals(){
m_d_prop*=sqrt(float(m_d_acc+1)/float(m_d_rej+1)); m_d_prop*=sqrt(float(m_d_acc+1)/float(m_d_rej+1));
m_d_prop=min(m_d_prop,maxfloat); m_d_prop=min(m_d_prop,maxfloat);
...@@ -831,16 +911,14 @@ namespace FIBRE{ ...@@ -831,16 +911,14 @@ namespace FIBRE{
} }
} }
//Function that performs a single MCMC iteration
void jump( bool can_use_ard=true ){ void jump( bool can_use_ard=true ){
// cout<<"modelnum "<<m_modelnum<<endl; // cout<<"modelnum "<<m_modelnum<<endl;
// cout<<"fmodelnum "<<m_fibres[0].get_modelnum()<<endl; // cout<<"fmodelnum "<<m_fibres[0].get_modelnum()<<endl;
// cout<<"fmodelnum2 "<<m_fibres[1].get_modelnum()<<endl; // cout<<"fmodelnum2 "<<m_fibres[1].get_modelnum()<<endl;
if(!propose_d()){ if(!propose_d()){ //Try d
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
compute_energy(); compute_energy();
...@@ -856,7 +934,7 @@ namespace FIBRE{ ...@@ -856,7 +934,7 @@ namespace FIBRE{
reject_d(); reject_d();
} }
if(m_modelnum==2){ if(m_modelnum==2){ //Try d_std
if(!propose_d_std()){ if(!propose_d_std()){
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
...@@ -873,7 +951,7 @@ namespace FIBRE{ ...@@ -873,7 +951,7 @@ namespace FIBRE{
reject_d_std(); reject_d_std();
} }
} }
if(!propose_S0()){ if(!propose_S0()){ //Try S0
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
compute_energy(); compute_energy();
...@@ -888,14 +966,14 @@ namespace FIBRE{ ...@@ -888,14 +966,14 @@ namespace FIBRE{
reject_S0(); reject_S0();
} }
for(unsigned int f=0;f<m_fibres.size();f++){ for(unsigned int f=0;f<m_fibres.size();f++){ //For each fibre
//cout<<"pre th "<<m_fibres[f].get_th()<<endl; //cout<<"pre th "<<m_fibres[f].get_th()<<endl;
//cout << m_fibres[f].getSignal()<<endl; //cout << m_fibres[f].getSignal()<<endl;
//cout<<"prior"<<m_prior_en<<endl; //cout<<"prior"<<m_prior_en<<endl;
//cout<<"lik"<<m_likelihood_en<<endl; //cout<<"lik"<<m_likelihood_en<<endl;
//cout<<"energy"<<m_energy<<endl; //cout<<"energy"<<m_energy<<endl;
//cout <<endl; //cout <<endl;
if(!m_fibres[f].propose_th()){ if(!m_fibres[f].propose_th()){ //Try theta
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
...@@ -923,7 +1001,7 @@ namespace FIBRE{ ...@@ -923,7 +1001,7 @@ namespace FIBRE{
m_fibres[f].reject_th(); m_fibres[f].reject_th();
} }
if(!m_fibres[f].propose_ph()){ if(!m_fibres[f].propose_ph()){ //Try phi
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
compute_energy(); compute_energy();
...@@ -938,7 +1016,7 @@ namespace FIBRE{ ...@@ -938,7 +1016,7 @@ namespace FIBRE{
m_fibres[f].reject_ph(); m_fibres[f].reject_ph();
} }
if(!m_fibres[f].propose_f( can_use_ard )){ if(!m_fibres[f].propose_f( can_use_ard )){ //Try f
if(!reject_f_sum()){ if(!reject_f_sum()){
compute_prior(); compute_prior();
compute_likelihood(); compute_likelihood();
...@@ -959,7 +1037,6 @@ namespace FIBRE{ ...@@ -959,7 +1037,6 @@ namespace FIBRE{
m_fibres[f].reject_f(); m_fibres[f].reject_f();
} }
// if(!m_fibres[f].propose_lam()){ // if(!m_fibres[f].propose_lam()){
// compute_prior(); // compute_prior();
// compute_energy(); // compute_energy();
...@@ -974,13 +1051,10 @@ namespace FIBRE{ ...@@ -974,13 +1051,10 @@ namespace FIBRE{
// else{ // else{
// m_fibres[f].reject_lam(); // m_fibres[f].reject_lam();
// } // }
} }
} }
Multifibre& operator=(const Multifibre& rhs){ Multifibre& operator=(const Multifibre& rhs){
m_fibres=rhs.m_fibres; m_fibres=rhs.m_fibres;
m_d=rhs.m_d; m_d=rhs.m_d;
...@@ -1011,13 +1085,14 @@ namespace FIBRE{ ...@@ -1011,13 +1085,14 @@ namespace FIBRE{
m_energy=rhs.m_energy; m_energy=rhs.m_energy;
m_old_energy=rhs.m_old_energy; m_old_energy=rhs.m_old_energy;
m_ardfudge=rhs.m_ardfudge; m_ardfudge=rhs.m_ardfudge;
m_iso_Signal=rhs.m_iso_Signal; //STAM: Equate the Iso_signal vectors as well
m_iso_Signal_old=rhs.m_iso_Signal_old;
return *this; return *this;
} }
}; };
} }
#endif #endif
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