Surface NMR processing and inversion GUI
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

logbarrier-nnls.py 13KB

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