Commit 762b292a authored by William Clarke's avatar William Clarke
Browse files

Clean up of utils misc.

parent a652450f
#!/usr/bin/env python
# misc.py - Various utils
#
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
......@@ -9,26 +7,26 @@
# SHBASECOPYRIGHT
import numpy as np
import scipy.fft
from scipy.signal import butter, lfilter
from scipy.interpolate import interp1d
import itertools as it
from .constants import H2O_PPM_TO_TMS
def ppm2hz(cf, ppm, shift=True, shift_amount=H2O_PPM_TO_TMS):
"""Convert ppm scale to frequency scale with optional shift."""
if shift:
return (ppm-shift_amount)*cf*1E-6
return (ppm - shift_amount) * cf * 1E-6
else:
return (ppm)*cf*1E-6
return (ppm) * cf * 1E-6
def hz2ppm(cf, hz, shift=True, shift_amount=H2O_PPM_TO_TMS):
"""Convert frequency scale to frequency scale with optional shift."""
if shift:
return 1E6 * hz/cf + shift_amount
return 1E6 * hz / cf + shift_amount
else:
return 1E6 * hz/cf
return 1E6 * hz / cf
def FIDToSpec(FID, axis=0):
......@@ -47,10 +45,10 @@ def FIDToSpec(FID, axis=0):
ss[axis] = slice(0, 1)
ss = tuple(ss)
FID[ss] *= 0.5
out = scipy.fft.fftshift(scipy.fft.fft(FID,
out = np.fft.fftshift(np.fft.fft(FID,
axis=axis,
norm='ortho'),
axes=axis)
axes=axis)
FID[ss] *= 2
return out
......@@ -66,9 +64,9 @@ def SpecToFID(spec, axis=0):
Returns:
x (np.array) : array of FIDs
"""
fid = scipy.fft.ifft(scipy.fft.ifftshift(spec,
axes=axis),
axis=axis, norm='ortho')
fid = np.fft.ifft(np.fft.ifftshift(spec,
axes=axis),
axis=axis, norm='ortho')
ss = [slice(None) for i in range(fid.ndim)]
ss[axis] = slice(0, 1)
ss = tuple(ss)
......@@ -77,12 +75,21 @@ def SpecToFID(spec, axis=0):
def calculateAxes(bandwidth, centralFrequency, points, shift):
dwellTime = 1/bandwidth
"""Generate time, frequency and ppm axes.
Input:
bandwidth: Bandwidth in Hz.
centralFrequency: Spectroscopy frequency in Hz.
points: Number of time domain points.
shift: Shift on ppm axis.
Returns:
Dict with 'time', 'freq', 'ppm', 'ppmshift' fields.
"""
dwellTime = 1 / bandwidth
timeAxis = np.linspace(dwellTime,
dwellTime * points,
points)
frequencyAxis = np.linspace(-bandwidth/2,
bandwidth/2,
frequencyAxis = np.linspace(-bandwidth / 2,
bandwidth / 2,
points)
ppmAxis = hz2ppm(centralFrequency,
frequencyAxis,
......@@ -101,117 +108,87 @@ def calculateAxes(bandwidth, centralFrequency, points, shift):
def checkCFUnits(cf, units='Hz'):
""" Check the units of central frequency and adjust if required."""
# Assume cf in Hz > 1E5, if it isn't assume that user has passed in MHz
if cf<1E5:
if units.lower()=='hz':
if cf < 1E5:
if units.lower() == 'hz':
cf *= 1E6
elif units.lower()=='mhz':
elif units.lower() == 'mhz':
pass
else:
raise ValueError('Only Hz or MHz defined')
else:
if units.lower()=='hz':
if units.lower() == 'hz':
pass
elif units.lower()=='mhz':
elif units.lower() == 'mhz':
cf /= 1E6
else:
raise ValueError('Only Hz or MHz defined')
return cf
def filter(mrs,FID,ppmlim,filter_type='bandpass'):
def filter(mrs, FID, ppmlim, filter_type='bandpass'):
"""
Filter in/out frequencies defined in ppm
Parameters
----------
mrs : MRS Object
FID : array-like
temporal signal to filter
ppmlim : float or tuple
ppmlim : float or tuple
filter_type: {'lowpass','highpass','bandpass', 'bandstop'}
default type is 'bandstop')
Outputs
-------
numpy array
numpy array
"""
# Sampling frequency (Hz)
fs = 1/mrs.dwellTime
fs = 1 / mrs.dwellTime
nyq = 0.5 * fs
#first,last = mrs.ppmlim_to_range(ppmlim)
#f1,f2 = np.abs(mrs.frequencyAxis[first]),np.abs(mrs.frequencyAxis[last])
#if f1>f2:
# f1,f2=f2,f1
#wn = [f1/nyq,f2/nyq]
f1 = np.abs(ppm2hz(mrs.centralFrequency, ppmlim[0]) / nyq)
f2 = np.abs(ppm2hz(mrs.centralFrequency, ppmlim[1]) / nyq)
f1 = np.abs(ppm2hz(mrs.centralFrequency,ppmlim[0])/ nyq)
f2 = np.abs(ppm2hz(mrs.centralFrequency,ppmlim[1])/ nyq)
if f1 > f2:
f1, f2 = f2, f1
wn = [f1, f2]
if f1>f2:
f1,f2=f2,f1
wn = [f1,f2]
#print(wn)
order = 6
b,a = butter(order, wn, btype=filter_type)
b, a = butter(order, wn, btype=filter_type)
y = lfilter(b, a, FID)
return y
def resample_ts(ts,dwell,new_dwell):
"""
Temporal resampling
Parameters
----------
ts : matrix (rows=time axis)
dwell and new_dwell : float
Returns
-------
matrix
def ts_to_ts(old_ts, old_dt, new_dt, new_n):
"""Temporal resampling where the new time series has a smaller number of points
Input:
old_ts: Input time-domain data
old_dt: Input dwelltime
new_dt: Output dwelltime
new_n: Output number of points
"""
numPoints = ts.shape[0]
t = np.linspace(dwell,dwell*numPoints,numPoints)-dwell
new_t = np.linspace(new_dwell,new_dwell*numPoints,numPoints)-new_dwell
f = interp1d(t,ts,axis=0)
new_ts = f(new_t)
return new_ts
old_n = old_ts.shape[0]
old_t = np.linspace(old_dt, old_dt * old_n, old_n) - old_dt
new_t = np.linspace(new_dt, new_dt * new_n, new_n) - new_dt
def ts_to_ts(old_ts,old_dt,new_dt,new_n):
"""
Temporal resampling where the new time series has a smaller number of points
"""
old_n = old_ts.shape[0]
old_t = np.linspace(old_dt,old_dt*old_n,old_n)-old_dt
new_t = np.linspace(new_dt,new_dt*new_n,new_n)-new_dt
f = interp1d(old_t,old_ts,axis=0)
f = interp1d(old_t, old_ts, axis=0)
new_ts = f(new_t)
return new_ts
# Numerical differentiation (light)
import numpy as np
#Gradient Function
# Gradient Function
def gradient(x, f):
"""
Calculate f'(x): the numerical gradient of a function
Parameters:
-----------
x : array-like
x : array-like
f : scalar function
Returns:
......@@ -221,24 +198,24 @@ def gradient(x, f):
N = len(x)
gradient = []
for i in range(N):
eps = abs(x[i]) * np.finfo(np.float32).eps
if eps==0:
eps = abs(x[i]) * np.finfo(np.float32).eps
if eps == 0:
eps = 1e-5
xl = np.array(x)
xu = np.array(x)
xl[i] -= eps
xu[i] += eps
fl = f(xl)
fu = f(xu)
gradient.append((fu-fl)/(2*eps))
fu = f(xu)
gradient.append((fu - fl) / (2 * eps))
return np.array(gradient)
#Hessian Matrix
def hessian (x, f):
# Hessian Matrix
def hessian(x, f):
"""
Calculate numerical Hessian of f at x
Parameters:
-----------
x : array-like
......@@ -250,25 +227,26 @@ def hessian (x, f):
"""
N = len(x)
hessian = []
gd_0 = gradient( x, f)
eps = np.linalg.norm(gd_0) * np.finfo(np.float32).eps
if eps==0:
gd_0 = gradient(x, f)
eps = np.linalg.norm(gd_0) * np.finfo(np.float32).eps
if eps == 0:
eps = 1e-5
for i in range(N):
hessian.append([])
xx0 = 1.*x[i]
xx0 = 1. * x[i]
x[i] = xx0 + eps
gd_1 = gradient(x, f)
gd_1 = gradient(x, f)
for j in range(N):
hessian[i].append((gd_1[j,:] - gd_0[j,:])/eps)
x[i] =xx0
hessian[i].append((gd_1[j, :] - gd_0[j, :]) / eps)
x[i] = xx0
return np.asarray(hessian)
def hessian_diag(x,f):
def hessian_diag(x, f):
"""
Calculate numerical second order derivative of f at x
(the diagonal of the Hessian)
Parameters:
-----------
x : array-like
......@@ -279,46 +257,50 @@ def hessian_diag(x,f):
array-like
"""
N = x.size
hess = np.zeros((N,1))
gd_0 = gradient( x, f)
hess = np.zeros((N, 1))
gd_0 = gradient(x, f)
eps = np.linalg.norm(gd_0) * np.finfo(np.float32).eps
if eps==0:
if eps == 0:
eps = 1e-5
for i in range(N):
xx0 = 1.*x[i]
xx0 = 1. * x[i]
x[i] = xx0 + eps
gd_1 = gradient(x, f)
hess[i] = ((gd_1[i] - gd_0[i])/eps)
x[i] =xx0
gd_1 = gradient(x, f)
hess[i] = ((gd_1[i] - gd_0[i]) / eps)
x[i] = xx0
return hess
# Little bit of code for checking the gradients
def check_gradients():
m = np.linspace(0,10,100)
cf = lambda p : np.sum(p[0]*np.exp(-p[1]*m))
x0 = np.random.randn(2)*.1
grad_num = gradient(x0,cf)
E = lambda x : np.sum(np.exp(-x[1]*m))
grad_anal = np.array([E(x0),-x0[0]*np.sum(m*np.exp(-x0[1]*m))])
hess_anal = np.zeros((2,2))
hess_anal[0,1] = -np.sum(m*np.exp(-x0[1]*m))
hess_anal[1,0] = -np.sum(m*np.exp(-x0[1]*m))
hess_anal[1,1] = x0[0]*np.sum(m**2*np.exp(-x0[1]*m))
hess_num = hessian(x0,cf)
hess_diag = hessian_diag(x0,cf)
print('x0 = {}, f(x0) = {}'.format(x0,cf(x0)))
"""Little bit of code for checking the gradients."""
m = np.linspace(0, 10, 100)
def cf(p):
return np.sum(p[0] * np.exp(-p[1] * m))
x0 = np.random.randn(2) * .1
grad_num = gradient(x0, cf)
def E(x):
return np.sum(np.exp(-x[1] * m))
grad_anal = np.array([E(x0), -x0[0] * np.sum(m * np.exp(-x0[1] * m))])
hess_anal = np.zeros((2, 2))
hess_anal[0, 1] = -np.sum(m * np.exp(-x0[1] * m))
hess_anal[1, 0] = -np.sum(m * np.exp(-x0[1] * m))
hess_anal[1, 1] = x0[0] * np.sum(m**2 * np.exp(-x0[1] * m))
hess_num = hessian(x0, cf)
hess_diag = hessian_diag(x0, cf)
print('x0 = {}, f(x0) = {}'.format(x0, cf(x0)))
print('Grad Analytic : {}'.format(grad_anal))
print('Grad Numerical : {}'.format(grad_num))
print('Hess Analytic : {}'.format(hess_anal))
print('Hess Numreical : {}'.format(hess_num))
print('Hess Diag : {}'.format(hess_diag))
def calculate_crlb(x,f,data):
def calculate_crlb(x, f, data):
"""
Calculate Cramer-Rao Lower Bound
This assumes a model of the form data = f(x) + noise
......@@ -335,52 +317,40 @@ def calculate_crlb(x,f,data):
array-like
"""
# estimate noise variance empirically
sig2 = np.var(data-f(x))
grad = gradient(x,f)
crlb = 1/(np.sum(np.abs(grad)**2,axis=1)/sig2)
sig2 = np.var(data - f(x))
grad = gradient(x, f)
crlb = 1 / (np.sum(np.abs(grad)**2, axis=1) / sig2)
return crlb
def calculate_lap_cov(x,f,data,sig2=None):
def calculate_lap_cov(x, f, data, sig2=None):
"""
Calculate approximate covariance using
Fisher information matrix
Assumes forward model is data=f(x)+N(0,sig^2)
Parameters:
x : array-like
f : function
data : array-like
sig2 : optional noise variance
sig2 : optional noise variance
Returns:
2D array
2D array
"""
x = np.asarray(x)
N = x.size
if sig2 is None:
sig2 = np.var(data-f(x))
grad = gradient(x,f)
J = np.concatenate((np.real(grad),np.imag(grad)),axis=1)
P0 = np.diag(np.ones(N)*1E-5)
P = np.dot(J,J.transpose()) / sig2
C = np.linalg.inv(P+P0)
return C
if sig2 is None:
sig2 = np.var(data - f(x))
grad = gradient(x, f)
J = np.concatenate((np.real(grad), np.imag(grad)), axis=1)
P0 = np.diag(np.ones(N) * 1E-5)
P = np.dot(J, J.transpose()) / sig2
C = np.linalg.inv(P + P0)
def calculate_lap_cov_with_grad(x,f,df,args,sig2=None):
grad = df(x,*args)
if sig2 is None:
sig2 = np.var(data-f(x,*args))
J = np.concatenate((np.real(grad),np.imag(grad)),axis=1)
P0 = np.diag(np.ones(N)*1E-5)
P = np.dot(J,J.transpose()) / sig2
C = np.linalg.inv(P+P0)
return C
......@@ -779,9 +749,3 @@ def smooth_FIDs(FIDlist,window):
fid = fid/n
sFIDlist.append(fid)
return sFIDlist
#### Dynamic MRS utils
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