Borehole NMR processing
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

logbarrier.py 10KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267
  1. from __future__ import division
  2. import numpy as np
  3. from scipy.sparse.linalg import iterative as iter
  4. import pylab
  5. import pprint
  6. from scipy.optimize import nnls
  7. import matplotlib.pyplot as plt
  8. def PhiB(mux, muy, minVal, maxVal, x):
  9. phib = mux * np.abs( np.sum(np.log( x-minVal)) )
  10. return phib
  11. def curvaturefd(x, y, t):
  12. x1 = np.gradient(x,t)
  13. x2 = np.gradient(x1,t)
  14. y1 = np.gradient(y,t)
  15. y2 = np.gradient(y1,t)
  16. return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
  17. def curvatureg(x, y):
  18. from scipy.ndimage import gaussian_filter1d
  19. #first and second derivative
  20. x1 = gaussian_filter1d(x, sigma=1, order=1)#, mode='constant', cval=x[-1])
  21. x2 = gaussian_filter1d(x1, sigma=1, order=1)#, mode='constant', cval=y[-1])
  22. y1 = gaussian_filter1d(y, sigma=1, order=1)#, mode='constant', cval=x1[-1])
  23. y2 = gaussian_filter1d(y1, sigma=1, order=1)#, mode='constant', cval=y1[-1])
  24. return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
  25. def logBarrier(A, b, T2Bins, lambdastar, x_0=0, xr=0, alpha=10, mu1=10, mu2=10, smooth=False, MAXITER=70, fignum=1000, sigma=1, callback=None):
  26. """Impliments a log barrier Tikhonov solution to a linear system of equations
  27. Ax = b s.t. x_min < x < x_max. A log-barrier term is used for the constraint
  28. """
  29. # TODO input
  30. minVal = 0.0
  31. maxVal = 1e8
  32. Wd = (np.eye(len(b)) / (sigma)) # Wd = eye( sigma )
  33. WdTWd = (np.eye(len(b)) / (sigma**2)) # Wd = eye( sigma )
  34. ATWdTWdA = np.dot(A.conj().transpose(), np.dot( WdTWd, A )) # TODO, implicit calculation instead?
  35. N = np.shape(A)[1] # number of model
  36. M = np.shape(A)[0] # number of data
  37. SIGMA = .25
  38. EPSILON = 1e-35
  39. # reference model
  40. if np.size(xr) == 1:
  41. xr = np.zeros(N)
  42. # initial guess
  43. if np.size(x_0) == 1:
  44. x = 1e-10 + np.zeros(N)
  45. else:
  46. x = 1e-10 + x_0
  47. # Construct model constraint base
  48. Phim_base = np.zeros( [N , N] )
  49. a1 = .05 # smallest too
  50. # calculate largest term
  51. D1 = 1./abs(T2Bins[1]-T2Bins[0])
  52. D2 = 1./abs(T2Bins[2]-T2Bins[1])
  53. #a2 = 1. #(1./(2.*D1+D2)) # smooth
  54. if smooth == "Both":
  55. #print ("Both small and smooth model")
  56. for ip in range(N):
  57. D1 = 0.
  58. D2 = 0.
  59. if ip > 0:
  60. #D1 = np.sqrt(1./abs(T2Bins[ip]-T2Bins[ip-1]))**.5
  61. D1 = (1./abs(T2Bins[ip]-T2Bins[ip-1])) #**2
  62. if ip < N-1:
  63. #D2 = np.sqrt(1./abs(T2Bins[ip+1]-T2Bins[ip]))**.5
  64. D2 = (1./abs(T2Bins[ip+1]-T2Bins[ip])) #**2
  65. if ip > 0:
  66. Phim_base[ip,ip-1] = -(D1)
  67. if ip == 0:
  68. Phim_base[ip,ip ] = 2.*(D1+D2)
  69. elif ip == N-1:
  70. Phim_base[ip,ip ] = 2.*(D1+D2)
  71. else:
  72. Phim_base[ip,ip ] = 2.*(D1+D2)
  73. if ip < N-1:
  74. Phim_base[ip,ip+1] = -(D2)
  75. Phim_base /= np.max(Phim_base) # normalize
  76. Phim_base += a1*np.eye(N)
  77. elif smooth == "Smooth":
  78. #print ("Smooth model")
  79. for ip in range(N):
  80. if ip > 0:
  81. Phim_base[ip,ip-1] = -1 # smooth in log space
  82. if ip == 0:
  83. Phim_base[ip,ip ] = 2.05 # Encourage a little low model
  84. elif ip == N-1:
  85. Phim_base[ip,ip ] = 2.5 # Penalize long decays
  86. else:
  87. Phim_base[ip,ip ] = 2.1 # Smooth and small
  88. if ip < N-1:
  89. Phim_base[ip,ip+1] = -1 # smooth in log space
  90. elif smooth == "Smallest":
  91. for ip in range(N):
  92. Phim_base[ip,ip ] = 1.
  93. else:
  94. print("non valid model constraint:", smooth)
  95. exit()
  96. Phi_m = alpha*Phim_base
  97. WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)
  98. b_pre = np.dot(A, x)
  99. phid = np.linalg.norm( np.dot(Wd, (b-b_pre)) )**2
  100. phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
  101. mu2 = phim
  102. phib = PhiB(mu1, mu2, 0, 1e8, x)
  103. mu1 = ((phid + alpha*phim) / phib)
  104. PHIM = []
  105. PHID = []
  106. MOD = []
  107. ALPHA = []
  108. ALPHA.append(alpha)
  109. for i in range(MAXITER):
  110. Phi_m = alpha*Phim_base
  111. # reset mu1 at each iteration
  112. # Calvetti -> No ; Li -> Yes
  113. # without this, non monotonic convergence occurs...which is OK if you really trust your noise
  114. mu1 = ((phid + alpha*phim) / phib)
  115. WmTWm = Phim_base # np.dot(Phim_base, Phim_base.T)
  116. phid_old = phid
  117. inner = 0
  118. First = True # guarantee entry
  119. xp = np.copy(x) # prior step x
  120. while ( (phib / (phid+alpha*phim)) > EPSILON or First==True ):
  121. First = False
  122. # Log barrier, keep each element above minVal
  123. X1 = np.eye(N) * (x-minVal)**-1
  124. X2 = np.eye(N) * (x-minVal)**-2
  125. # Log barrier, keep sum below maxVal TODO normalize by component. Don't want to push all down
  126. Y1 = np.eye(N) * (maxVal - np.sum(x))**-1
  127. Y2 = np.eye(N) * (maxVal - np.sum(x))**-2
  128. AA = ATWdTWdA + mu1*X2 + mu2*Y2 + Phi_m
  129. M = np.eye( N ) * (1./np.diag(ATWdTWdA + mu1*X2 + mu2*Y2 + Phi_m))
  130. # Solve system (newton step) (Li)
  131. 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))
  132. ztilde = iter.cg(AA, b2, M = M)
  133. h = (ztilde[0].real)
  134. # Solve system (direct solution) (Calvetti)
  135. #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))
  136. #ztilde = iter.cg(AA, b2, M=M, x0=x)
  137. #h = (ztilde[0].real - x)
  138. # step size
  139. d = np.min( (1, 0.95 * np.min(x/np.abs(h+1e-120))) )
  140. ##########################################################
  141. # Update and fix any over/under stepping
  142. x += d*h
  143. # Determine mu steps to take
  144. s1 = mu1 * (np.dot(X2, ztilde[0].real) - 2.*np.diag(X1))
  145. s2 = mu2 * (np.dot(Y2, ztilde[0].real) - 2.*np.diag(Y1))
  146. # determine mu for next step
  147. mu1 = SIGMA/N * np.abs(np.dot(s1, x))
  148. mu2 = SIGMA/N * np.abs(np.dot(s2, x))
  149. b_pre = np.dot(A, x)
  150. phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
  151. phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
  152. phib = PhiB(mu1, mu2, minVal, maxVal, x)
  153. inner += 1
  154. PHIM.append(phim)
  155. PHID.append(phid)
  156. MOD.append(np.copy(x))
  157. # determine alpha
  158. scale = (len(b)/phid)
  159. #alpha *= np.sqrt(scale)
  160. alpha *= min(scale, .95) # was .85...
  161. ALPHA.append(alpha)
  162. #alpha = ALPHA[i+1]
  163. # if np.sqrt(phid/len(b)) < 0.97:
  164. # ibreak = -1
  165. # print ("------------overshot--------------------", alpha, np.sqrt(phid/len(b)), ibreak)
  166. # alpha *= 2. #0
  167. # x -= d*h
  168. # b_pre = np.dot(A, x)
  169. # phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
  170. # phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )#**2
  171. # mu1 = ((phid + alpha*phim) / phib)
  172. if lambdastar == "discrepency":
  173. if np.sqrt(phid/len(b)) < 1.00 or alpha < 1e-5:
  174. ibreak = 1
  175. print ("optimal solution found", alpha, np.sqrt(phid/len(b)), ibreak)
  176. break
  177. # slow convergence, bail and use L-curve
  178. # TI- only use L-curve. Otherwise results for perlin noise are too spurious for paper.
  179. if lambdastar == "lcurve":
  180. if i>20 and ((np.sqrt(phid_old/len(b))-np.sqrt(phid/len(b))) < 1e-4):
  181. #if np.sqrt(phid/len(b)) < 3.0 and ((np.sqrt(phid_old/len(b))-np.sqrt(phid/len(b))) < 1e-3):
  182. ibreak = 1
  183. MOD = np.array(MOD)
  184. print ("###########################") #slow convergence", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
  185. print ("Using L-curve criteria")
  186. kappa = curvaturefd(np.log(np.array(PHIM)), np.log(np.array(PHID)), ALPHA[0:-1])
  187. #kappa2 = curvatureg(np.log(np.array(PHIM)), np.log(np.array(PHID)))
  188. #kappa = curvature( np.array(PHIM), np.array(PHID))
  189. x = MOD[ np.argmax(kappa) ]
  190. b_pre = np.dot(A, x)
  191. phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
  192. phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
  193. mu1 = ((phid + alpha*phim) / phib)
  194. print ("L-curve selected", alpha, "phid_old", np.sqrt(phid_old/len(b)), "phid", np.sqrt(phid/len(b)), ibreak)
  195. print ("###########################")
  196. if np.sqrt(phid/len(b)) <= 1:
  197. ibreak=0
  198. #plt.figure()
  199. #plt.plot( (np.array(PHIM)), np.log(np.array(PHID)/len(b)), '.-')
  200. #plt.plot( ((np.array(PHIM))[np.argmax(kappa)]) , np.log( (np.array(PHID)/len(b))[np.argmax(kappa)] ), '.', markersize=12)
  201. #plt.axhline()
  202. #plt.plot( np.array(PHIM), np.array(PHID), '.-')
  203. #plt.plot( np.array(PHIM)[np.argmax(kappa)], np.array(PHID)[np.argmax(kappa)], '.', markersize=12)
  204. #plt.savefig('lcurve.pdf')
  205. break
  206. PHIM = np.array(PHIM)
  207. PHID = np.array(PHID)
  208. if (i == MAXITER-1 ):
  209. ibreak = 2
  210. print("Reached max iterations!!", alpha, np.sqrt(phid/len(b)), ibreak)
  211. kappa = curvaturefd(np.log(np.array(PHIM)), np.log(np.array(PHID)), ALPHA[0:-1])
  212. x = MOD[ np.argmax(kappa) ]
  213. b_pre = np.dot(A, x)
  214. phid = np.linalg.norm( np.dot(Wd, (b-b_pre)))**2
  215. phim = np.linalg.norm( np.dot(Phim_base, (x-xr)) )**2
  216. mu1 = ((phid + alpha*phim) / phib)
  217. if lambdastar == "lcurve":
  218. return x, ibreak, np.sqrt(phid/len(b)), PHIM, PHID/len(b), np.argmax(kappa)
  219. else:
  220. return x, ibreak, np.sqrt(phid/len(b))
  221. if __name__ == "__main__":
  222. print("Test")