AkvoFluo/pwctimeWhite.py

#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()