modified: psava/Msmopick.py

user/psava/Msmopick.py

@@ -50,8 +50,6 @@ def lapop2D(nx,nz):
 n1 = Fin.int  ("n1")
 o1 = Fin.float("o1")
 d1 = Fin.float("d1")
-l1 = 0.5 * n1*d1              # xcorr time window
-z1 = int( abs(o1)/d1)
 
 nd = Fin.size(1);             # number of traces
 nm = nd
@@ -62,103 +60,89 @@ def lapop2D(nx,nz):
 Fou.put("o1",0.0)
 Fou.put('d1',1.0)
 
-Fou.put("n2",3) # temporary
-
+Fou.put("n2",3)
 Fou.put("n3",1)
 Fou.put("n4",1)
 
 # ------------------------------------------------------------
 # pick max value and weight
 # ------------------------------------------------------------
-din = np.zeros(n1,'f')
+din = np.zeros(n1,'f') # xcorrelation
 pck = np.zeros(nd,'f') # picks
 wgh = np.zeros(nd,'f') # weights
 
+# ------------------------------------------------------------
+l1 = 0.5 * n1*d1              # xcorr time window
+
 for i in range(nd):
     Fin.read(din)
 
-    # raw picks
-    #pck[i] = o1 + np.argmax(din) * d1
-    #wgh[i] = np.max(din)
-
-    din += 1.0 # avoid zero weights
-    pck[i] = o1 + z1 * d1
-    wgh[i] = din[z1]
-    for j in range(n1):
-        if  wgh[i] < din[j]:
-            pck[i] = o1 + j * d1
-            wgh[i] = din[j]
+    pck[i] = o1 + np.argmax(din) * d1  # raw picks
+    wgh[i] = np.max(din)               # data weight
 
-# bias weights
+# pick bias
 pckmed = np.median(pck)
-for i in range(nd):
-    #wgh[i]*= np.exp(- abs(pck[i]) )
-    wgh[i] *= np.cos( 0.5*np.pi * abs((pck[i]-pckmed)/l1) )
 
+# data weight
+for i in range(nd):
+    wgh[i] *= np.cos( 0.25*np.pi * abs((pck[i]-pckmed)/l1) )
 wgh /= np.max(wgh)
+wgh = 1 * np.power(wgh-np.min(wgh),wpo)
 
 # ------------------------------------------------------------
-# smooth picks
+# setup inverse problem
 # ------------------------------------------------------------
-
-xx = np.arange(nd)
-xS = xx[0]
-xE = xx[nd-1]
-
-nb = np.max([1,nd//10])
-#mbar = pck*0                          # reference model
-mbar = pck*0 + pckmed
+mbar = pck*0 + pckmed                 # reference model
 dbar = pck                            # rough picks
 
-wgh = np.power(wgh-np.min(wgh),wpo)
-
 WDop = spdiags(wgh,[0],nd,nd)         # data weight operator
 Gop  = idnop1D(nd)                    # mapping operator
-Rop  = lapop1D(nm)                    # Laplacian operator
+Rop  = lapop1D(nm)                    # regularization operator
 
 # ------------------------------------------------------------
 # L-curve
-ne = par.int('ne',1)
+# ------------------------------------------------------------
+ne = par.int  ('ne',1)
 oe = par.float('oe',0.0)
 de = par.float('de',+0.1)
 ee = np.arange(oe, oe+ne*de, de)
-#print(ee,file=sys.stderr)
 se = np.power(10,ee)
 
 ME = np.zeros( (nm,ne), dtype='float')
 rd = np.zeros(     ne,  dtype='float')
 rm = np.zeros(     ne,  dtype='float')
 
 for ie in range(ne):
-    WMop =  Rop / se[ie] # setup regularization
 
-    # solve IP
+    # scale the regularization operator
+    WMop =  Rop / se[ie] 
+
+    # solve the IP
     modE = spsolve( (WDop*Gop).T * (WDop*Gop)  + WMop.T * WMop , \
                     (WDop*Gop).T * WDop * dbar + WMop.T * WMop * mbar)
-    ME[:,ie] = modE # store model
+
+    # store the model
+    ME[:,ie] = modE 
 
     # compute residual norms
     rd[ie] = np.linalg.norm( WDop * (Gop * modE - dbar))
     rm[ie] = np.linalg.norm(  Rop * (      modE - mbar))
 
-if rd > 0:
-    rdn = rd / np.max(rd) # normalized residual norm
-else:
-    rdn = 1.0
-if rm > 0:
-    rmn = rm / np.max(rm)
-else:
-    rmn = 1.0
-rrn = rdn**2 + rmn**2 # distance from origin
+rdn = rd / np.max(rd) # normalized the  data residual norm
+rmn = rm / np.max(rm) # normalized the model residual norm
+rrn = rdn**2 + rmn**2 # compute the distance from origin
+
 je = np.argmin(rrn)   # index of the optimal model
-spk = ME[:,je]        # smooth picks
+spk = ME[:,je]        #              optimal model
+
+# print(je,ee[je], file=sys.stderr)
 
 # ------------------------------------------------------------
-# write picks
+# write picks and weights
 # ------------------------------------------------------------
 Fou.write(spk) # smooth picks
 Fou.write(pck) # rough picks
-Fou.write(wgh)
+Fou.write(wgh) # xcor pick weight
 
 Fin.close()
 Fou.close()