Surface NMR processing and inversion GUI
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

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