updates for inversionwq
authorTrevor Irons <tirons@egi.utah.edu>
Fri, 8 Dec 2017 20:53:36 +0000 (13:53 -0700)
committerTrevor Irons <tirons@egi.utah.edu>
Fri, 8 Dec 2017 20:53:36 +0000 (13:53 -0700)
Schillerslage/Ch1.yaml
Schillerslage/SlidesPlot.py [new file with mode: 0644]
Schillerslage/akvodata.yaml
Schillerslage/invert.yaml [new file with mode: 0644]
Schillerslage/invertTA.py [new file with mode: 0644]
Schillerslage/lemma_yaml.py
Schillerslage/logbarrier.py [new file with mode: 0644]
Schillerslage/plotKernel.py

index dfca20c..26be07a 100644 (file)
@@ -97,5 +97,4 @@ Freqs: !<VectorXr>
 minDipoleRatio: 0.15
 maxDipoleMoment: 10
 minDipoleMoment: 1e-06
----
 
diff --git a/Schillerslage/SlidesPlot.py b/Schillerslage/SlidesPlot.py
new file mode 100644 (file)
index 0000000..bf89c5c
--- /dev/null
@@ -0,0 +1,34 @@
+#################################################################################
+# GJI final pub specs                                                           #
+import matplotlib                                                               #
+from matplotlib import rc                                                           #
+matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{timet,amsmath,amssymb}"]  #
+rc('font',**{'family':'sans-serif','serif':['timet']})                                   #
+rc('font',**{'size':11})                                                             #
+rc('text', usetex=True)                                                             #
+# converts pc that GJI is defined in to inches                                      # 
+# In GEOPHYSICS \textwidth = 42pc                                               #
+#        \columnwidth = 20pc                                                    #
+#        one column widthe figures are 20 picas                                 #
+#        one and one third column figures are 26 picas                          #
+def pc2in(pc):                                                                  #
+    return pc*12/72.27                                                          #
+#################################################################################
+import numpy as np
+light_grey = np.array([float(248)/float(255)]*3)
+
+def fixLeg(legend):
+    rect = legend.get_frame()
+    #rect.set_color('None')
+    rect.set_facecolor(light_grey)
+    rect.set_linewidth(0.0)
+    rect.set_alpha(0.5)
+
+def deSpine(ax1):
+    spines_to_remove = ['top', 'right']
+    for spine in spines_to_remove:
+        ax1.spines[spine].set_visible(False)
+    #ax1.xaxis.set_ticks_position('none')
+    #ax1.yaxis.set_ticks_position('none')
+    ax1.get_xaxis().tick_bottom()
+    ax1.get_yaxis().tick_left()
index 863fde2..04cb49b 100644 (file)
@@ -1,4 +1,4 @@
-!AkvoData
+!<AkvoData>
 Akvo_VERSION: 1.0.5
 Import:
   GMR Header: /home/tirons/Data/Liege/All multi-central/central_loop
