Skip to content
GitLab
Commit 33708a85 authored by Saad Jbabdi's avatar Saad Jbabdi
Browse files

Changes related to dynMRS and VariableMapping

parent 768dca61
Hide whitespace changes
Inline Side-by-side
fsl_mrs/core/MRS.py
View file @ 33708a85
......@@ -169,8 +169,8 @@ class MRS(object):
(cf_MHz > sevent_range[0] and cf_MHz < sevent_range[1]) or \
(cf_MHz > ninefourt_range[0] and cf_MHz < ninefourt_range[1]) or \
(cf_MHz > elevensevent_range[0] and cf_MHz < elevensevent_range[1]):
print(f'Identified as {key} nucleus data.'
f' Esitmated field: {cf_MHz/GYRO_MAG_RATIO[key]} T.')
#print(f'Identified as {key} nucleus data.'
# f' Esitmated field: {cf_MHz/GYRO_MAG_RATIO[key]} T.')
return key
raise ValueError(f'Unidentified nucleus,'
......
fsl_mrs/utils/example.py
View file @ 33708a85
......@@ -37,7 +37,16 @@ def simulated(ID=1):
return mrs
def dMRS():
def dMRS(mouse='mouse1'):
"""
Load Diffusion MRS data from hard-coded location on disk
Args:
mouse: 'mouse1' ... 'mouse10'
Returns:
mrs Object list
list (time variable)
"""
from scipy.io import loadmat
from fsl_mrs.utils.preproc.phasing import phaseCorrect
from fsl_mrs.utils.preproc.align import phase_freq_align
......@@ -52,7 +61,7 @@ def dMRS():
bandwidth = 5000
currentdir = os.path.join(dataPath,'mouse1')
currentdir = os.path.join(dataPath,mouse)
fidList = []
......@@ -70,11 +79,134 @@ def dMRS():
mrsList = []
for fid,b in zip(alignedFids,blist):
fid,_ = shiftToRef(fid,3.027,bandwidth,centralFrequency,ppmlim=(2.9,3.1))
mrs = MRS(FID=fid,cf=centralFrequency,bw=bandwidth,basis=basis,names=names,basis_hdr=header[0])
mrs = MRS(FID=fid,cf=centralFrequency,bw=bandwidth,basis=basis,names=names,basis_hdr=header[0],nucleus='1H')
mrs.check_FID(repair=True)
mrs.check_Basis(repair=True)
mrs.ignore(['Gly'])
mrsList.append(mrs)
mrsList[0].rescaleForFitting()
for i, mrs in enumerate(mrsList):
if i > 0:
mrs.FID *= mrsList[0].scaling['FID']
mrs.basis *= mrsList[0].scaling['basis']
blist = [b / 1e3 for b in blist]
return mrsList,blist
def dMRS_SNR(avg=1):
"""
Load DMRS data with different numbers of averages
Args:
avg: int (1,2,4,8,16,32,64,128)
Returns:
list of MRS objects
bvals
"""
from fsl_mrs.utils import mrs_io
basispath = '/Users/saad/Desktop/Spectroscopy/DMRS/basis_STE_LASER_8_25_50_LacZ.BASIS'
basis,names,Bheader = mrs_io.read_basis(basispath)
# Bheader['ResonantNucleus'] = '1H'
centralFrequency = 500.30
bandwidth = 5000
FIDpath = f'/Users/saad/Desktop/Spectroscopy/DMRS/WT_multi/{avg:03}_avg'
bvals = [20,3020,6000,10000,20000,30000,50000]
MRSlist = []
for b in bvals:
FID,FIDheader = mrs_io.read_FID(FIDpath+f'/b_{b:05}.nii.gz')
MRSArgs = {'header': FIDheader, 'basis': basis, 'names': names, 'basis_hdr': Bheader[0]}
mrs = MRS(FID=FID, **MRSArgs)
MRSlist.append(mrs)
MRSlist[0].rescaleForFitting()
for i,mrs in enumerate(MRSlist):
if i>0:
mrs.FID *= MRSlist[0].scaling['FID']
mrs.basis *= MRSlist[0].scaling['basis']
bvals = [b / 1e3 for b in bvals]
return MRSlist,bvals
def FMRS(smooth=False):
"""
Load Functional MRS data from hard-coded location on disk
Args:
smooth: bool
Returns:
mrs Object list
list (time variable)
"""
import glob
from numpy import repeat
folder = '/Users/saad/Desktop/Spectroscopy/FMRS/Jacob_data'
filelist = glob.glob(os.path.join(folder, 'RAWFORMAT', 'run_???.RAW'))
FIDlist = []
for file in filelist:
FID, FIDheader = mrs_io.read_FID(file)
FIDlist.append(FID)
basisfile = os.path.join(folder, '7T_slaser36ms_2013_oxford_tdcslb1_ivan.BASIS')
basis, names, basisheader = mrs_io.read_basis(basisfile)
# # Resample basis
from fsl_mrs.utils import misc
basis = misc.ts_to_ts(basis, basisheader[0]['dwelltime'], FIDheader['dwelltime'], FID.shape[0])
if smooth:
sFIDlist = misc.smooth_FIDs(FIDlist, window=5)
else:
sFIDlist = FIDlist
MRSargs = {'names': names, 'basis': basis, 'basis_hdr': basisheader[0], 'bw': FIDheader['bandwidth'],
'cf': FIDheader['centralFrequency']}
mrsList = []
for fid in sFIDlist:
mrs = MRS(FID=fid, **MRSargs)
mrsList.append(mrs)
# Stimulus variable
stim = repeat([0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], 16, axis=0).flatten()
return mrsList,stim
def MPRESS(noise=1):
"""
Load MegaPress edited MRS data from hard-coded location on disk
Args:
noise: float (noise std)
Returns:
mrs Object list
list (time variable)
"""
from fsl_mrs.utils import mrs_io
from fsl_mrs.utils.synthetic.synthetic_from_basis import syntheticFromBasisFile
mpress_on = '/Users/saad/Desktop/Spectroscopy/mpress_basis/ON'
mpress_off = '/Users/saad/Desktop/Spectroscopy/mpress_basis/OFF'
basis,names,basis_hdr = mrs_io.read_basis(mpress_on)
FIDs,header,conc = syntheticFromBasisFile(mpress_on,noisecovariance=[[noise]])
mrs1 = MRS(FID=FIDs,header=header,basis=basis,basis_hdr=basis_hdr[0],names=names)
mrs1.check_FID(repair=True)
mrs1.Spec = misc.FIDToSpec(mrs1.FID)
mrs1.check_Basis(repair=True)
basis,names,basis_hdr = mrs_io.read_basis(mpress_off)
FIDs,header,conc = syntheticFromBasisFile(mpress_off,noisecovariance=[[noise]])
mrs2 = MRS(FID=FIDs,header=header,basis=basis,basis_hdr=basis_hdr[0],names=names)
mrs2.check_FID(repair=True)
mrs2.Spec = misc.FIDToSpec(mrs2.FID)
mrs2.check_Basis(repair=True)
return [mrs1,mrs2],[0,1]
\ No newline at end of file
fsl_mrs/utils/fitting.py
View file @ 33708a85
......@@ -330,12 +330,12 @@ def fit_FSLModel(mrs,
results.loadResults(mrs,x0)
elif method == 'MH':
forward_mh = lambda p : forward(p,freq,time,basis,B,metab_groups,g)
forward_mh = lambda p : forward(p,freq,time,basis,B,metab_groups,g)[first:last]
numPoints_over_2 = (last-first)/2.0
y = data[first:last]
y = data[first:last]
def loglik(p):
return np.log(np.linalg.norm(y-forward_mh(p)[first:last]))*numPoints_over_2
return np.log(np.linalg.norm(y-forward_mh(p)))*numPoints_over_2
if disable_mh_priors:
......@@ -397,13 +397,13 @@ def fit_FSLModel(mrs,
res = fit_FSLModel(mrs,method='Newton',ppmlim=ppmlim,
metab_groups=metab_groups,baseline_order=baseline_order,model=model)
baseline_order = 0 if disableBaseline else baseline_order
# Create maks and bounds for MH fit
# Create masks and bounds for MH fit
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=fit_baseline_mh)
# Check that the values initilised by the newton
# Check that the values initialised by the newton
# method don't exceed these bounds (unlikely but possible with bad data)
for i,(p, u, l) in enumerate(zip(p0, UB, LB)):
if p>u:
......
fsl_mrs/utils/misc.py
View file @ 33708a85
......@@ -13,6 +13,8 @@ import scipy.fft
from scipy.signal import butter, lfilter
from scipy.interpolate import interp1d
import itertools as it
from scipy.optimize import minimize
from functools import partial
from .constants import H2O_PPM_TO_TMS
......@@ -48,8 +50,8 @@ def FIDToSpec(FID, axis=0):
ss = tuple(ss)
FID[ss] *= 0.5
out = scipy.fft.fftshift(scipy.fft.fft(FID,
axis=axis,
norm='ortho'),
axis=axis,
norm='ortho'),
axes=axis)
FID[ss] *= 2
return out
......@@ -67,8 +69,8 @@ def SpecToFID(spec, axis=0):
x (np.array) : array of FIDs
"""
fid = scipy.fft.ifft(scipy.fft.ifftshift(spec,
axes=axis),
axis=axis, norm='ortho')
axes=axis),
axis=axis, norm='ortho')
ss = [slice(None) for i in range(fid.ndim)]
ss[axis] = slice(0, 1)
ss = tuple(ss)
......@@ -779,3 +781,329 @@ def smooth_FIDs(FIDlist,window):
fid = fid/n
sFIDlist.append(fid)
return sFIDlist
#### Dynamic MRS utils
# THE MAIN MAPPING CLASS
# Class responsible for variable mapping
class VariableMapping(object):
def __init__(self,
param_names,
param_sizes,
time_variable,
config_file):
"""
Variable Mapping Class Constructor
Mapping betwee free and mapped:
Mapped = TxN matrix
Mapped[i,j] = float or 1D-array of floats with size param_sizes[j]
Parameters
----------
param_names : list
param_sizes : list
time_variale : array-like
config_file : string
"""
self.time_variable = np.asarray(time_variable)
self.ntimes = self.time_variable.shape[0]
self.mapped_names = param_names
self.mapped_nparams = len(self.mapped_names)
self.mapped_sizes = param_sizes
from runpy import run_path
settings = run_path(config_file)
self.Parameters = settings['Parameters']
for name in self.mapped_names:
if name not in self.Parameters:
self.Parameters[name] = 'fixed'
self.fcns = {}
for key in settings:
if callable(settings[key]):
self.fcns[key] = settings[key]
if 'Bounds' in settings:
self.Bounds = self.create_constraints(settings['Bounds'])
else:
self.Bounds = self.create_constraints(None)
self.nfree = self.calc_nfree()
def __str__(self):
OUT = '-----------------------\n'
OUT += 'Variable Mapping Object\n'
OUT += '-----------------------\n'
OUT += f'Number of Mapped param groups = {len(self.mapped_names)}\n'
OUT += f'Number of Mapped params = {sum(self.mapped_sizes)}\n'
OUT += f'Number of Free params = {self.nfree}\n'
OUT += f'Number of params if all indep = {sum(self.mapped_sizes)*self.ntimes}\n'
OUT += 'Dynamic functions\n'
for param_name in self.mapped_names:
beh = self.Parameters[param_name]
OUT += f'{param_name} \t {beh}\n'
return OUT
def calc_nfree(self):
"""
Calculate number of free parameters based on mapped behaviour
Returns
-------
int
"""
N = 0
for index,param in enumerate(self.mapped_names):
beh = self.Parameters[param]
if (beh == 'fixed'):
N += self.mapped_sizes[index]
elif (beh == 'variable'):
N += self.ntimes*self.mapped_sizes[index]
else:
if 'dynamic' in beh:
N+= len(beh['params'])*self.mapped_sizes[index]
return N
def create_constraints(self,bounds):
"""
Create list of constraints to be used in optimization
Parameters:
-----------
bounds : dict {param:bounds}
Returns
-------
list
"""
if bounds is None:
return [(None,None)]*self.calc_nfree()
if not isinstance(bounds,dict):
raise(Exception('Input should either be a dict or None'))
b = [] # list of bounds
for index,name in enumerate(self.mapped_names):
psize = self.mapped_sizes[index]
if (self.Parameters[name] == 'fixed'):
# check if there are bound on this param
if name in bounds:
for s in range(psize):
b.append(bounds[name])
else:
for s in range(psize):
b.append((None,None))
elif (self.Parameters[name] == 'variable'):
for t in range(self.ntimes):
for s in range(psize):
if name in bounds:
b.append(bounds[name])
else:
b.append((None,None))
else:
if 'dynamic' in self.Parameters[name]:
pnames = self.Parameters[name]['params']
for s in range(psize):
for p in pnames:
if p in bounds:
b.append(bounds[p])
else:
b.append((None,None))
return b
def mapped_from_list(self,p):
"""
Converts list of params into Mapped by repeating over time
Parameters
----------
p : list
Returns
-------
2D array
"""
if isinstance(p,list):
p = np.asarray(p)
if (p.ndim==1):
p = np.repeat(p[None,:],self.ntimes,0)
return p
def create_free_names(self):
"""
create list of names for free params
Returns
-------
list of strings
"""
names = []
for index,param in enumerate(self.mapped_names):
beh = self.Parameters[param]
if (beh == 'fixed'):
name = [f'{param}_{x}' for x in range(self.mapped_sizes[index])]
names.extend(name)
elif (beh == 'variable'):
name = [f'{param}_{x}_t{t}' for x in range(self.mapped_sizes[index]) for t in range(self.ntimes)]
names.extend(name)
else:
if 'dynamic' in beh:
dyn_name = self.Parameters[param]['params']
name = [f'{param}_{y}_{x}' for x in range(self.mapped_sizes[index]) for y in dyn_name]
names.extend(name)
return names
def free_to_mapped(self,p):
"""
Convert free into mapped params over time
fixed params get copied over time domain
variable params are indep over time
dynamic params are mapped using dyn model
Parameters
----------
p : 1D array
Returns
-------
2D array (time X params)
"""
# Check input
if (p.size != self.nfree):
raise(Exception(f'Input free params does not have expected number of entries. Found {p.size}, expected {self.nfree}'))
# Mapped params is time X nparams (each param is an array of params)
mapped_params = np.empty((self.ntimes,self.mapped_nparams),dtype=object)
counter = 0
for index,name in enumerate(self.mapped_names):
nmapped = self.mapped_sizes[index]
if (self.Parameters[name] == 'fixed'): # repeat param over time
for t in range(self.ntimes):
mapped_params[t,index] = p[counter:counter+nmapped]
counter += nmapped
elif (self.Parameters[name] == 'variable'): # copy one param for each time point
for t in range(self.ntimes):
mapped_params[t,index] = p[counter+t*nmapped:counter+t*nmapped+nmapped]
counter += nmapped
else:
if 'dynamic' in self.Parameters[name]:
# Generate time courses
func_name = self.Parameters[name]['dynamic']
nfree = len(self.Parameters[name]['params'])
mapped = np.zeros((self.ntimes,nmapped))
for i in range(nmapped):
params = p[counter:counter+nfree]
mapped[:,i] = self.fcns[func_name](params,self.time_variable)
counter += nfree
for t in range(self.ntimes):
mapped_params[t,index] = mapped[t,:]
else:
raise(Exception("Unknown Parameter type - should be one of 'fixed', 'variable', {'dynamic'}"))
return mapped_params
def print_free(self,x):
"""
Print free params and their names
"""
print(dict(zip(vm.create_free_names(),x)))
def check_bounds(self,x,tol=1e-10):
"""
Check that bounds apply and return corrected x
"""
if self.Bounds is None:
return x
for i,b in enumerate(self.Bounds):
LB = b[0] if b[0] is not None else -np.inf
UB = b[1] if b[1] is not None else np.inf
if (x[i] < LB):
x[i] = LB+tol
if (x[i] > UB):
x[i] = UB-tol
return x
# This function may 'invert' the dynamic mapping
# if the input params are from a single timepoint it assumes constant
def mapped_to_free(self,p):
"""
Convert mapped params to free (e.g. to initialise the free params)
fixed and variable params are simply copied
dynamic params are converted by inverting the dyn model with Scipy optimize
Parameters
----------
p : 2D array (time X params)
Returns
-------
1D array
"""
# Check input
p = self.mapped_from_list(p)
if (p.shape != (self.ntimes,self.mapped_nparams)):
raise(Exception(f'Input mapped params does not have expected number of entries. Found {p.shape}, expected {(self.ntimes,self.mapped_nparams)}'))
free_params = np.empty(self.nfree)
counter = 0
for index,name in enumerate(self.mapped_names):
psize = self.mapped_sizes[index]
if (self.Parameters[name] == 'fixed'):
free_params[counter:counter+psize] = p[0,index]
counter += psize
elif (self.Parameters[name] == 'variable'):
for t in range(self.ntimes):
free_params[counter:counter+psize] = p[t,index]
counter += psize
else:
if 'dynamic' in self.Parameters[name]:
func_name = self.Parameters[name]['dynamic']
time_var = self.time_variable
func = partial(self.fcns[func_name],t=time_var)
nfree = len(self.Parameters[name]['params'])
pp = np.stack(p[:,index][:],axis=0)
for ppp in range(pp.shape[1]):
def loss(x):
pred = func(x)
return np.mean((pp[:,ppp]-pred)**2)
bounds = self.Bounds[counter:counter+nfree]
vals = minimize(loss,
np.zeros(len(self.Parameters[name]['params'])),
method='TNC',bounds=bounds).x
free_params[counter:counter+nfree] = vals
counter += nfree
else:
raise(Exception("Unknown Parameter type - should be one of 'fixed', 'variable', {'dynamic'}"))
return free_params
fsl_mrs/utils/models.py
View file @ 33708a85
......@@ -15,6 +15,22 @@ from fsl_mrs.utils.misc import FIDToSpec,SpecToFID
# ##################### FSL MODEL
def FSLModel_vars(model='voigt'):
"""
Print out parameter names as a list of strings
Args:
model: str (either 'lorientzian' or 'voigt'
Returns:
list of strings
"""
if model == 'lorentzian':
var_names = ['conc', 'gamma', 'eps', 'Phi_0', 'Phi_1', 'baseline']
elif model == 'voigt':
var_names = ['conc', 'gamma', 'sigma', 'eps', 'Phi_0', 'Phi_1', 'baseline']
else:
raise(Exception('model must be either "voigt" or "lorentzian"'))
return var_names
def FSLModel_x2param(x,n,g):
"""
......@@ -83,9 +99,10 @@ def FSLModel_forward(x,nu,t,m,B,G,g):
E[:,gg] = np.exp(-(1j*eps[gg]+gamma[gg])*t).flatten()
# E = np.exp(-(1j*eps+gamma)*t) # THis is actually slower! But maybe more optimisable longterm with numexpr or numba
tmp = np.zeros(m.shape,dtype=np.complex)
for i,gg in enumerate(G):
tmp[:,i] = m[:,i]*E[:,gg]
#tmp = np.zeros(m.shape,dtype=np.complex)
#for i,gg in enumerate(G):
# tmp[:,i] = m[:,i]*E[:,gg]
tmp = m*E[:,G]
M = FIDToSpec(tmp)
S = np.exp(-1j*(phi0+phi1*nu)) * (M@con[:,None])
......