Commit a5dd7d18 authored by William Clarke's avatar William Clarke
Browse files

Priors now stored in constants.

parent 75725368
......@@ -5,48 +5,52 @@
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
# William Clarke <william.clarke@ndcn.ox.ac.uk>
#
# Copyright (C) 2019 University of Oxford
# Copyright (C) 2019 University of Oxford
# SHBASECOPYRIGHT
H2O_PPM_TO_TMS = 4.65 # Shift of water to Tetramethylsilane
H1_gamma = 42.576 # MHz/tesla
H2O_PPM_TO_TMS = 4.65 # Shift of water to Tetramethylsilane
H1_gamma = 42.576 # MHz/tesla
# Concentration scaling parameters
TISSUE_WATER_DENSITY = {'GM':0.78,'WM':0.65,'CSF':0.97}
#TISSUE_WATER_DENSITY reference: Ernst T, Kreis R, Ross BD. Absolute Quantitation of Water and Metabolites in the Human Brain. I. Compartments and Water. Journal of Magnetic Resonance, Series B 1993;102:1–8 doi: 10.1006/jmrb.1993.1055.
TISSUE_WATER_DENSITY = {'GM': 0.78, 'WM': 0.65, 'CSF': 0.97}
# TISSUE_WATER_DENSITY reference: Ernst T, Kreis R, Ross BD.
# Absolute Quantitation of Water and Metabolites in the Human Brain. I.
# Compartments and Water.
# Journal of Magnetic Resonance, Series B 1993;102:1–8
# doi: 10.1006/jmrb.1993.1055.
H2O_MOLECULAR_MASS = 18.01528 # g/mol
H2O_MOLALITY = 55.51E3 # mmol/kg
H2O_MOLALITY = 55.51E3 # mmol/kg
H2O_PROTONS = 2
# T1 parameters
# T1 parameters
# Derived from a survey of the literature. References listed below.
# Metabolite values derived from an average of NAA, Cr and Cho peaks.
STANDARD_T1 = { '3T':{'H2O_WM':0.97, # Ref: 1-6
'H2O_GM':1.50, # Ref: 1-6
'H2O_CSF':4.47, # Ref: 4
'METAB':1.29}, # Ref: 2, 7-9
'7T':{'H2O_WM':1.21, # Ref: 1-6
'H2O_GM':2.05, # Ref: 1-6
'H2O_CSF':4.43, # Ref: 4
'METAB':1.43}} # Ref: 2, 7-9
STANDARD_T1 = {'3T': {'H2O_WM': 0.97, # Ref: 1-6
'H2O_GM': 1.50, # Ref: 1-6
'H2O_CSF': 4.47, # Ref: 4
'METAB': 1.29}, # Ref: 2, 7-9
'7T': {'H2O_WM': 1.21, # Ref: 1-6
'H2O_GM': 2.05, # Ref: 1-6
'H2O_CSF': 4.43, # Ref: 4
'METAB': 1.43}} # Ref: 2, 7-9
# T2 parameters
STANDARD_T2 = { '3T':{'H2O_WM':0.073, # Ref: 1,3,10-11
'H2O_GM':0.088, # Ref: 1,3,10-11
'H2O_CSF':2.030, # Ref: 12
'METAB':0.194}, # Ref: 7-9,13-15
'7T':{'H2O_WM':0.055, # Ref: 1,3,10-11
'H2O_GM':0.050, # Ref: 1,3,10-11
'H2O_CSF':1.050, # Ref: 12
'METAB':0.131}} # Ref: 7-9,13-15
STANDARD_T2 = {'3T': {'H2O_WM': 0.073, # Ref: 1,3,10-11
'H2O_GM': 0.088, # Ref: 1,3,10-11
'H2O_CSF': 2.030, # Ref: 12
'METAB': 0.194}, # Ref: 7-9,13-15
'7T': {'H2O_WM': 0.055, # Ref: 1,3,10-11
'H2O_GM': 0.050, # Ref: 1,3,10-11
'H2O_CSF': 1.050, # Ref: 12
'METAB': 0.131}} # Ref: 7-9,13-15
'''
T1 & T2 References:
1. Stanisz GJ et al. doi: 10.1002/mrm.20605.
2. Ethofer T et al. doi: 10.1002/mrm.10640.
3. Wansapura JP et al. doi: 10.1002/(SICI)1522-2586(199904)9:4<531::AID-JMRI4>3.0.CO;2-L.
3. Wansapura JP et al. doi: 10.1002/
(SICI)1522-2586(199904)9:4<531::AID-JMRI4>3.0.CO;2-L.
4. Rooney WD et al. doi: 10.1002/mrm.21122.
5. Dieringer MA et al. doi: 10.1371/journal.pone.0091318.
6. Wright PJ et al. doi: 10.1007/s10334-008-0104-8.
......@@ -59,4 +63,28 @@ T1 & T2 References:
13. Marjańska M et al. doi: 10.1002/nbm.1754.
14. Träber F et al. doi: 10.1002/jmri.20053.
15. Wyss PO et al. doi: 10.1002/mrm.27067.
'''
\ No newline at end of file
'''
'''
MCMC PRIORS
Modify these values to change the centre or width (SD) of
the Gaussian priors applied to the MCMC optimised fitting routine.
Priors are defined for the lorentzian and voigt (default) models.
Use the disable_mh_priors flag to fit_FSLModel or the
disable_MH_priors command line argument to set all priors to uniform.
'_loc' values are the centre and '_scale' values are the standard
deviation of the Guassian prior.
'''
MCMC_PRIORS = {'lorentzian': {'conc_loc': 0.0, 'conc_scale': 1E0,
'gamma_loc': 5.0, 'gamma_scale': 2.5, # Hz
'eps_loc': 0.0, 'eps_scale': 0.005, # ppm
'phi0_loc': 0.0, 'phi0_scale': 5.0, # degrees
'phi1_loc': 0.0, 'phi1_scale': 1E-5}, # seconds
'voigt': {'conc_loc': 0.0, 'conc_scale': 1E0,
'gamma_loc': 5.0, 'gamma_scale': 2.5, # Hz
'sigma_loc': 5.0, 'sigma_scale': 2.5, # Hz
'eps_loc': 0.0, 'eps_scale': 0.005, # ppm
'phi0_loc': 0.0, 'phi0_scale': 5.0, # degrees
'phi1_loc': 0.0, 'phi1_scale': 1E-5}} # seconds
......@@ -11,12 +11,10 @@
import numpy as np
from fsl_mrs.utils import models, misc
from fsl_mrs.utils.stats import mh,vb,dist
from fsl_mrs.utils.constants import *
from fsl_mrs.core import MRS
from fsl_mrs.utils.stats import mh, vb, dist
from fsl_mrs.utils.results import FitRes
from scipy.optimize import minimize,nnls
from scipy.optimize import minimize
def print_params(x,mrs,metab_groups,ref_metab='Cr',scale_factor=1):
......@@ -344,30 +342,53 @@ def fit_FSLModel(mrs,
def logpr(p):
return np.sum(dist.gauss_logpdf(p,loc=np.zeros_like(p),scale=np.ones_like(p)*1E2))
else:
from fsl_mrs.utils.constants import MCMC_PRIORS
def logpr(p):
def make_prior(param,loc,scale):
return np.sum(dist.gauss_logpdf(param,
loc=loc*np.ones_like(param),
scale=scale*np.ones_like(param)))
prior = 0
if model.lower()=='lorentzian':
con,gamma,eps,phi0,phi1,b = x2p(p,mrs.numBasis,g)
prior += np.sum(dist.gauss_logpdf(con,loc=np.zeros_like(con),scale=np.ones_like(con)*1E0))
prior += np.sum(dist.gauss_logpdf(gamma,loc=np.ones_like(gamma)*5*np.pi,scale=np.ones_like(gamma)*2.5*np.pi))
prior += np.sum(dist.gauss_logpdf(eps,loc=np.zeros_like(eps),scale=np.ones_like(eps)*0.005*(2*np.pi*mrs.centralFrequency/1E6)))
prior += np.sum(dist.gauss_logpdf(phi0,loc=np.zeros_like(phi0),scale=np.ones_like(phi0)*(np.pi*10/180)))
prior += np.sum(dist.gauss_logpdf(phi1,loc=np.zeros_like(phi1),scale=np.ones_like(phi1)*(1E-5*2*np.pi)))
prior += 0
PRIORS = MCMC_PRIORS['lorentzian']
prior += make_prior(con, PRIORS['conc_loc'], PRIORS['conc_scale'])
prior += make_prior(gamma,
PRIORS['gamma_loc']*np.pi,
PRIORS['gamma_scale']*np.pi)
prior += make_prior(eps,
PRIORS['eps_loc']*2*np.pi*mrs.centralFrequency/1E6,
PRIORS['eps_scale']*2*np.pi*mrs.centralFrequency/1E6)
prior += make_prior(phi0,
PRIORS['phi0_loc']*np.pi/180,
PRIORS['phi0_scale']*np.pi/180)
prior += make_prior(phi1,
PRIORS['phi1_loc']*2*np.pi,
PRIORS['phi1_scale']*2*np.pi)
elif model.lower()=='voigt':
con,gamma,sigma,eps,phi0,phi1,b = x2p(p,mrs.numBasis,g)
prior += np.sum(dist.gauss_logpdf(con,loc=np.zeros_like(con),scale=np.ones_like(con)*1E0))
prior += np.sum(dist.gauss_logpdf(gamma,loc=np.ones_like(gamma)*5*np.pi,scale=np.ones_like(gamma)*2.5*np.pi))
prior += np.sum(dist.gauss_logpdf(sigma,loc=np.ones_like(sigma)*5*np.pi,scale=np.ones_like(sigma)*2.5*np.pi))
prior += np.sum(dist.gauss_logpdf(eps,loc=np.zeros_like(eps),scale=np.ones_like(eps)*0.005*(2*np.pi*mrs.centralFrequency/1E6)))
prior += np.sum(dist.gauss_logpdf(phi0,loc=np.zeros_like(phi0),scale=np.ones_like(phi0)*(np.pi*5/180)))
prior += np.sum(dist.gauss_logpdf(phi1,loc=np.zeros_like(phi1),scale=np.ones_like(phi1)*(1E-5*2*np.pi)))
prior += 0
PRIORS = MCMC_PRIORS['voigt']
prior += make_prior(con, PRIORS['conc_loc'], PRIORS['conc_scale'])
prior += make_prior(gamma,
PRIORS['gamma_loc']*np.pi,
PRIORS['gamma_scale']*np.pi)
prior += make_prior(sigma,
PRIORS['sigma_loc']*np.pi,
PRIORS['sigma_scale']*np.pi)
prior += make_prior(eps,
PRIORS['eps_loc']*2*np.pi*mrs.centralFrequency/1E6,
PRIORS['eps_scale']*2*np.pi*mrs.centralFrequency/1E6)
prior += make_prior(phi0,
PRIORS['phi0_loc']*np.pi/180,
PRIORS['phi0_scale']*np.pi/180)
prior += make_prior(phi1,
PRIORS['phi1_loc']*2*np.pi,
PRIORS['phi1_scale']*2*np.pi)
return prior
#loglik = lambda p : np.log(np.linalg.norm(y-forward_mh(p)[first:last]))*numPoints_over_2
#logpr = lambda p : 0
# Setup the fitting
# Init with nonlinear fit
......
Supports Markdown
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