Commit 590c0c38 authored by William Clarke's avatar William Clarke
Browse files

Add docs and tests.

parent 578ddbb4
from fsl_mrs.utils import synthetic as syn
from fsl_mrs.utils.misc import FIDToSpec
from fsl_mrs.utils import mrs_io
from fsl_mrs.core import MRS
import numpy as np
from pathlib import Path
testsPath = Path(__file__).parent
basis_path = testsPath / 'testdata/fsl_mrs/steam_basis'
dynamic_model_path = testsPath / 'testdata/dynamic/fmrs_sin_model.py'
def test_noisecov():
# Create positive semi-definite noise covariance
inputnoisecov = np.random.random((2,2))
inputnoisecov= np.dot(inputnoisecov,inputnoisecov.T)
inputnoisecov = np.random.random((2, 2))
inputnoisecov = np.dot(inputnoisecov, inputnoisecov.T)
testFID,hdr = syn.syntheticFID(coilamps=[1.0,1.0],
coilphase=[0.0,0.0],
noisecovariance=inputnoisecov,
amplitude=[0.0,0.0],
points= 32768)
testFID, hdr = syn.syntheticFID(coilamps=[1.0, 1.0],
coilphase=[0.0, 0.0],
noisecovariance=inputnoisecov,
amplitude=[0.0, 0.0],
points=32768)
outcov = np.cov(np.asarray(testFID))
# Noise cov is for both real and imag, so multiply by 2
assert np.isclose(outcov,2*inputnoisecov,atol=1E-1).all()
assert np.allclose(outcov, 2 * inputnoisecov, atol=1E-1)
def test_syntheticFID():
testFID,hdr = syn.syntheticFID(noisecovariance=[[0.0]],points=16384)
testFID, hdr = syn.syntheticFID(noisecovariance=[[0.0]], points=16384)
# Check FID is sum of lorentzian lineshapes
# anlytical solution
T2 = 1/(hdr['inputopts']['damping'][0])
T2 = 1 / (hdr['inputopts']['damping'][0])
M0 = hdr['inputopts']['amplitude'][0]
f0 = hdr['inputopts']['centralfrequency']*hdr['inputopts']['chemicalshift'][0]
f1 = hdr['inputopts']['centralfrequency']*hdr['inputopts']['chemicalshift'][1]
f0 = hdr['inputopts']['centralfrequency'] * hdr['inputopts']['chemicalshift'][0]
f1 = hdr['inputopts']['centralfrequency'] * hdr['inputopts']['chemicalshift'][1]
f = hdr['faxis']
spec = (M0*T2)/(1+4*np.pi**2*(f0-f)**2*T2**2) +1j*(2*np.pi*M0*(f0-f)*T2**2)/(1+4*np.pi**2*(f0-f)**2*T2**2)
spec += (M0*T2)/(1+4*np.pi**2*(f1-f)**2*T2**2) +1j*(2*np.pi*M0*(f1-f)*T2**2)/(1+4*np.pi**2*(f1-f)**2*T2**2)
spec = (M0 * T2) \
/ (1 + 4 * np.pi**2 * (f0 - f)**2 * T2**2) \
+ 1j * (2 * np.pi * M0 * (f0 - f) * T2**2) \
/ (1 + 4 * np.pi**2 * (f0 - f)**2 * T2**2)
spec += (M0 * T2) \
/ (1 + 4 * np.pi**2 * (f1 - f)**2 * T2**2) \
+ 1j * (2 * np.pi * M0 * (f1 - f) * T2**2) \
/ (1 + 4 * np.pi**2 * (f1 - f)**2 * T2**2)
# Can't quite get the scaling right here.
testSpec = FIDToSpec(testFID[0])
spec /= np.max(np.abs(spec))
testSpec /= np.max(np.abs(testSpec))
assert np.isclose(spec,FIDToSpec(testFID[0]),atol = 1E-2,rtol = 1E0).all()
assert np.allclose(spec, FIDToSpec(testFID[0]), atol=1E-2, rtol=1E0)
def test_syntheticFromBasis():
# TO DO
pass
fid, header, _ = syn.syntheticFromBasisFile(str(basis_path),
ignore=['Scyllo'],
baseline=[0.0, 0.0],
concentrations={'Mac': 2.0},
coilamps=[1.0, 1.0],
coilphase=[0.0, np.pi],
noisecovariance=[[0.1, 0.0], [0.0, 0.1]])
assert fid.shape == (2048, 2)
def test_syntheticFromBasis_baseline():
fid, header, _ = syn.syntheticFromBasisFile(str(basis_path),
baseline=((0.0, 0.0),),
baseline=[0.0, 0.0],
concentrations={'Mac': 2.0},
noisecovariance=[[0.0]])
......@@ -61,7 +77,7 @@ def test_syntheticFromBasis_baseline():
mrs.conj_FID()
fid, header, _ = syn.syntheticFromBasisFile(str(basis_path),
baseline=((1.0, 1.0),),
baseline=[1.0, 1.0],
concentrations={'Mac': 2.0},
noisecovariance=[[0.0]])
......@@ -69,3 +85,48 @@ def test_syntheticFromBasis_baseline():
mrs2.conj_FID()
assert np.allclose(mrs2.get_spec(), mrs.get_spec() + np.complex(1.0, -1.0))
def test_synthetic_spectra_from_model():
_, names, _ = mrs_io.read_basis(str(basis_path))
time_var = np.arange(0, 10)
period = 10.0
time_var = time_var / period
camp = [0] * len(names)
camp[names.index('NAA')] = 2
defined_vals = {'c_0': 'conc',
'c_amp': camp,
'gamma': (20, 0),
'sigma': (20, 0),
'b_intercept': [1, 1, 1, 1, 1, 1],
'b_slope': [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]}
p_noise = {'eps': (0, [0.01])}
p_rel_noise = {'conc': (0, [0.05])}
mrs_list, _, _ = syn.synthetic_spectra_from_model(
str(dynamic_model_path),
time_var,
str(basis_path),
ignore=None,
metab_groups=['Mac'],
ind_scaling=['Mac'],
baseline_order=2,
concentrations={'Mac': 15},
defined_vals=defined_vals,
bandwidth=6000,
points=2048,
baseline_ppm=(0, 5),
param_noise=p_noise,
param_rel_noise=p_rel_noise,
coilamps=[1.0, 1.0],
coilphase=[0.0, np.pi],
noisecovariance=[[0.1, 0.0], [0.0, 0.1]])
assert len(mrs_list) == 10
assert len(mrs_list[0]) == 2
assert mrs_list[0][0].numBasis == len(names)
assert mrs_list[0][0].numPoints == 2048
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
from fsl_mrs.utils.synthetic.synthetic import syntheticFID
from fsl_mrs.utils.synthetic.synthetic_from_basis import syntheticFromBasisFile
from fsl_mrs.utils.synthetic.synthetic_from_dynamic import synthetic_spectra_from_model
......@@ -96,27 +96,32 @@ def syntheticFromBasisFile(basisFile,
ind_scaling=None,
concentrations=None,
baseline=None,
baseline_ppm=None,
broadening=(9.0, 0.0),
shifting=0.0,
coilamps=[1.0],
coilphase=[0.0],
noisecovariance=[[0.1]],
bandwidth=4000.0,
points=2048,
baseline_ppm=None):
points=2048):
""" Create synthetic data from a set of FSL-MRS basis files.
Args:
basisFile (str): path to directory containg basis spectra json files
ignore (list of str, optional): Ignore metabolites in basis set.
metab_groups (list of str, optional): Group metabolites to apply different model parameters to groups.
ind_scaling (list of str, optional): Independently scale basis spectra in the basis set.
concentrations (list or dict or None, optional ): If None, standard concentrations will be used.
If list of same length as basis spectra, then these will be
used.
Pass dict to overide standard values for specific metabolites.
Key should be metabolite name.
baseline (list of floats, optional): Baseline coeeficients, 2 (real, imag) needed per order.
e.g. [1,1,0.1, 0.1] to specifiy a 1st order baseline.
baseline_ppm (tuple, optional): Specify ppm range over which baseline is calculated.
broadening (list of tuples or tuple:floats, optional): Tuple containg a gamma and sigma or a list of tuples
for each basis.
shifting (list of floats or float, optional): Eps shift value or a list for each basis.
baseline (list of floats, optional): Baseline parameters. Not yet implemented
coilamps (list of floats, optional): If multiple coils, specify magnitude scaling.
coilphase (list of floats, optional): If multiple coils, specify phase.
noisecovariance (list of floats, optional): N coils x N coils array of noise variance/covariance.
......
......@@ -20,17 +20,53 @@ def synthetic_spectra_from_model(config_file,
ignore=None,
metab_groups=None,
ind_scaling=None,
baseline_order=0,
concentrations=None,
baseline_order=0,
baseline_ppm=None,
defined_vals=None,
param_noise=None,
bandwidth=4000,
points=2048,
baseline_ppm=None,
param_rel_noise=None,
coilamps=[1.0],
coilphase=[0.0],
noisecovariance=[[0.1]]):
noisecovariance=[[0.1]],
bandwidth=4000,
points=2048):
""" Create synthetic dynamic data from FSL-MRS basis file and dynamic configuration file.
Model parameters may be specified using the defined_vals argument. Otheriwse values are randomly set between
the model defined bounds (or -1 and 1 if bounds not set).
Args:
config_file (str): path to dynamic configuration file.
time_var (list of floats): Dynamic 'time' axis.
basisFile (str): path to directory containg basis spectra json files
ignore (list of str, optional): Ignore metabolites in basis set.
metab_groups (list of str, optional): Group metabolites to apply different model parameters to groups.
ind_scaling (list of str, optional): Independently scale basis spectra in the basis set.
concentrations (list or dict or None, optional ): If None, standard concentrations will be used.
If list of same length as basis spectra, then these will be used.
Pass dict to overide standard values for specific metabolites.
Key should be metabolite name.
baseline_order (int, optional): Order of baseline to simulate.
baseline_ppm (tuple, optional): Specify ppm range over which baseline is calculated.
defined_vals (dict, optional): Model parameters can be specified by adding a key of the same name.
If the value is a single value it will be coppied for all relavent values e.g. each metabolite group.
If the value is the same length as the parameter group then then each value will be assigned.
If set to one of 'conc', 'gamma', 'sigma', 'eps', 'baseline' then standard values will be used.
param_noise (dict, optional): Add fixed Gaussian noise to 'conc', 'gamma', 'sigma', 'eps', 'baseline' at
each timepoint. Specified as (mean, SD) for each key.
param_rel_noise (dict, optional): Add relative Gaussian noise to 'conc', 'gamma', 'sigma', 'eps', 'baseline'
at each timepoint. Specified as (mean, SD) for each key.
coilamps (list of floats, optional): If multiple coils, specify magnitude scaling.
coilphase (list of floats, optional): If multiple coils, specify phase.
noisecovariance (list of floats, optional): N coils x N coils array of noise variance/covariance.
bandwidth (float,optional): Bandwidth of output spectrum in Hz
points (int,optional): Number of points in output spectrum.
Returns:
mrs_list: Numpy array of synthetic MRS objects at each timepoint.
vm: fsl_mrs.utils.dynamic.VariableMapping object.
syn_free_params: Value of model free parameters used in the simulation.
"""
empty_mrs, concentrations = prep_mrs_for_synthetic(basis_file,
points,
bandwidth,
......@@ -41,7 +77,7 @@ def synthetic_spectra_from_model(config_file,
mg = parse_metab_groups(empty_mrs, metab_groups)
model = 'voigt'
_, _, forward, _, p2x = getModelFunctions(model)
_, _, forward, x2p, p2x = getModelFunctions(model)
varNames, sizes = FSLModel_vars(model,
n_basis=empty_mrs.numBasis,
n_groups=max(mg) + 1,
......@@ -145,15 +181,33 @@ def synthetic_spectra_from_model(config_file,
for jdx, k in enumerate(varNames):
if param_noise is None \
or k not in param_noise:
mapped_noise[idx, jdx] = np.zeros(sizes[jdx])
fixed_noise = np.zeros(sizes[jdx])
else:
fixed_noise = rng.normal(loc=param_noise[k][0],
scale=param_noise[k][1],
size=sizes[jdx])
if param_rel_noise is None \
or k not in param_rel_noise:
rel_noise = np.zeros(sizes[jdx])
else:
mapped_noise[idx, jdx] = rng.normal(loc=param_noise[k][0],
scale=param_noise[k][1],
size=sizes[jdx])
rel_noise = rng.normal(loc=param_rel_noise[k][0],
scale=param_rel_noise[k][1],
size=sizes[jdx]) \
* np.asarray(mapped[idx, jdx])
mapped_noise[idx, jdx] = fixed_noise + rel_noise
mrs_list = []
for mm, mn in zip(mapped, mapped_noise):
x = np.concatenate(mm) + np.concatenate(mn)
con, gamma, sigma, eps, phi0, phi1, b = x2p(x,
empty_mrs.numBasis,
max(mg) + 1)
con[con < 0] = 0
gamma[gamma < 0] = 0
sigma[sigma < 0] = 0
x = p2x(con, gamma, sigma, eps, phi0, phi1, b)
syn_fid = synthetic_from_fwd_model(forward,
x,
empty_mrs,
......@@ -163,16 +217,29 @@ def synthetic_spectra_from_model(config_file,
coilamps=coilamps,
coilphase=coilphase,
noisecovariance=noisecovariance)
mrs_out = MRS(FID=syn_fid,
cf=empty_mrs.centralFrequency,
bw=bandwidth,
nucleus='1H',
basis=empty_mrs.basis,
names=empty_mrs.names,
basis_hdr={'centralFrequency': empty_mrs.centralFrequency,
'bandwidth': bandwidth})
mrs_list.append(mrs_out)
if syn_fid.ndim > 1:
coils_mrs = []
for fid in syn_fid.T:
mrs_out = MRS(FID=fid,
cf=empty_mrs.centralFrequency,
bw=bandwidth,
nucleus='1H',
basis=empty_mrs.basis,
names=empty_mrs.names,
basis_hdr={'centralFrequency': empty_mrs.centralFrequency,
'bandwidth': bandwidth})
coils_mrs.append(mrs_out)
mrs_list.append(coils_mrs)
else:
mrs_out = MRS(FID=syn_fid,
cf=empty_mrs.centralFrequency,
bw=bandwidth,
nucleus='1H',
basis=empty_mrs.basis,
names=empty_mrs.names,
basis_hdr={'centralFrequency': empty_mrs.centralFrequency,
'bandwidth': bandwidth})
mrs_list.append(mrs_out)
return mrs_list, vm, syn_free_params
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