Stuff
authorT-bone <trevorirons@gmail.com>
Mon, 11 Dec 2017 15:00:49 +0000 (08:00 -0700)
committerT-bone <trevorirons@gmail.com>
Mon, 11 Dec 2017 15:00:49 +0000 (08:00 -0700)
Schillerslage/invertEAGE.py [new file with mode: 0644]
Schillerslage/invertTA.py
Schillerslage/plotKernel.py
Schillerslage/setup.sh

diff --git a/Schillerslage/invertEAGE.py b/Schillerslage/invertEAGE.py
new file mode 100644 (file)
index 0000000..db44a8b
--- /dev/null
@@ -0,0 +1,249 @@
+from SlidesPlot import *
+import numpy as np
+import sys
+import matplotlib.pyplot as plt
+from pylab import meshgrid 
+from logbarrier import *
+import yaml,os
+
+from matplotlib.colors import LogNorm
+from matplotlib.colors import LightSource
+from matplotlib.ticker import ScalarFormatter
+from matplotlib.ticker import MaxNLocator
+from matplotlib.ticker import AutoMinorLocator
+from matplotlib.ticker import LogLocator
+from matplotlib.ticker import FormatStrFormatter
+
+import cmocean
+
+class VectorXr(yaml.YAMLObject):
+    yaml_tag = r'VectorXr'
+    def __init__(self, array):
+        self.size = np.shape(array)[0]
+        self.data = array.tolist()
+    def __repr__(self):
+        #return np.array(self.data) # % (self.data)
+        return "%r" % (self.data)
+
+class MatrixXr(yaml.YAMLObject):
+    yaml_tag = u'MatrixXr'
+    def __init__(self, rows, cols, data):
+        self.rows = rows
+        self.cols = cols
+        self.data = np.zeros((rows,cols))
+
+class DataFID(yaml.YAMLObject):
+    yaml_tag = u'DataFID'
+    def __init__(self):
+        self.Serialized = "now"
+        self.Lemma_VERSION = "0.0"
+        self.Merlin_VERSION = "0.0"
+
+def buildKQT(K0,tg,T2Bins):
+    """ 
+        Constructs a QT inversion kernel from an initial amplitude one.  
+    """
+    nlay, nq = np.shape(K0)
+    nt2 = len(T2Bins)
+    nt = len(tg)
+    KQT = np.zeros( ( nq*nt,nt2*nlay) )
+    for iq in range(nq):
+        for it in range(nt):
+            for ilay in range(nlay):
+                for it2 in range(nt2):
+                    #KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
+                    #KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-((15+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
+                    KQT[iq*nt + it,ilay*nt2+it2] = K0[ilay,iq]*np.exp(-(tg[it]*1e-3)/(1e-3*T2Bins[it2])) # no RDP 
+    return KQT
+
+def loadAkvoData(fnamein, chan):
+    """ Loads data from an Akvo YAML file. The 0.02 is hard coded as the pulse length. This needs to be 
+        corrected in future kernel calculations. The current was reported but not the pulse length. 
+    """
+    fname = (os.path.splitext(fnamein)[0])
+    with open(fnamein, 'r') as stream:
+        try:
+            AKVO = (yaml.load(stream))
+        except yaml.YAMLError as exc:
+            print(exc)
+
+        #print( yaml.dump(AKVO) )
+
+        Z = 1e9*np.array(AKVO.FID_REAL.data).T
+        ZS = 1e9*np.array(AKVO.NoiseEstimate.data).T
+        T = np.array(AKVO.WindowCentres.data)
+
+        #QT = 1e9* np.array(AKVO.Qt_MOD.data)
+        #TruMod = np.array(AKVO.TrueModel.data)
+
+        #plt.matshow(QT)
+        #plt.colorbar()
+        #plt.title("c++")
+
+#     Z = np.zeros( (AKVO["nPulseMoments"], AKVO["GATED"]["Pulse 1"]["abscissa"].size ) )
+#     ZS = np.zeros( (AKVO["nPulseMoments"], AKVO["GATED"]["Pulse 1"]["abscissa"].size ) )
+#     for q in range(AKVO["nPulseMoments"]):
+#         Z[q] = AKVO["GATED"]["Pulse 1"][chan]["Q-"+str(q)].data
+#         if chan == "Chan. 1":
+#             ZS[q] = AKVO["GATED"]["Pulse 1"][chan]["Q-"+str(q) + " STD"].data
+#         elif chan == "Chan. 2":
+#             ZS[q] = AKVO["GATED"]["Pulse 1"][chan]["Q-"+str(q) + " STD"].data
+#         elif chan == "Chan. 4":
+#             ZS[q] = AKVO["GATED"]["Pulse 1"][chan]["Q-"+str(q) + " STD"].data
+#         else:
+#             print("DOOM!!!")
+#             exit()
+#     J = AKVO["Pulse 1"]["current"].data 
+#     J = np.append(J,J[-1]+(J[-1]-J[-2]))
+#     Q = AKVO["pulseLength"][0]*J
+
+    return Z, ZS, T#, QT, TruMod #, Q
+
+def loadK0(fname):
+    """ Loads in initial amplitude kernel
+    """
+    print("loading K0", fname)
+    ifaces = np.genfromtxt(fname, comments="#", max_rows=1 ) 
+    #q = 0.03*np.genfromtxt(fname, comments="#", skip_header=4, max_rows=1 ) # TODO pulse moment (0.02) hard coded! 
+    K1 = 1e9*np.genfromtxt(fname, comments="#", skip_header=6) 
+    nx, ny = np.shape(K1)
+    K1 = K1[0:nx//2] + 1j*K1[nx//2:]
+    #print(K1) 
+    return ifaces,np.abs(K1)
+
+if __name__ == "__main__":
+
+    if (len (sys.argv) < 3):
+        print ("python invertTA.py  K0.dat  Data.yaml")
+    print (sys.argv[1])
+    with open(sys.argv[1], 'r') as stream:
+        try:
+            cont = (yaml.load(stream))
+        except yaml.YAMLError as exc:
+            print(exc)
+    ###############################################
+    # Load in data
+    ###############################################
+    V = []
+    VS = []
+    tg = 0
+    for dat in cont['data']:
+        for ch in dat['channels']:
+            v,vs,tg = loadAkvoData(dat['file'], ch)
+            V.append(v)
+            VS.append(vs)
+    for iv in range(1, len(V)):
+        V[0] = np.concatenate( (V[0], V[iv]) )
+        VS[0] = np.concatenate( (VS[0], VS[iv]) )
+    V = V[0]
+    VS = VS[0]
+
+    #plt.matshow(V, cmap='RdBu', vmin=-np.max(np.abs(V)), vmax=np.max(np.abs(V)))
+    #plt.title("data")
+    #plt.colorbar()
+    
+    #plt.matshow(VS, cmap='inferno_r', vmin=0, vmax=np.max(np.abs(VS)))
+    #plt.title("error")
+    #plt.colorbar()
+    
+#     trudat = np.dot(qtt, tru)
+#     TRU = np.reshape( trudat, np.shape(V)  )
+#     plt.matshow(TRU, cmap='inferno')
+#     plt.gca().set_title("trudatt")
+#     plt.colorbar()
+#     plt.show()
+#     exit()
+
+    ###############################################
+    # Load in kernels
+    ###############################################
+    K0 = []
+    for kern in cont["K0"]:
+        ifaces,k0 = loadK0( kern )
+        K0.append(k0)
+    for ik in range(1, len(K0)):
+        K0[0] = np.concatenate( (K0[0].T, K0[ik].T) ).T
+    K0 = K0[0]
+    #plt.matshow(K0)
+
+    ###############################################
+    # Build full kernel
+    ############################################### 
+    T2Bins = np.logspace( np.log10(cont["T2Bins"]["low"]), np.log10(cont["T2Bins"]["high"]), cont["T2Bins"]["number"], endpoint=True, base=10)  
+    KQT = buildKQT(K0,tg,T2Bins)
+    #print("Constructed KQT")
+    #plt.matshow(KQT)
+    #plt.colorbar()
+    #plt.title("py")
+    
+
+    # Chan. 1    reg = 1.5   
+    # Chan. 2    reg = 1.35
+    # Chan. 4    reg = 1.95  
+    ###############################################
+    # Invert
+    ############################################### 
+    inv = logBarrier(KQT, np.ravel(V, order='C'), MAXITER=50, sigma=np.ravel(2.25*VS, order='C'), alpha=1e15, smooth=False, nT2=len(T2Bins) ) #, smooth=True) # 
+
+    ###############################################
+    # Appraise
+    ############################################### 
+    pre = np.dot(KQT,inv) 
+    PRE = np.reshape( pre, np.shape(V)  )
+    plt.matshow(PRE, cmap='Blues')
+    plt.colorbar()
+
+    DIFF = (PRE-V) #/ (1.5*VS)
+    plt.matshow(DIFF, cmap='inferno')
+    plt.gca().set_title("misfit / $\widehat{\sigma}$")
+    plt.colorbar()
+    
+    DIFF = (V) #/ (1.5*VS)
+    plt.matshow(DIFF, cmap='inferno')
+    plt.gca().set_title("V")
+    plt.colorbar()
+
+
+    T2Bins = np.append( T2Bins, T2Bins[-1] + (T2Bins[-1]-T2Bins[-2]) )
+    
+    print( np.shape(inv), len(ifaces)-1,cont["T2Bins"]["number"] )
+    INV = np.reshape(inv, (len(ifaces)-1,cont["T2Bins"]["number"]), order='C' )
+    #INV = np.reshape(tru, (len(ifaces)-1, cont["T2Bins"]["number"]) ) # TRUE MODEL
+    #INV = np.reshape(tru, (len(ifaces)-1, cont["T2Bins"]["number"]), order='F' ) # TRUE MODEL
+    Y,X = meshgrid( ifaces, T2Bins )
+    fig = plt.figure( figsize=(pc2in(17.0),pc2in(18.)) )
+    ax1 = fig.add_axes( [.2,.15,.6,.7] )
+    im = ax1.pcolor(X, Y, INV.T, cmap=cmocean.cm.tempo)#'viridis')
+    im.set_edgecolor('face')
+    ax1.set_xlim( T2Bins[0], T2Bins[-1] )
+    ax1.set_ylim( ifaces[-1], ifaces[0] )
+    cb = plt.colorbar(im, label = u"PWC (m$^3$/m$^3$)") #, format='%1.1f')
+    cb.locator = MaxNLocator( nbins = 4)
+    cb.ax.yaxis.set_offset_position('left')                         
+    cb.update_ticks()
+    ax1.set_xlabel(u"$T_2^*$ (ms)")
+    ax1.set_ylabel(u"depth (m)$")
+    
+    ax1.get_xaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+    ax1.get_yaxis().set_major_formatter(FormatStrFormatter('%1.0f'))
+    ax1.xaxis.set_major_locator( MaxNLocator(nbins = 4) )   
+
+    ax1.set_xscale('log')
+    #ax1.xaxis.set_label_position('top') 
+
+    ax2 = ax1.twiny()
+    ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='black' )
+    ax2.set_xlabel(u"total water (m$^3$/m$^3$)")
+    ax2.set_ylim( ifaces[-1], ifaces[0] )
+    ax2.xaxis.set_major_locator( MaxNLocator(nbins = 3) )   
+    ax2.get_xaxis().set_major_formatter(FormatStrFormatter('%0.2f'))
+    #ax2.xaxis.set_label_position('bottom') 
+
+    plt.savefig(sys.argv[2]+".pdf")
+    
+    plt.show()
+
index 5f837cd..5861bca 100644 (file)
@@ -70,12 +70,15 @@ def loadK0(fname):
     """ Loads in initial amplitude kernel
     """
     print("loading K0", fname)
