#import matplotlib #matplotlib.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) #matplotlib.rc('text', usetex = True) #matplotlib.rc('font', family='times') from matplotlib.colors import LogNorm from math import modf import matplotlib import numpy as np import pylab as plt from scipy.sparse.linalg import cg #from scipy.linalg.iterative import bicgstab from scipy.sparse.linalg import iterative as iter # import bicgstab #from iterative import bicgstab from matplotlib.ticker import FuncFormatter #from smooth import smooth from logbarrier import logBarrier from cmath import phase from decay import * class pwcTime: """Simple class for evaluating pwc extimation in the time domain""" def __init__(self): self.T2Limits = np.zeros(2) self.T2Bins = np.zeros(1) self.Larmor = 0. self.wL = 0. self.dt = 1e-4 # TODO this is critical and strange in FD, maybe error in mapping self.T = 3. #self.dead = 0.01327 # Wurks self.dead = 0.0 #15 #1e-3 #031415926 #.05# (np.pi*2)/(2000.) #0.01 #13 self.phi = 0. # phase self.t = np.arange(self.dead, self.T, self.dt) self.nt = len(self.t) self.G = np.zeros((2,2)) self.Gw = np.zeros((2,2)) def setPhi(self, phiIn): self.phi = phiIn def setT2(self, low, high, number, spacing): """ Sets the low and high limits of the T2 bins""" # Log spaced self.T2Limits[0] = low self.T2Limits[1] = high #self.T2Bins = np.arange(np.log(low), np.log(high), (np.log(high)-np.log(low))/number, endpoint=True) #self.T2Bins = np.linspace(np.log(low), np.log(high), (np.log(high)-np.log(low))/number, endpoint=True) #self.T2Bins = np.exp( self.T2Bins ) #print ("T2Bins", self.T2Bins) if spacing == "Log_10": self.T2Bins = np.logspace( np.log10(low), np.log10(high), num=number, endpoint=True, base=10 ) # base=10 elif spacing == "Log_2": self.T2Bins = np.logspace( np.log2(low), np.log2(high), num=number, endpoint=True, base=2 ) # base=10 elif spacing == "Log_e": self.T2Bins = np.logspace( np.log(low), np.log(high), num=number, endpoint=True, base=np.e ) # base=10 elif spacing == "Linear": self.T2Bins = np.linspace( low, high, num=number, endpoint=True ) # base=10 else: print("In setT2, spacing not recognized") exit(1) #print ("T2Bins", self.T2Bins) #exit() def setLarmor(self, freq): """ Set the Larmor frequency, in Hz""" self.Larmor = freq self.wL = np.pi * 2. * freq #def setSampling (self, dtin, Tin): # self.T = Tin # self.dt = dtin # self.t = np.arange(self.dead, self.T, self.dt) # self.nt = len(self.t) def setSampling (self, times): self.T = times[-1] self.dt = times[1] - times[0] self.t = times self.nt = len(self.t) def getSampling(self): return self.dt def getNumT2(self): """Returns the number of T2Bins""" return len(self.T2Bins) def getNumTime(self): """ Returns the number of times """ return self.nt def generateGenv(self): """ Fills in the matrix G, this is for spin-spin decay only """ #nt = (int) ( (self.T-self.dead) / self.dt) nT2 = len(self.T2Bins) # rows are times, cols are T2 bins self.Genv = np.zeros((self.nt, nT2)) # could do this w/o loops but who cares for this for iT2 in range(len(self.T2Bins)): for it in range(len(self.t)): self.Genv[it, iT2] = np.exp(-self.t[it]/self.T2Bins[iT2]) def generateJGKernel(self, TauE, Times, D_w, G , DTau): gamma_H = 26752. #/ (2*np.pi) # rad G^-1 s^-1 (wL = 2 pi gamma B0) TODO be aware of 2 pi term. #print( "TauE", TauE ) #print( "D_w", D_w ) #print( "G", G ) #print( "Times", Times ) ntt = 0 for time in TauE: ntt += len(Times[time]) self.JG = np.zeros( (ntt, len(G)) ) #print ("JG", self.JG) for gg in range(len(G)): it = 0 for TE in TauE: #print("TE", TE) NOTE in Dunn, this should be \tau == TE/2 tau = TE/2. #T2D = ( 3. * (gamma_H**2) * (G[gg]**2) * D_w * ((tau)**2) ) # Dunn #T2D = (D_w/3.) * tau * gamma_H**2 * G[gg]**2 # Keating 2009 T2D = (D_w * (tau**DTau) * (gamma_H**2) * (G[gg]**2) ) / 3. # Keating 2009 for t in Times[TE]: self.JG[it, gg] = np.exp( -t * T2D ) it += 1 #plt.matshow(self.JG) #plt.colorbar() #plt.show() def generateGDenv(self, TauE, Times, D, Gbins, Jg, DTau): """ Generates a diffusion sensitivity matrix TauE is an array of the TE echo times Times is a Dictionary of the Times so that Times[TauE[0]]...Times[TauE[-1]] are all 1D arrays of times D is an array of the D times. """ #G = 10. # G/cm ORS Tool #G = 6. # G/cm VC Mid field #G = .6 # G/cm VC Mid field ? ################################################## # Vista Clara with water # G = 100 Iteration 10 5.05285635664 # G = 60 Iteration 10 4.75446648427 # G = 10 Iteration 10 1.8 # G = 6 Iteration 10 1.1 # G = 3 Iteration 10 0.870649020301 #print( "TauE", TauE ) #print( "D", D ) #print( "Gbins", Gbins ) #print( "Times", Times ) #G = .8 # G/cm VC Mid field eta = 1 # unitless # TODO check units of gamma, and D. Are they both in si or cgs? gamma_H = 26752. # rad G^-1 s^-1 (wL = 2 pi gamma B0) TODO be aware of 2 pi term. #D_w = 1e-6 #2.3e-9 # diffusion of bulk water (Hurlimann 1995) diffusion limit # 2.34e-9 # [m^2/s] # Spyrou.pdf # 2.34e-5 # [cm^2/s] # Spyrou.pdf # # D_w = .140 m^2/s # Hayden --> .0000140 in cgs # So D *T_2 in cgs yields cm^2 -- which is the same as permeability ?? # D_0 # Hurliman1995 D_0 = 2.3e-9 [m^2/s] # Hayden 2004 D_0 = .140 [m^2/s] # extrude Genv --> len(Tau_E) times more data and len(D) times more model #self.GD = np.tile(self.Genv, (len(TauE)-1, len(D) )) ntt = 0 nT2 = len(self.T2Bins) #print (Times) #print (np.concatenate(Times)) for time in TauE: ntt += len(Times[time]) self.GD = np.zeros( (ntt , nT2*len(D)) ) #print (np.shape(self.GD)) """ [ {T2[0]D[0]} {T2[1]D[0]} ... {T2[nt]D[nd]} G = [ {T[0]TE[0]} [ {T[1]TE[0]} ... [ {T[0]TE[1]} [ {T[nt]TE[nTE]} """ for d in range(len(D)): for iT2 in range(len(self.T2Bins)): #for it in range(len(self.t)): #print ( "T2bins", self.T2Bins[iT2] ) it = 0 for TE in TauE: tau = TE#/2. for t in Times[TE]: ig = 0 for gg in Gbins: #print ( "Gradient info, gg=", gg, "J(g)", Jg[ig]) #T2D = (3* (gamma_H**2) * (gg**2) * D[d] * (tau**2) ) # (D[d]/12.) * ((gamma_H*G*(TE/2.))**2) # inverse tau = TE/2. #T2D = ( 3. * (gamma_H**2) * (G[gg]**2) * D_w * ((tau)**2) ) # Dunn #T2D = (D[d]/3.) * tau * gamma_H**2 * gg**2 # Keating 2009 T2D = (D[d] * tau**DTau * gamma_H**2 * gg**2) / 3. # Keating 2009 self.GD[it, d*nT2+iT2] += Jg[ig] * np.exp(-t/self.T2Bins[iT2]) * np.exp( -t * T2D ) ig += 1 it += 1 # now scale this bastard #TEs = np.ones( len(TauE) * ) #plt.figure() #plt.matshow(self.GD) #plt.colorbar() #plt.show() def plotGenv(self): # Make plot for SEG def tformatter(x, pos): return '%1.2f' %(x*self.dt) def t2formatter(x, pos): return '%0.3f' %(self.T2Bins[0] + x*(self.T2Bins[1]-self.T2Bins[0])) fig = plt.figure(997, facecolor='white') #self.G += 1. + 1e-3 plt.imshow(np.abs(self.Genv), aspect='auto', cmap='jet', norm=LogNorm(vmin=1e-4, vmax = 1.) ) cb = plt.colorbar() ax = plt.gca() for t in cb.ax.get_yticklabels(): t.set_color("black") t.set_fontsize(16) for t in cb.ax.get_yticklines(): t.set_color("black") cb.ax.spines['top'].set_color('black') cb.ax.spines['bottom'].set_color('black') cb.ax.spines['left'].set_color('black') cb.ax.spines['right'].set_color('black') cb.set_label(r"amplitude", color=('black')) formatter = FuncFormatter(tformatter) ax.yaxis.set_major_formatter(formatter) wformatter = FuncFormatter(t2formatter) ax.xaxis.set_major_formatter(wformatter) ax.set_xlabel(r"$T_2$ [s]", color = 'black') ax.set_ylabel(r"time [s]", color = 'black') ax.spines['top'].set_color('black') ax.spines['bottom'].set_color('black') ax.spines['left'].set_color('black') ax.spines['right'].set_color('black') for line in ax.yaxis.get_ticklines(): line.set_color('black') for line in ax.yaxis.get_ticklabels(): line.set_color('black') line.set_fontsize(16) for line in ax.xaxis.get_ticklines(): line.set_color('black') for line in ax.xaxis.get_minorticklines(): line.set_color('black') for line in ax.xaxis.get_ticklabels(): line.set_color('black') line.set_fontsize(16) plt.savefig('timeenvmatrix.png', dpi=200, facecolor='white', edgecolor='white') if __name__ == "__main__": Time = pwcTime() Time.setT2(.01, .6, 60) Time.setSampling(1e-3, 2) Time.generateGenv() # Generate some noisy data model = np.zeros(Time.getNumT2()) #model[5] = .05 model[16] = .1 model[17] = .15 model[18] = .15 model[19] = .1 #model[10] = .05 model[52] = .1 model[53] = .2 model[54] = .1 sigma = .02 * .4 noise = np.random.normal(0, sigma, Time.getNumTime()) dataenv = np.dot(Time.Genv, model) + noise cleanenv = np.dot(Time.Genv, model) #+ noise) # monoexponential fit [a0,b0,rt20] = regressCurve(dataenv, Time.t, intercept=True) # print(a0, b0, rt20) T2R = a0 + b0*np.exp(-np.array(Time.t)/rt20) guess = np.zeros( len(Time.T2Bins ) ) ib = 0 for bins in Time.T2Bins: if bins > rt20: guess[ib] = a0 + b0 break ib += 1 #exit() #timemodenv = logBarrier(Time.Genv, dataenv, x_0=guess, xr=guess, MAXITER=500, sigma=sigma, alpha=1e8) #, smooth=True) timemodenv = logBarrier(Time.Genv, dataenv, MAXITER=500, sigma=sigma, alpha=1e10, smooth=True) #, smooth=True) preenv = np.dot(Time.Genv, timemodenv) #+ noise) # Data plot fig = plt.figure(facecolor='white') axmine = fig.add_axes([.125,.1,.8,.8]) axmine.plot(Time.t, dataenv, color='red', alpha=.5, label='noisy data') axmine.plot(Time.t, cleanenv, color='blue', linewidth=4, label='true') axmine.plot(Time.t, preenv, color='green', linewidth=2, label='multi') axmine.plot(Time.t, T2R, color='black', linewidth=2, label='mono') #axmine.plot(Time.t, noise, color='purple', linewidth=1, label='noise') #axmine.set_xscale("log", nonposx='clip') plt.savefig("data2Dfits.pdf", facecolor='white', edgecolor='white', dpi=200) plt.legend() # Model plot plt.figure(1000) plt.plot(Time.T2Bins, model, linewidth=3, color='green', label='true') axmine = plt.gca() plt.semilogx(Time.T2Bins, timemodenv, color='purple', linewidth=2, label="multi-recovered") plt.semilogx(Time.T2Bins, guess, color='red', linewidth=2, label="mono-recovered") leg = plt.legend() axmine.set_xlabel(r"$T_2$ [s]", color='black') axmine.set_ylabel(r"$A_0$ [i]", color='black') plt.savefig("modelreconstruct.pdf", facecolor='white', edgecolor='white', dpi=200) plt.show()