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

Merge branch 'revisions' into 'master'

V1.0.5

See merge request fsl/fsl_mrs!11
parents 6f08010d d7f2361f
stages:
- Static Analysis
- Test
flake8:
image: python:3.7-slim-buster
stage: Static Analysis
before_script:
- python --version
- pip install flake8
script:
- flake8 --max-line-length=120 fsl_mrs
allow_failure: true
pytest:
image: continuumio/miniconda3:latest
stage: Test
variables:
GIT_SUBMODULE_STRATEGY: normal
before_script:
- conda init bash
- source ~/.bashrc
- conda create -n fsl_mrs -y python=3.7
- conda activate fsl_mrs
- conda install -y -c conda-forge -c https://fsl.fmrib.ox.ac.uk/fsldownloads/fslconda/channel/ --file requirements.txt
- conda install -y -c conda-forge pytest
- pip install .
script:
- pytest fsl_mrs/tests
\ No newline at end of file
This document contains the FSL-MRS release history in reverse chronological order.
1.0.5 (Friday 9th October 2020)
-------------------------------
- Extended documentation of hardcoded constants, including MCMC priors.
- Extended documentation of synthetic macromolecules.
- Added flag to MCMC optimise baseline parameters.
1.0.4 (Friday 14th August 2020)
-------------------------------
- Fixed bug in automatic conjugation facility of fsl_mrs_preproc
......
.. _constants:
Constants
=========
The following constants are hardcoded into FSL-MRS and are used in the quantisation module.
Water to tetramethylsilane (TMS) chemical shift ::
H2O_PPM_TO_TMS = 4.65
Gyromagnetic ratio of proton::
H1_gamma = 42.576 # MHz/tesla
Concentration scaling parameters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Tissue water density [0]_::
TISSUE_WATER_DENSITY = {'GM':0.78,'WM':0.65,'CSF':0.97}
Mollecular mass of water::
H2O_MOLECULAR_MASS = 18.01528 # g/mol
Molality of pure water::
H2O_MOLALITY = 55.51E3 # mmol/kg
Number of protons in water::
H2O_PROTONS = 2
Relaxation parameters
*********************
Values are derived from a survey of the literature. References listed below. Values for metabolites are derived from an average of NAA, Cr and Cho peaks.
T1 (all values in seconds)
__________________________
.. csv-table::
:header: "Field strength (T)", "Water: WM", "Water: GM", "Water: CSF" , "Metabolites"
:widths: 10, 10, 10, 10, 10
3, 0.97, 1.50, 4.47, 1.29
7, 1.21, 2.05, 4.43, 1.43
References, 1-6, 1-6, 4, "2, 7-9"
Code definitions::
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_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
References
----------
.. [0] Kreis R, Ernst T, Ross BD. Absolute quantitation of water and metabolites in the human brain. II. Metabolite concentrations. J Magn Reson B. 1993;102:9-19.
.. [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.
.. [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.
.. [7] Mlynárik V et al. doi: 10.1002/nbm.713.
.. [8] Li Y. doi: 10.4172/2155-9937.S1-002.
.. [9] An L et al. doi: 10.1002/mrm.26612.
.. [10] Gelman N et al. doi: 10.1148/radiology.210.3.r99fe41759.
.. [11] Bartha R et al. doi: 10.1002/mrm.10112.
.. [12] Spijkerman JM et al. doi: 10.1007/s10334-017-0659-3.
.. [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
......@@ -167,7 +167,7 @@ Below are detailed explanations of some of the optional arguments in the wrapper
:code:`--ind_scale`
Allow independent scaling of specified basis spectra before fitting. For example this can be used to independently scale empirically measured macromolecules combined with simulated metabolite spectra.
:code:`--disable_MH_priors`
Disable the priors on the MH fitting. The priors are tuned for *in vivo* human brain spectroscopy. Use this option if your spectra has significantly different line widths, phases or large shifts. E.g. in liquid phase phantom or (potentially) pre-clinical systems.
Disable the priors on the MH fitting. The priors are tuned for *in vivo* human brain spectroscopy. Use this option if your spectra has significantly different line widths, phases or large shifts. E.g. in liquid phase phantom or (potentially) pre-clinical systems. Priors can be fine tuned by altering the values in :code:`fsl_mrs.utils.constants`.
The wrapper scripts can also take a configuration file as an input. For example, say we have a text file called :code:`config.txt` which contains the below:
......
......@@ -22,7 +22,9 @@ If this is your first experience with FSL-MRS, get started with the :ref:`Quick
processing
fitting
quantitation
macromolecules
simulation
visualisation
constants
troubleshooting
changelog
.. _macromolecules:
Macromolecules
==============
It is recommended to use empirically derived macromolecular basis spectra for the fitting step. A measured MM spectrum described in JSON format can be added to the basis spectra by following the instructions on the :ref:`basis spectra simulation <sim_mm>` page.
A collection of MM basis spectra is available for a variety of field strengths, species and anatomies at `MRSHub <https://mrshub.org/datasets_mm/>`_.
In situations where an empirically measured macromolecular spectrum is not available FSL-MRS includes methods for quickly generating synthetic MM basis spectra. For details see `Synthetic MM`_. Both methods should not be used simultaneously.
For an in depth discussion of the effects of MM basis spectra choice on fitting performance see [CUDA12]_ and [GIAP19]_.
Synthetic MM
~~~~~~~~~~~~
Syntheic MM basis spectra may be added using the :code:`--add_MM` flag with :code:`fsl_mrs` or :code:`fsl_mrsi`. In the interactive environment the same can be achieved by calling the method :code:`add_MM_peaks` of :code:`fsl_mrs.core.MRS`.
By default this option will add the following basis spectra (in separate metabolite groups) to the basis sets.
.. csv-table::
:header: Peak name, Peak location(s) (ppm), Peak relative amplitude(s), Peak broadening (gamma/sigma)
:widths: 10, 10, 10, 10
MM09, 0.9, 3, 10/20
MM12, 1.2, 2, 10/20
MM14, 1.4, 2, 10/20
MM17, 1.7, 2, 10/20
MM21, "2.08, 2.25, 1.95, 3.0", "1.33, 0.22, 0.33, 0.4", 10/20
Additional peaks may be added int he interactive environment by calling :code:`add_MM_peaks` with optional arguments to override the defaults.
References
==========
.. [CUDA12] Cudalbu C, Mlynárik V, Gruetter R. Handling Macromolecule Signals in the Quantification of the Neurochemical Profile. Journal of Alzheimer’s Disease 2012;31:S101–S115 doi: 10.3233/JAD-2012-120100.
.. [GIAP19] Giapitzakis I-A, Borbath T, Murali‐Manohar S, Avdievich N, Henning A. Investigation of the influence of macromolecules and spline baseline in the fitting model of human brain spectra at 9.4T. Magnetic Resonance in Medicine 2019;81:746–758 doi: 10.1002/mrm.27467.
......@@ -74,11 +74,13 @@ glycine Gly
Metabolites to simulate can be specified on the command line using the :code:`–m` option with a list typed on the command line, with the :code:`–b` option specifying a text file with one metabolite listed per line, or the :code:`–s` option pointing to a spin system json file for custom spin systems.
It is **not recommended** to simulate and use all of the metabolites. A typical list to start with for short echo time spectroscopy might be
:
It is **not recommended** to simulate and use all of the metabolites. A typical list to start with for short echo time spectroscopy might be::
Ala, Asp, GPC, PCh, Cr, PCr, GABA, Glc, Gln, Glu, GSH, Ins, Lac, NAA, NAAG, PE, Tau
.. _sim_mm:
Including macromolecules in your basis set
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
......
......@@ -23,7 +23,7 @@ Troubleshooting hints
1) Basis spectra are inconsistently scaled. For example empirically derived macromolecular basis spectra can be orders of magnitude larger than the other basis spectra. Before fitting, fsl_mrs(i) scales the magnitude of the data and basis spectra to a known range. Relative scales are preserved within the basis spectra. To permit fsl_mrs(i) to apply different scales to individual basis spectra use the :code:`--ind_scale` option with a list of basis names.
2) The data might have parameters unlike a 7T or 3T human *in vivo* brain spectrum. I.e. the spectrum originates from a pre-clinical system or from phantom. In this case the MCMC priors which are suitable for *in vivo* human case can be disabled using the :code:`--disable_MH_priors` option.
2) The data might have parameters unlike a 7T or 3T human *in vivo* brain spectrum. I.e. the spectrum originates from a pre-clinical system or from phantom. In this case the MCMC priors which are suitable for *in vivo* human case can be disabled using the :code:`--disable_MH_priors` option. Priors can be fine tuned by altering the values in :code:`fsl_mrs.utils.constants`.
3. Identifying the correct files for conversion
Raw data files, especially DICOM files can have obscure naming conventions. It can be difficult to determine which files should be converted for use in FSL-MRS. Tools such as gdcmdump from `GDCM <http://gdcm.sourceforge.net/>`_ can help in identifying the scans by giving you access to the DICOM headers.
......@@ -32,5 +32,6 @@ Troubleshooting hints
4. Data looks 'wrong' after conversion
If when using :code:`mrs_vis` you see no signal and just noise try conjugating the data using :code:`fsl_mrs_proc conj` or try expanding the ppm range plotted :code:`--ppmlim -10 10`. If you see a flat line, then conversion failed. The data might be corrupted - did the acquisition complete successfully?
.. image:: data/bad_data.png
:width: 600
\ No newline at end of file
from fsl_mrs.utils.synthetic import syntheticFID
from fsl_mrs.utils.synthetic.synthetic_from_basis import syntheticFromBasisFile
from fsl_mrs.core import MRS
from fsl_mrs.utils.fitting import fit_FSLModel
from fsl_mrs.utils import mrs_io
from pytest import fixture
import numpy as np
from pathlib import Path
testsPath = Path(__file__).parent
basis_path = testsPath / 'testdata/fsl_mrs/steam_basis'
@fixture
def data():
......@@ -55,6 +60,25 @@ def test_fit_FSLModel_Newton(data):
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_fit_FSLModel_lorentzian_Newton(data):
mrs = data[0]
amplitudes = data[1]
metab_groups = [0]*mrs.numBasis
Fitargs = {'ppmlim': [0.2,4.2],
'method': 'Newton',
'baseline_order': -1,
'metab_groups': metab_groups,
'model': 'lorentzian'}
res = fit_FSLModel(mrs,**Fitargs)
fittedconcs = res.getConc(metab = mrs.names)
fittedRelconcs = res.getConc(scaling='internal',metab = mrs.names)
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_fit_FSLModel_MH(data):
mrs = data[0]
......@@ -72,4 +96,53 @@ def test_fit_FSLModel_MH(data):
fittedRelconcs = res.getConc(scaling='internal',metab = mrs.names)
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
\ No newline at end of file
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_fit_FSLModel_lorentzian_MH(data):
mrs = data[0]
amplitudes = data[1]
metab_groups = [0]*mrs.numBasis
Fitargs = {'ppmlim':[0.2,4.2],
'method':'MH',
'baseline_order':-1,
'metab_groups':metab_groups,
'MHSamples':100,
'model': 'lorentzian'}
res = fit_FSLModel(mrs,**Fitargs)
fittedconcs = res.getConc(metab = mrs.names)
fittedRelconcs = res.getConc(scaling='internal',metab = mrs.names)
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_fit_FSLModel_on_invivo_sim():
FIDs,hdr,trueconcs = syntheticFromBasisFile(basis_path,
noisecovariance=[[1E-3]],
broadening=(9.0, 9.0),
concentrations={'Mac':2.0})
basis,names,header = mrs_io.read_basis(basis_path)
mrs = MRS(FID=FIDs,header=hdr,basis=basis,basis_hdr=header[0],names=names)
mrs.processForFitting()
metab_groups = [0]*mrs.numBasis
Fitargs = {'ppmlim':[0.2,4.2],
'method':'MH',
'baseline_order':-1,
'metab_groups':metab_groups,
'MHSamples':100}
res = fit_FSLModel(mrs,**Fitargs)
fittedRelconcs = res.getConc(scaling='internal',metab = mrs.names)
answers = np.asarray(trueconcs)
answers /= (answers[names.index('Cr')]+trueconcs[names.index('PCr')])
assert np.allclose(fittedRelconcs,answers,atol=5E-2)
\ No newline at end of file
......@@ -11,7 +11,7 @@ import numpy as np
# Set up some synthetic data to use
@fixture(scope='module')
def data():
noiseCov = 0.01
noiseCov = 0.001
amplitude = np.asarray([0.5, 0.5, 1.0])*10
chemshift = np.asarray([3.0, 3.05, 2.0])-4.65
lw = [10, 10, 10]
......
......@@ -5,28 +5,86 @@
# 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} # Kreis R, Ernst T, Ross BD. Absolute quantitation of water and metabolites in the human brain. II. Metabolite concentrations. J Magn Reson B. 1993;102:9-19.
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
# 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
# T2 parameters
# 3T from https://onlinelibrary.wiley.com/doi/epdf/10.1002/nbm.3914
# 7T from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4549223/ - except copied CSF
STANDARD_T2 = { '7T':{'H2O_GM':0.05,
'H2O_WM':0.055,
'H2O_CSF':2.55,
'METAB':0.160},
'3T':{'H2O_GM':0.11,
'H2O_WM':0.08,
'H2O_CSF':2.55,
'METAB':0.271}}
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.
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.
7. Mlynárik V et al. doi: 10.1002/nbm.713.
8. Li Y. doi: 10.4172/2155-9937.S1-002.
9. An L et al. doi: 10.1002/mrm.26612.
10. Gelman N et al. doi: 10.1148/radiology.210.3.r99fe41759.
11. Bartha R et al. doi: 10.1002/mrm.10112.
12. Spijkerman JM et al. doi: 10.1007/s10334-017-0659-3.
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.
'''
'''
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):
......@@ -273,10 +271,11 @@ def fit_FSLModel(mrs,
ppmlim=None,
baseline_order=5,
metab_groups=None,
model='lorentzian',
model='voigt',
x0=None,
MHSamples=500,
disable_mh_priors = False,
fit_baseline_mh = False,
vb_iter=50):
"""
A simplified version of LCModel
......@@ -343,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
......@@ -379,7 +401,7 @@ def fit_FSLModel(mrs,
p0 = res.params
LB,UB = get_bounds(mrs.numBasis,g,baseline_order,model,method,disableBaseline=disableBaseline)
mask = get_fitting_mask(mrs.numBasis,g,baseline_order,model,fit_baseline=False)
mask = get_fitting_mask(mrs.numBasis,g,baseline_order,model,fit_baseline=fit_baseline_mh)
# Check that the values initilised by the newton
# method don't exceed these bounds (unlikely but possible with bad data)
......
......@@ -74,7 +74,7 @@ def syntheticFromBasisFile(basisFile,
Returns:
FIDs: Numpy array of synthetic FIDs
outHeader: Header suitable for loading FIDs into MRS object.
concentrations: Final concentration scalinfs
concentrations: Final concentration scalings
"""
basis,names,header = mrs_io.read_basis(basisFile)
......
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