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

logbarrier.py 12KB

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