Skip to content
Snippets Groups Projects
Commit f600d0a9 authored by Tim Behrens's avatar Tim Behrens
Browse files

made fibre.h stand alone so can be librarified

parent 4616d0c7
No related branches found
No related tags found
No related merge requests found
...@@ -11,12 +11,11 @@ ...@@ -11,12 +11,11 @@
#include "libprob.h" #include "libprob.h"
#include <cmath> #include <cmath>
#include "miscmaths/miscprob.h" #include "miscmaths/miscprob.h"
#include "xfibresoptions.h"
using namespace std; using namespace std;
using namespace NEWMAT; using namespace NEWMAT;
using namespace MISCMATHS; using namespace MISCMATHS;
using namespace Xfibres;
const float maxfloat=1e10; const float maxfloat=1e10;
const float minfloat=1e-10; const float minfloat=1e-10;
const float maxlogfloat=23; const float maxlogfloat=23;
...@@ -47,6 +46,7 @@ namespace FIBRE{ ...@@ -47,6 +46,7 @@ namespace FIBRE{
float m_lam_old_prior; float m_lam_old_prior;
float m_prior_en; float m_prior_en;
float m_old_prior_en; float m_old_prior_en;
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;
...@@ -62,12 +62,11 @@ namespace FIBRE{ ...@@ -62,12 +62,11 @@ namespace FIBRE{
const ColumnVector& m_alpha; const ColumnVector& m_alpha;
const ColumnVector& m_beta; const ColumnVector& m_beta;
const Matrix& m_bvals; const Matrix& m_bvals;
xfibresOptions& opts;
public: public:
//constructors:: //constructors::
Fibre( const ColumnVector& alpha, const ColumnVector& beta, Fibre( const ColumnVector& alpha, const ColumnVector& beta,
const Matrix& bvals,const float& d): const Matrix& bvals,const float& d,const float& ardfudge=1):
m_d(d), m_alpha(alpha), m_beta(beta), m_bvals(bvals),opts(xfibresOptions::getInstance()){ m_ardfudge(ardfudge), m_d(d), m_alpha(alpha), m_beta(beta), m_bvals(bvals){
m_th=M_PI/2; m_th=M_PI/2;
m_th_old=m_th; m_th_old=m_th;
...@@ -109,11 +108,11 @@ namespace FIBRE{ ...@@ -109,11 +108,11 @@ namespace FIBRE{
} }
Fibre(const ColumnVector& alpha, Fibre(const ColumnVector& alpha,
const ColumnVector& beta, const Matrix& bvals, const float& d, const ColumnVector& beta, const Matrix& bvals, const float& d, const float& ardfudge,
const float& th, const float& ph, const float& f, const float& th, const float& ph, const float& f,
const float& lam, const bool lam_jump=true) : const float& lam, const bool lam_jump=true) :
m_th(th), m_ph(ph), m_f(f), m_lam(lam),m_lam_jump(lam_jump), m_d(d), m_th(th), m_ph(ph), m_f(f), m_lam(lam), m_ardfudge(ardfudge),m_lam_jump(lam_jump), m_d(d),
m_alpha(alpha), m_beta(beta), m_bvals(bvals),opts(xfibresOptions::getInstance()) m_alpha(alpha), m_beta(beta), m_bvals(bvals)
{ {
m_th_old=m_th; m_th_old=m_th;
m_ph_old=m_ph; m_ph_old=m_ph;
...@@ -151,7 +150,7 @@ namespace FIBRE{ ...@@ -151,7 +150,7 @@ namespace FIBRE{
} }
Fibre(const Fibre& rhs): Fibre(const Fibre& rhs):
m_d(rhs.m_d), m_alpha(rhs.m_alpha), m_beta(rhs.m_beta), m_bvals(rhs.m_bvals),opts(xfibresOptions::getInstance()){ m_d(rhs.m_d), m_alpha(rhs.m_alpha), m_beta(rhs.m_beta), m_bvals(rhs.m_bvals){
*this=rhs; *this=rhs;
} }
...@@ -265,7 +264,7 @@ namespace FIBRE{ ...@@ -265,7 +264,7 @@ namespace FIBRE{
return true; return true;
else{ else{
//m_f_prior=-(log(m_lam) + (m_lam-1)*log(1-m_f)); //m_f_prior=-(log(m_lam) + (m_lam-1)*log(1-m_f));
if(opts.no_ard.value() || !can_use_ard ){ if(!can_use_ard ){
m_f_prior=0; m_f_prior=0;
} }
else{ else{
...@@ -275,7 +274,7 @@ namespace FIBRE{ ...@@ -275,7 +274,7 @@ namespace FIBRE{
else else
m_f_prior=0; m_f_prior=0;
} }
m_f_prior=opts.fudge.value()*m_f_prior; m_f_prior=m_ardfudge*m_f_prior;
return false; return false;
} }
} }
...@@ -459,7 +458,8 @@ namespace FIBRE{ ...@@ -459,7 +458,8 @@ namespace FIBRE{
float m_old_likelihood_en; float m_old_likelihood_en;
float m_energy; float m_energy;
float m_old_energy; float m_old_energy;
float m_ardfudge;
const ColumnVector& m_data; const ColumnVector& m_data;
const ColumnVector& m_alpha; const ColumnVector& m_alpha;
const ColumnVector& m_beta; const ColumnVector& m_beta;
...@@ -467,8 +467,8 @@ namespace FIBRE{ ...@@ -467,8 +467,8 @@ namespace FIBRE{
public: public:
Multifibre(const ColumnVector& data,const ColumnVector& alpha, Multifibre(const ColumnVector& data,const ColumnVector& alpha,
const ColumnVector& beta, const Matrix& b, int N ): const ColumnVector& beta, const Matrix& b, int N ,float ardfudge=1):
m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b){ m_ardfudge(ardfudge),m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b){
// initialise(Amat,N); // initialise(Amat,N);
} }
...@@ -480,12 +480,12 @@ namespace FIBRE{ ...@@ -480,12 +480,12 @@ namespace FIBRE{
} }
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,th,ph,f,lam,lam_jump); Fibre fib(m_alpha,m_beta,m_bvals,m_d,m_ardfudge,th,ph,f,lam,lam_jump);
m_fibres.push_back(fib); m_fibres.push_back(fib);
} }
void addfibre(){ void addfibre(){
Fibre fib(m_alpha,m_beta,m_bvals,m_d); Fibre fib(m_alpha,m_beta,m_bvals,m_d,m_ardfudge);
m_fibres.push_back(fib); m_fibres.push_back(fib);
} }
......
...@@ -437,7 +437,7 @@ class xfibresVoxelManager{ ...@@ -437,7 +437,7 @@ class xfibresVoxelManager{
opts(xfibresOptions::getInstance()), opts(xfibresOptions::getInstance()),
m_samples(samples),m_voxelnumber(voxelnumber),m_data(data), m_samples(samples),m_voxelnumber(voxelnumber),m_data(data),
m_alpha(alpha), m_beta(beta), m_bvals(b), m_alpha(alpha), m_beta(beta), m_bvals(b),
m_multifibre(m_data,m_alpha,m_beta,m_bvals,opts.nfibres.value()){ } m_multifibre(m_data,m_alpha,m_beta,m_bvals,opts.nfibres.value(),opts.fudge.value()){ }
void initialise(const Matrix& Amat){ void initialise(const Matrix& Amat){
...@@ -519,7 +519,7 @@ class xfibresVoxelManager{ ...@@ -519,7 +519,7 @@ class xfibresVoxelManager{
void runmcmc(){ void runmcmc(){
int count=0, recordcount=0,sample=1;//sample will index a newmat matrix int count=0, recordcount=0,sample=1;//sample will index a newmat matrix
for( int i =0;i<opts.nburn.value();i++){ for( int i =0;i<opts.nburn.value();i++){
m_multifibre.jump( i > opts.nburn_noard.value() ); m_multifibre.jump( !opts.no_ard.value() );
count++; count++;
if(count==opts.updateproposalevery.value()){ if(count==opts.updateproposalevery.value()){
m_multifibre.update_proposals(); m_multifibre.update_proposals();
...@@ -528,7 +528,7 @@ class xfibresVoxelManager{ ...@@ -528,7 +528,7 @@ class xfibresVoxelManager{
} }
for( int i =0;i<opts.njumps.value();i++){ for( int i =0;i<opts.njumps.value();i++){
m_multifibre.jump(); m_multifibre.jump(!opts.no_ard.value());
count++; count++;
if(opts.verbose.value()) if(opts.verbose.value())
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment