Commit a9fa0d04 by William Clarke

### Fixed voigt lineshapes.

parent 8c54e555
 ... ... @@ -40,7 +40,7 @@ def print_params(x,mrs,metab_groups,ref_metab='Cr',scale_factor=1): # New strategy for init def init_params(mrs,baseline,ppmlim): first,last = mrs.ppmlim_to_range(ppmlim) y = misc.FIDToSpec(mrs.FID)[first:last] y = mrs.getSpectrum(ppmlim=ppmlim) y = np.concatenate((np.real(y),np.imag(y)),axis=0).flatten() B = baseline[first:last,:].copy() B = np.concatenate((np.real(B),np.imag(B)),axis=0) ... ... @@ -77,32 +77,6 @@ def init_params(mrs,baseline,ppmlim): return g,e,con,b # def init_FSLModel_old(mrs,metab_groups,baseline_order): # """ # Initialise params of FSLModel # """ # # 1. Find gamma and eps # # 2. Use those to init concentrations # # 3. How about phi0 and phi1? # # 1. gamma/eps # gamma,eps,con = init_gamma_eps(mrs) # # Append # x0 = con # concentrations # g = max(metab_groups)+1 # number of metab groups # x0 = np.append(x0,[gamma]*g) # gamma[0].. # x0 = np.append(x0,[eps]*g) # eps[0].. # x0 = np.append(x0,[0,0]) # phi0 and phi1 # x0 = np.append(x0,[0]*2*(baseline_order+1)) # baseline # return x0 def init_FSLModel(mrs,metab_groups,baseline,ppmlim): """ Initialise params of FSLModel ... ... @@ -120,63 +94,59 @@ def init_FSLModel(mrs,metab_groups,baseline,ppmlim): return x0 # THE BELOW NEEDS TO BE REVISTED IN LIGHT OF THE LORENTZIAN INITIALISATION def init_gamma_sigma_eps(mrs): """ Initialise gamma/sigma/epsilon parameters This is done by summing all the basis FIDs and maximizing the correlation with the data FID after shifting and blurring correlation is calculated in the range [.2,4.2] ppm """ target = mrs.FID[:,None] target = extract_spectrum(mrs,target) b = np.sum(mrs.basis,axis=1)[:,None] def cf(p): gamma = p[0] sigma = p[1] eps = p[2] bs = blur_FID_Voigt(mrs,b,gamma,sigma) bs = shift_FID(mrs,bs,eps) bs = extract_spectrum(mrs,bs) xx = 1-correlate(bs,target) return xx x0 = np.array([1,0,0]) res = minimize(cf, x0, method='Powell') g = res.x[0] s = res.x[1] e = res.x[2] def init_params_voigt(mrs,baseline,ppmlim): first,last = mrs.ppmlim_to_range(ppmlim) y = mrs.getSpectrum(ppmlim=ppmlim) y = np.concatenate((np.real(y),np.imag(y)),axis=0).flatten() B = baseline[first:last,:].copy() B = np.concatenate((np.real(B),np.imag(B)),axis=0) def modify_basis(mrs,gamma,sigma,eps): bs = mrs.basis * np.exp(-(gamma+(sigma**2*mrs.timeAxis)+1j*eps)*mrs.timeAxis) bs = misc.FIDToSpec(bs,axis=0) bs = bs[first:last,:] return np.concatenate((np.real(bs),np.imag(bs)),axis=0) def loss(p): gamma,sigma,eps = np.exp(p[0]),np.exp(p[1]),p[2] basis = modify_basis(mrs,gamma,sigma,eps) desmat = np.concatenate((basis,B),axis=1) beta = np.real(np.linalg.pinv(desmat)@y) beta[:mrs.numBasis] = np.clip(beta[:mrs.numBasis],0,None) # project onto >0 concentration pred = np.matmul(desmat,beta) val = np.mean(np.abs(pred-y)**2) return val return g,s,e x0 = np.array([np.log(1e-5),np.log(1e-5),0]) res = minimize(loss, x0, method='Powell') g,s,e = np.exp(res.x[0]),np.exp(res.x[1]),res.x[2] def init_FSLModel_Voigt(mrs,metab_groups,baseline_order): # get concentrations and baseline params basis = modify_basis(mrs,g,s,e) desmat = np.concatenate((basis,B),axis=1) beta = np.real(np.linalg.pinv(desmat)@y) con = np.clip(beta[:mrs.numBasis],0,None) #con = beta[:mrs.numBasis] b = beta[mrs.numBasis:] return g,s,e,con,b def init_FSLModel_Voigt(mrs,metab_groups,baseline,ppmlim): """ Initialise params of FSLModel Initialise params of FSLModel for Voigt linesahapes """ # 1. Find theta, k and eps # 2. Use those to init concentrations # 3. How about phi0 and phi1? # 1. theta/k/eps gamma,sigma,eps = init_gamma_sigma_eps(mrs) new_basis = mrs.basis*np.exp(-(1j*eps+gamma+mrs.timeAxis*sigma**2)*mrs.timeAxis) gamma,sigma,eps,con,b0 = init_params_voigt(mrs,baseline,ppmlim) data = np.append(np.real(mrs.FID),np.imag(mrs.FID),axis=0) desmat = np.append(np.real(new_basis),np.imag(new_basis),axis=0) con = np.real(np.linalg.pinv(desmat)@data) # Append x0 = con g = max(metab_groups)+1 # number of metab groups x0 = np.append(x0,[gamma]*g) # gamma[0].. x0 = np.append(x0,[sigma]*g) # sigma[0].. x0 = np.append(x0,[eps]*g) # eps[0].. x0 = np.append(x0,[0,0]) # phi0 and phi1 x0 = np.append(x0,[0]*2*(baseline_order+1)) # baseline x0 = np.append(x0,b0) # baseline return x0 ... ... @@ -227,7 +197,7 @@ def get_bounds(num_basis,num_metab_groups,baseline_order,model,method,disableBas bnds = [(0,None)]*num_basis # gamma/sigma/eps bnds.extend([(0,None)]*num_metab_groups) if model == 'Voigt': if model.lower() == 'voigt': bnds.extend([(0,None)]*num_metab_groups) bnds.extend([(None,None)]*num_metab_groups) # phi0,phi1 ... ... @@ -249,7 +219,7 @@ def get_bounds(num_basis,num_metab_groups,baseline_order,model,method,disableBas # gamma/sigma/eps LB.extend([0]*num_metab_groups) UB.extend([MAX]*num_metab_groups) if model == 'Voigt': if model.lower() == 'voigt': LB.extend([0]*num_metab_groups) UB.extend([MAX]*num_metab_groups) LB.extend([MIN]*num_metab_groups) ... ... @@ -280,7 +250,7 @@ def get_fitting_mask(num_basis,num_metab_groups,baseline_order,model, else: mask = [0]*num_basis n = 2*num_metab_groups if model == 'Voigt': if model.lower() == 'voigt': n += num_metab_groups if fit_shape: mask.extend([1]*n) ... ... @@ -310,9 +280,9 @@ def fit_FSLModel(mrs, A simplified version of LCModel """ err_func,grad_func,forward,x2p,p2x = models.getModelFunctions(model) if model == 'lorentzian': if model.lower() == 'lorentzian': init_func = init_FSLModel # initialisation of params elif model == 'voigt': elif model.lower() == 'voigt': init_func = init_FSLModel_Voigt # initialisation of params data = mrs.Spec.copy() # data copied to keep it safe ... ... @@ -399,15 +369,27 @@ def fit_FSLModel(mrs, metab_groups=metab_groups,baseline_order=baseline_order,model=model) x0 = res.params # log-transform positive params con,gamma,eps,phi0,phi1,b = x2p(x0,mrs.numBasis,g) con[con<=0] = 1e-10 gamma[gamma<=0] = 1e-10 logcon,loggamma = np.log(con),np.log(gamma) vbx0 = p2x(logcon,loggamma,eps,phi0,phi1,b) # run VB fitting vbmodel = vb.NonlinVB(models.FSLModel_forward_vb) if model.lower()=='lorentzian': # log-transform positive params con,gamma,eps,phi0,phi1,b = x2p(x0,mrs.numBasis,g) con[con<=0] = 1e-10 gamma[gamma<=0] = 1e-10 logcon,loggamma = np.log(con),np.log(gamma) vbx0 = p2x(logcon,loggamma,eps,phi0,phi1,b) vbmodel = vb.NonlinVB(models.FSLModel_forward_vb) elif model.lower()=='voigt': # log-transform positive params con,gamma,sigma,eps,phi0,phi1,b = x2p(x0,mrs.numBasis,g) con[con<=0] = 1e-10 gamma[gamma<=0] = 1e-10 sigma[sigma<=0] = 1e-10 logcon,loggamma,logsigma = np.log(con),np.log(gamma),np.log(sigma) vbx0 = p2x(logcon,loggamma,logsigma,eps,phi0,phi1,b) vbmodel = vb.NonlinVB(models.FSLModel_forward_vb_voigt) datasplit = np.concatenate((np.real(data[first:last]),np.imag(data[first:last]))) args = [freq,time,basis,B,metab_groups,g,first,last] ... ... @@ -418,8 +400,12 @@ def fit_FSLModel(mrs, args=args) # de-log logcon,loggamma,eps,phi0,phi1,b = x2p(res_vb.x,mrs.numBasis,g) x = p2x(np.exp(logcon),np.exp(loggamma),eps,phi0,phi1,b) if model.lower()=='lorentzian': logcon,loggamma,eps,phi0,phi1,b = x2p(res_vb.x,mrs.numBasis,g) x = p2x(np.exp(logcon),np.exp(loggamma),eps,phi0,phi1,b) elif model.lower()=='voigt': logcon,loggamma,logsigma,eps,phi0,phi1,b = x2p(res_vb.x,mrs.numBasis,g) x = p2x(np.exp(logcon),np.exp(loggamma),np.exp(logsigma),eps,phi0,phi1,b) # collect results results.loadResults(mrs,x) ... ...
 ... ... @@ -13,146 +13,8 @@ from fsl_mrs.utils.misc import FIDToSpec,SpecToFID # Helper functions for LCModel fitting # faster than utils.misc.FIDtoSpec #def FIDToSpec(FID): # return np.fft.fft(FID,axis=0) #def SpecToFID(FID): # return np.fft.ifft(FID,axis=0) ################################################################################## # ###### First implementation of LCModel-style model. This is already obsolete! # # ###### Use FSLModel instead! # ################################ ################ TIME DOMAIN FUNCTIONS def LCModel_forward(x,nu,t,m): """ x = [con[0],...,con[n-1],gamma,eps,phi0,phi1] """ n = m.shape[1] con = x[:n] gamma = x[n] eps = x[n+1] phi0 = x[n+2] phi1 = x[n+3] M = np.fft.fft(m*np.exp(-(gamma+1j*eps)*t),axis=0) Y_nu = np.exp(-1j*(phi0+phi1*nu)) * (M @ con[:,None]) Y_t = np.fft.ifft(Y_nu,axis=0) return Y_t.flatten() def LCModel_jac(x,nu,t,m,FID,first=None,last=None): n = m.shape[1] con = x[:n] gamma = x[n] eps = x[n+1] phi0 = x[n+2] phi1 = x[n+3] m_term = m*np.exp(-(gamma+1j*eps)*t) phi_term = np.exp(-1j*(phi0+phi1*nu)) Fmet = np.fft.fft(m_term,axis=0) cFmet = Fmet@con[:,None] Y = LCModel_forward(x,nu,t,m) Y = Y[:,None] dYdc = np.fft.ifft(phi_term*Fmet,axis=0) dYdgamma = np.fft.ifft(phi_term*(np.fft.fft( -t*m_term,axis=0)@con[:,None]),axis=0) dYdeps = np.fft.ifft(phi_term*(np.fft.fft(-1j*t*m_term,axis=0)@con[:,None]),axis=0) dYdphi0 = np.fft.ifft(-1j* phi_term*cFmet,axis=0) dYdphi1 = np.fft.ifft(-1j*nu*phi_term*cFmet,axis=0) dY = np.concatenate((dYdc,dYdgamma,dYdeps,dYdphi0,dYdphi1),axis=1) jac = np.real(np.sum(Y*np.conj(dY)+np.conj(Y)*dY - np.conj(FID[:,None])*dY - FID[:,None]*np.conj(dY),axis=0)) return jac def LCModel_err(x,nu,t,m,y,first,last): pred = LCModel_forward(x,nu,t,m) return np.sum(np.absolute(y[first:last]-pred[first:last])**2) def LCModel_approxCovariance(x,n,nu,t,m,FID): noise = error_fun(x,n,nu,t,m,y) J = LCModel_jac(x,n,nu,t,m,FID) H = J[:,None]@J[:,None].T C = np.diag(np.linalg.inverse(H))*noise return C ################ FREQUENCY DOMAIN FUNCTIONS def LCModel_forward_freq(x,nu,t,m): """ x = [con[0],...,con[n-1],gamma,eps,phi0,phi1] """ n = m.shape[1] con = x[:n] gamma = x[n] eps = x[n+1] phi0 = x[n+2] phi1 = x[n+3] M = FIDToSpec(m*np.exp(-(gamma+1j*eps)*t)) Y_nu = np.exp(-1j*(phi0+phi1*nu)) * (M@con[:,None]) return Y_nu.flatten() def LCModel_jac_freq(x,nu,t,m,Spec,first,last): n = m.shape[1] con = x[:n] gamma = x[n] eps = x[n+1] phi0 = x[n+2] phi1 = x[n+3] m_term = m*np.exp(-(gamma+1j*eps)*t) phi_term = np.exp(-1j*(phi0+phi1*nu)) Fmet = FIDToSpec(m_term) Ftmet = FIDToSpec(t*m_term) cFmet = Fmet@con[:,None] Y = LCModel_forward_freq(x,nu,t,m) dYdc = phi_term*Fmet dYdgamma = phi_term*(-Ftmet@con[:,None]) dYdeps = phi_term*(-1j*Ftmet@con[:,None]) dYdphi0 = -1j* phi_term*cFmet dYdphi1 = -1j*nu*phi_term*cFmet # Only compute within a range Y = Y[first:last,None] Spec = Spec[first:last,None] dYdc = dYdc[first:last,:] dYdgamma = dYdgamma[first:last] dYdeps = dYdeps[first:last] dYdphi0 = dYdphi0[first:last] dYdphi1 = dYdphi1[first:last] dY = np.concatenate((dYdc,dYdgamma,dYdeps,dYdphi0,dYdphi1),axis=1) jac = np.real(np.sum(Y*np.conj(dY)+np.conj(Y)*dY - np.conj(Spec)*dY - Spec*np.conj(dY),axis=0)) # bit quicker? #err = Y-Spec #jac = np.sum(dY*np.conj(err) + np.conj(dY)*err,axis=0) return jac #np.real(jac) def LCModel_err_freq(x,nu,t,m,data,first,last): pred = LCModel_forward_freq(x,nu,t,m) err = data[first:last]-pred[first:last] sse = np.real(np.sum(err*np.conj(err))) return sse # ##################### FSL MODEL # Modifications on LCModel to be listed here # def FSLModel_x2param(x,n,g): """ ... ... @@ -601,3 +463,29 @@ def FSLModel_forward_vb(x,nu,t,m,B,G,g,first,last): return np.concatenate((np.real(S),np.imag(S))) def FSLModel_forward_vb_voigt(x,nu,t,m,B,G,g,first,last): n = m.shape[1] # get number of basis functions logcon,loggamma,logsigma,eps,phi0,phi1,b = FSLModel_x2param_Voigt(x,n,g) con = np.exp(logcon) gamma = np.exp(loggamma) sigma = np.exp(logsigma) E = np.zeros((m.shape[0],g),dtype=np.complex) for gg in range(g): E[:,gg] = np.exp(-(1j*eps[gg]+gamma[gg]+t*sigma[gg]**2)*t).flatten() tmp = np.zeros(m.shape,dtype=np.complex) for i,gg in enumerate(G): tmp[:,i] = m[:,i]*E[:,gg] M = FIDToSpec(tmp) S = np.exp(-1j*(phi0+phi1*nu)) * (M@con[:,None]) # add baseline if B is not None: S += B@b[:,None] S = S.flatten()[first:last] return np.concatenate((np.real(S),np.imag(S)))
 ... ... @@ -91,7 +91,7 @@ def calcQCOnResults(mrs,res,resparams,ppmlim): for basemrs in basisMRS: #FWHM baseFWHM = res.getLineShapeParams() fwhm_curr,_,_ = idPeaksCalcFWHM(basemrs,estimatedFWHM=baseFWHM[0],ppmlim=ppmlim) fwhm_curr,_,_ = idPeaksCalcFWHM(basemrs,estimatedFWHM=np.max(baseFWHM[0]),ppmlim=ppmlim) fwhm.append(fwhm_curr) #Basis SNR ... ...
 ... ... @@ -21,9 +21,9 @@ class FitRes(object): def __init__(self,model,method,names,metab_groups,baseline_order,B,ppmlim): """ Short class initilisation """ # Initilise some basic parameters - includes fitting options # Populate parameter names # Populate parameter names self.model = model self.fill_names(names,baseline_order=baseline_order,metab_groups=metab_groups) self.model = model self.method = method self.ppmlim = ppmlim self.baseline_order = baseline_order ... ... @@ -200,6 +200,10 @@ class FitRes(object): for i in range(g): self.params_names.append(f'gamma_{i}') if self.model.lower() == 'voigt': for i in range(g): self.params_names.append(f'sigma_{i}') for i in range(g): self.params_names.append(f"eps_{i}") ... ... @@ -346,11 +350,35 @@ class FitRes(object): gamma *= 1/(np.pi*self.hzperppm) else: raise ValueError('Units must be Hz or ppm.') return gamma combined = gamma return combined,gamma elif self.model=='voigt': raise Exception('TO DO.') if function is None: gamma = np.zeros([self.fitResults.shape[0],self.g]) sigma = np.zeros([self.fitResults.shape[0],self.g]) for g in range(self.g): gamma[:,g] = self.fitResults[f'gamma_{g}'].values sigma[:,g] = self.fitResults[f'sigma_{g}'].values else: gamma = [] sigma = [] for g in range(self.g): gamma.append(self.fitResults[f'gamma_{g}'].apply(function)) sigma.append(self.fitResults[f'sigma_{g}'].apply(function)) gamma = np.asarray(gamma) sigma = np.asarray(sigma) if units.lower() =='hz': gamma *= 1/(np.pi) sigma *= 1/(np.pi) elif units.lower() == 'ppm': gamma *= 1/(np.pi*self.hzperppm) sigma *= 1/(np.pi*self.hzperppm) else: raise ValueError('Units must be Hz or ppm.') # return gamma,sigma combined = gamma/2 + np.sqrt((gamma**2/4.0)+(sigma*2*np.log(2))**2) return combined,gamma,sigma def getBaselineParams(self): """ Return normalised complex baseline parameters.""" ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!