-    ifaces = np.genfromtxt(fname, comments="#", max_rows=1 ) 
-    q = 0.03*np.genfromtxt(fname, comments="#", skip_header=6, max_rows=1 ) # TODO pulse moment (0.02) hard coded! 
-    K1 = 1e9*np.genfromtxt(fname, comments="#", skip_header=7) 
-    nx, ny = np.shape(K1)
-    K1 = K1[0:nx//2] + 1j*K1[nx//2:] 
-    return ifaces,np.abs(K1)
+#    ifaces = np.genfromtxt(fname, comments="#", max_rows=1 ) 
+#    q = 0.03*np.genfromtxt(fname, comments="#", skip_header=6, max_rows=1 ) # TODO pulse moment (0.02) hard coded! 
+#    K1 = 1e9*np.genfromtxt(fname, comments="#", skip_header=7) 
+#    nx, ny = np.shape(K1)
+#    K1 = K1[0:nx//2] + 1j*K1[nx//2:] 
+#    return ifaces,np.abs(K1)
+    K0 = yaml.load(fname)
+    ifaces = K0.Interfaces.data
+
 
 if __name__ == "__main__":
 
index 44f0734..ea989da 100644 (file)
@@ -24,8 +24,6 @@ if __name__ == "__main__":
         # use safe_load instead load
         K0 = yaml.load(f)
 
-    
-
     K = catLayers(K0.K0)
     print(K)
     plt.pcolor(np.abs(K))
index fd2607e..65602bd 100755 (executable)
@@ -4,10 +4,10 @@
 ./CircularLoop 10 0 0 5 21 > Ch3.yaml
 
 # create dummy kernel 
-./KernelAligner 3loop.yaml
+#./KernelAligner 3loop.yaml
 
 # calculate kernels 
-./KernelV0-2 Kern.yaml  Ch1  100m
+./KernelV0-2 Kern.yaml  Ch1  Ch1 false  
 
 # align kernel with data and set T2* bins
-./ModelAligner 
+#./ModelAligner