diff --git a/Schillerslage/invert.yaml b/Schillerslage/invert.yaml
new file mode 100644 (file)
index 0000000..1fb7c85
--- /dev/null
@@ -0,0 +1,11 @@
+data:
+  - file : akvodata.yaml 
+    channels: 
+      - Chan. 1
+# These are the initial amplitude kernels 
+K0:
+  - "TxCh1Rx-Ch1.dat" 
+T2Bins:
+  low : 10 
+  high : 500 
+  number : 20
diff --git a/Schillerslage/invertTA.py b/Schillerslage/invertTA.py
new file mode 100644 (file)
index 0000000..5f837cd
--- /dev/null
@@ -0,0 +1,193 @@
+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
+
+from lemma_yaml import * 
+
+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(-((10+tg[it])*1e-3)/(1e-3*T2Bins[it2]))
+    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)
+            exit()
+
+    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) + " CA"].data
+        if chan == "Chan. 1":
+            ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
+            ZS[q] *= 1.5
+        elif chan == "Chan. 2":
+            ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
+            ZS[q] *= 1.35
+        elif chan == "Chan. 4":
+            ZS[q] = AKVO.Gated["Pulse 1"][chan]["STD"].data
+            ZS[q] *= 1.95
+        else:
+            print("DOOM!!!")
+            exit()
+    J = AKVO.Pulses["Pulse 1"]["current"].data 
+    J = np.append(J,J[-1]+(J[-1]-J[-2]))
+    Q = AKVO.pulseLength[0]*J
+
+    return Z, ZS, AKVO.Gated["Pulse 1"]["abscissa"].data  #, 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=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)
+
+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', vmin=-np.max(np.abs(VS)), vmax=np.max(np.abs(VS)))
+    #plt.title("error")
+    #plt.colorbar()
+
+    ###############################################
+    # 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)
+    #plt.matshow(KQT)
+
+    # Chan. 1    reg = 1.5   
+    # Chan. 2    reg = 1.35
+    # Chan. 4    reg = 1.95  
+    ###############################################
+    # Invert
+    ############################################### 
+    inv = logBarrier(KQT, np.ravel(V), MAXITER=50, sigma=np.ravel(0.6*VS), alpha=1e15, smooth=True, 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()
+
+    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"]) )
+    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='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.xaxis.set_label_position('top') 
+
+    ax2 = ax1.twiny()
+    ax2.plot( np.sum(INV, axis=1), (ifaces[1:]+ifaces[0:-1])/2 ,  color='red' )
+    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("inversion-chtempa.pdf")
+    
+    plt.show()
+
index b2d4ba7..bdb6e31 100644 (file)
@@ -109,3 +109,13 @@ class PolygonalWireAntenna(yaml.YAMLObject):
     def __init__(self, val):
         self.val = val
 
+class AkvoData(yaml.YAMLObject):
+    yaml_tag = u'AkvoData'
+    def __init__(self, obj): #nPulseMoments, pulseLength):
+    #def __init__(self, rows, cols, data):
+        #self.nPulseMoments = nPulseMoments
+        #self.pulseLength = pulseLength 
+        #for key in obj.keys:
+        #    self[key] = obj.key
+        pass
+
diff --git a/Schillerslage/logbarrier.py b/Schillerslage/logbarrier.py
new file mode 100644 (file)
index 0000000..fe30286
--- /dev/null
@@ -0,0 +1,223 @@
+from __future__ import division
+import numpy as np
+from scipy.sparse.linalg import iterative  as iter
+import pylab 
+import pprint 
+from scipy.optimize import nnls 
+
+def PhiB(mux, muy, minVal, maxVal, x):
+    phib = mux * np.abs( np.sum(np.log( x-minVal)) )
+    return phib
+
+def curvature(x, y):
+    from scipy.ndimage import gaussian_filter1d
+    #first and second derivative
+    x1 = gaussian_filter1d(x, sigma=1, order=1)#, mode='constant', cval=x[-1])
+    x2 = gaussian_filter1d(x1, sigma=1, order=1)#, mode='constant', cval=y[-1])
+    y1 = gaussian_filter1d(y, sigma=1, order=1)#, mode='constant', cval=x1[-1])
+    y2 = gaussian_filter1d(y1, sigma=1, order=1)#, mode='constant', cval=y1[-1])
+    return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
+
+def logBarrier(A, b, T2Bins, x_0=0, xr=0, alpha=10, mu1=10, mu2=10, smooth=False, MAXITER=70, fignum=1000, sigma=1, callback=None):
+    """Impliments a log barrier Tikhonov solution to a linear system of equations 
+        Ax = b  s.t.  x_min < x < x_max. A log-barrier term is used for the constraint
+    """
+    # TODO input
+    minVal = 0.0
+    maxVal = 1e8
+
+    Wd =    (np.eye(len(b)) / (sigma))        # Wd = eye( sigma )
+    WdTWd = (np.eye(len(b)) / (sigma**2))     # Wd = eye( sigma )
+
+    ATWdTWdA = np.dot(A.conj().transpose(), np.dot( WdTWd, A ))     # TODO, implicit calculation instead?
+    N = np.shape(A)[1]                        # number of model
+    M = np.shape(A)[0]                        # number of data
+    SIGMA = .25 
+    EPSILON = 1e-35 
+
+    # reference model
+    if np.size(xr) == 1:
+        xr =  np.zeros(N)     
+    
+    # initial guess
+    if np.size(x_0) == 1:
+        x = 1e-10 + np.zeros(N)     
+    else:
+        x = 1e-10 + x_0
+        
+    # Construct model constraint base   
+    Phim_base = np.zeros( [N , N] ) 
+    a1 = 1.0     # smallest
+    
+    # calculate largest term            
+    D1 = 1./abs(T2Bins[1]-T2Bins[0])
+    D2 = 1./abs(T2Bins[2]-T2Bins[1])
+    #a2 = 1. #(1./(2.*D1+D2))    # smooth
+    
+    if smooth == "Both":
+        #print ("Both small and smooth model")
+        for ip in range(N):
+            D1 = 0.
+            D2 = 0.
+            if ip > 0:
+                #D1 = np.sqrt(1./abs(T2Bins[ip]-T2Bins[ip-1]))**.5
+                D1 = (1./abs(T2Bins[ip]-T2Bins[ip-1])) #**2
+            if ip < N-1:
+                #D2 = np.sqrt(1./abs(T2Bins[ip+1]-T2Bins[ip]))**.5
+                D2 = (1./abs(T2Bins[ip+1]-T2Bins[ip])) #**2
+            if ip > 0:
+                Phim_base[ip,ip-1] =   -(D1)      
+            if ip == 0:
+                Phim_base[ip,ip  ] = 1.*(D1+D2)  
+            elif ip == N-1:
+                Phim_base[ip,ip  ] = 2.*(D1+D2) 
+            else:
+                Phim_base[ip,ip  ] = 2.*(D1+D2)
+            if ip < N-1:
+                Phim_base[ip,ip+1] =   -(D2)  
+        Phim_base /= np.max(Phim_base)            # normalize 
+        Phim_base += a1*np.eye(N)
+
+    elif smooth == "Smooth":
+        #print ("Smooth model")
+        for ip in range(N):
+            if ip > 0:
+                Phim_base[ip,ip-1] = -1    # smooth in log space
+            if ip == 0:
+                Phim_base[ip,ip  ] = 2.5   # Encourage a little low model
+            elif ip == N-1:
+                Phim_base[ip,ip  ] = 2.5   # Penalize long decays
+            else:
+                Phim_base[ip,ip  ] = 2.5   # Smooth and small
+            if ip < N-1:
+                Phim_base[ip,ip+1] = -1    # smooth in log space
+
+    else: 
+        #print ("SMALLEST model norm", smooth)
+        # Smallest model
+        for ip in range(N):
+            Phim_base[ip,ip  ] = 1.
+    
+    Phi_m =  alpha*Phim_base
+    WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)            
+    b_pre = np.dot(A, x)
+    phid = np.linalg.norm( np.dot(Wd, (b-b_pre)) )**2
+    phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
+
+    mu2 = phim
+    phib = PhiB(mu1, mu2, 0, 1e8, x) 
+    mu1 = ((phid + alpha*phim) / phib)
+
+    PHIM = []
+    PHID = []
+    MOD = []
+
+    for i in range(MAXITER):
+        
+        Phi_m =  alpha*Phim_base
+        
+        # reset mu1 at each iteration 
+        # Calvetti -> No ; Li -> Yes   
+        # without this, non monotonic convergence occurs...which is OK if you really trust your noise 
+        mu1 = ((phid + alpha*phim) / phib) 
+
+        WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)            
+        phid_old = phid
+        inner = 0
+
+        First = True # guarantee entry 
+
+        xp = np.copy(x) # prior step x 
+
+        while ( (phib / (phid+alpha*phim)) > EPSILON  or First==True ):
+
+            First = False
+            # Log barrier, keep each element above minVal
+            X1 = np.eye(N) * (x-minVal)**-1           
+            X2 = np.eye(N) * (x-minVal)**-2         
+            
+            # Log barrier, keep sum below maxVal TODO normalize by component. Don't want to push all down  
+            Y1 = np.eye(N) * (maxVal - np.sum(x))**-1           
+            Y2 = np.eye(N) * (maxVal - np.sum(x))**-2         
+            
+            AA = ATWdTWdA + mu1*X2 + mu2*Y2 + Phi_m 
+            M = np.eye( N ) * (1./np.diag(ATWdTWdA + mu1*X2 + mu2*Y2 + Phi_m))
+        
+            # Solve system (newton step) (Li)
+            b2 = np.dot(A.conj().transpose(), np.dot(WdTWd, b-b_pre) ) + 2.*mu1*np.diag(X1) + 2.*mu2*np.diag(Y1) - alpha*np.dot(WmTWm,(x-xr))
+            ztilde = iter.cg(AA, b2, M = M) 
+            h = (ztilde[0].real) 
+            
+            # Solve system (direct solution) (Calvetti) 
+            #b2 = np.dot(A.conj().transpose(), np.dot(WdTWd, b)) + 2.*mu1*np.diag(X1) + 2.*mu2*np.diag(Y1) - alpha*np.dot(WmTWm,(x-xr))
+            #ztilde = iter.cg(AA, b2, M=M, x0=x) 
+            #h = (ztilde[0].real - x) 
+
+            # step size
+            d = np.min( (1, 0.95 * np.min(x/np.abs(h+1e-120))) )
+            
+            ##########################################################
+            # Update and fix any over/under stepping
+            x += d*h
+        
+            # Determine mu steps to take
+            s1 = mu1 * (np.dot(X2, ztilde[0].real) - 2.*np.diag(X1))
+            s2 = mu2 * (np.dot(Y2, ztilde[0].real) - 2.*np.diag(Y1))
+
+            # determine mu for next step
+            mu1 = SIGMA/N * np.abs(np.dot(s1, x))
+            mu2 = SIGMA/N * np.abs(np.dot(s2, x))
+            
+            b_pre = np.dot(A, x)
+            phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+            phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
+            phib = PhiB(mu1, mu2, minVal, maxVal, x)
+            inner += 1
+        
+        PHIM.append(phim)      
+        PHID.append(phid)      
+        MOD.append(np.copy(x))  
+
+        # determine alpha
+        scale = (len(b)/phid)
+        #alpha *= np.sqrt(scale)
+        alpha *= min(scale, .95)
+
+#         if np.sqrt(phid/len(b)) < 0.97: 
+#             ibreak = -1
+#             print ("------------overshot--------------------", alpha, np.sqrt(phid/len(b)), ibreak)
+#             alpha *= 2. #0
+#             x -= d*h
+#             b_pre = np.dot(A, x)
+#             phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+#             phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
+#             mu1 = ((phid + alpha*phim) / phib) 
+        #if np.sqrt(phid/len(b)) < 1.00: 
+        #    ibreak = 0
+        #    print ("optimal solution found", alpha, np.sqrt(phid/len(b)), ibreak)
+        #    break
+        # slow convergence, bail and use L-curve 
+        # TI- only use L-curve. Otherwise results for perlin noise are too spurious for paper.  
+        if np.sqrt(phid/len(b)) < 2. and ( (np.sqrt(phid_old/len(b))-np.sqrt(phid/len(b))) < 1e-4): 
+            ibreak = 1
+            MOD = np.array(MOD)
+            print ("###########################") #slow convergence", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
+            print ("Using L-curve criteria") 
+            kappa = curvature(np.log(np.array(PHIM)), np.log(np.array(PHID)))
+            x = MOD[ np.argmax(kappa) ]
+            b_pre = np.dot(A, x)
+            phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
+            phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
+            mu1 = ((phid + alpha*phim) / phib) 
+            print ("L-curve selected", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
+            print ("###########################")
+            break
+
+    if (i == MAXITER-1 ):
+        ibreak = 2
+        print("Reached max iterations!!", alpha, np.sqrt(phid/len(b)), ibreak)
+
+    return x, ibreak, np.sqrt(phid/len(b))
+
+
+
index 599cd5a..44f0734 100644 (file)
@@ -12,7 +12,7 @@ import yaml
 from  lemma_yaml import *
 
 def catLayers(K0):
-    #print( type(K0))
+    print( type(K0))
     K = np.zeros( (len(K0.keys()), len(K0["layer-0"].data)) , dtype=complex )
     for lay in range(len(K0.keys())):
         #print(K0["layer-"+str(lay)].data) #    print (lay)
@@ -24,6 +24,8 @@ if __name__ == "__main__":
         # use safe_load instead load
         K0 = yaml.load(f)
 
+    
+
     K = catLayers(K0.K0)
     print(K)
     plt.pcolor(np.abs(K))