Surface NMR processing and inversion GUI
選択できるのは25トピックまでです。 トピックは、先頭が英数字で、英数字とダッシュ('-')を使用した35文字以内のものにしてください。

logbarrier.py 11KB

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