Surface NMR processing and inversion GUI
您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符

decay.py 27KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720
  1. import numpy, array #,rpy2
  2. from matplotlib import pyplot as plt
  3. import numpy as np
  4. from scipy.optimize import least_squares
  5. from rpy2.robjects.packages import importr
  6. import rpy2.robjects as robjects
  7. import rpy2.robjects.numpy2ri
  8. #import notch
  9. from numpy.fft import fft, fftfreq
  10. # We know/can calculate frequency peak, use this to guess where picks will be.
  11. # maybe have a sliding window that reports peak values.
  12. def peakPicker(data, omega, dt):
  13. # compute window based on omega and dt
  14. # make sure you are not aliased, grab every other peak
  15. window = (2*numpy.pi) / (omega*dt)
  16. data = numpy.array(data)
  17. peaks = []
  18. troughs = []
  19. times = []
  20. times2 = []
  21. indices = []
  22. ws = 0
  23. we = window
  24. ii = 0
  25. for i in range((int)(len(data)/window)):
  26. # initially was just returning this I think avg is better
  27. #times.append( (ws + numpy.abs(data[ws:we]).argmax()) * dt )
  28. peaks.append(numpy.max(data[ws:we]))
  29. times.append( (ws + data[ws:we].argmax()) * dt )
  30. indices.append( ii + data[ws:we].argmax() )
  31. troughs.append(numpy.min(data[ws:we]))
  32. times2.append( (ws + (data[ws:we]).argmin()) * dt )
  33. indices.append( ii + data[ws:we].argmin() )
  34. ws += window
  35. we += window
  36. ii += (int)(we-ws)
  37. #return numpy.array(peaks), numpy.array(times)
  38. # Averaging peaks does a good job of removing bias in noise
  39. return (numpy.array(peaks)-numpy.array(troughs))/2., \
  40. (numpy.array(times)+numpy.array(times2))/2., \
  41. indices
  42. #################################################
  43. # Regress for T2 using rpy2 interface
  44. def regressCurve(peaks,times,sigma2=1,intercept=True):
  45. # TODO, if regression fails, it might be because there is no exponential
  46. # term, maybe do a second regression then on a linear model.
  47. b1 = 0 # Bias
  48. b2 = 0 # Linear
  49. rT2 = 0.3 # T2 regressed
  50. r = robjects.r
  51. # Variable shared between R and Python
  52. robjects.globalenv['b1'] = b1
  53. robjects.globalenv['b2'] = b2
  54. robjects.globalenv['rT2'] = rT2
  55. robjects.globalenv['sigma2'] = sigma2
  56. value = robjects.FloatVector(peaks)
  57. times = robjects.FloatVector(numpy.array(times))
  58. # my_weights = robjects.RVector(value/sigma2)
  59. # robjects.globalenv['my_weigts'] = my_weights
  60. # if sigma2 != 0:
  61. # print "weighting"
  62. # tw = numpy.array(peaks)/sigma2
  63. # my_weights = robjects.RVector( tw/numpy.max(tw) )
  64. # else:
  65. # my_weights = robjects.RVector(numpy.ones(len(peaks)))
  66. # robjects.globalenv['my_weights'] = my_weights
  67. if (intercept):
  68. my_list = robjects.r('list(b1=50, b2=1e2, rT2=0.03)')
  69. my_lower = robjects.r('list(b1=0, b2=0, rT2=.005)')
  70. my_upper = robjects.r('list(b1=20000, b2=2000, rT2=.700)')
  71. else:
  72. my_list = robjects.r('list(b2=1e2, rT2=0.3)')
  73. my_lower = robjects.r('list(b2=0, rT2=.005)')
  74. my_upper = robjects.r('list(b2=2000, rT2=.700)')
  75. my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
  76. if (intercept):
  77. #fmla = robjects.RFormula('value ~ b1 + exp(-times/rT2)')
  78. fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2)')
  79. #fmla = robjects.RFormula('value ~ b1 + b2*times + exp(-times/rT2)')
  80. else:
  81. fmla = robjects.Formula('value ~ b2*exp(-times/rT2)')
  82. env = fmla.getenvironment()
  83. env['value'] = value
  84. env['times'] = times
  85. # ugly, but I get errors with everything else I've tried
  86. my_weights = robjects.r('rep(1,length(value))')
  87. for ii in range(len(my_weights)):
  88. my_weights[ii] *= peaks[ii]/sigma2
  89. Error = False
  90. #fit = robjects.r.nls(fmla,start=my_list,control=my_cont,weights=my_weights)
  91. if (sigma2 != 1):
  92. print("SIGMA 2")
  93. #fit = robjects.r.tryCatch(robjects.r.suppressWarnings(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port", \
  94. # weights=my_weights)), 'silent=TRUE')
  95. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont))#, \
  96. # weights=my_weights))
  97. else:
  98. try:
  99. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port"))#,lower=my_lower,upper=my_upper))
  100. except:
  101. print("regression issue pass")
  102. Error = True
  103. # If failure fall back on zero regression values
  104. if not Error:
  105. #Error = fit[3][0]
  106. report = r.summary(fit)
  107. b1 = 0
  108. b2 = 0
  109. rT2 = 1
  110. if (intercept):
  111. if not Error:
  112. b1 = r['$'](report,'par')[0]
  113. b2 = r['$'](report,'par')[1]
  114. rT2 = r['$'](report,'par')[2]
  115. #print report
  116. #print r['$'](report,'convergence')
  117. #print r['convergence'] #(report,'convergence')
  118. #print r['$'](report,'par')[13]
  119. #print r['$'](report,'par')[14]
  120. else:
  121. print("ERROR DETECTED, regressed values set to default")
  122. b1 = 1e1
  123. b2 = 1e-2
  124. rT2 = 1e-2
  125. #print r['$'](report,'par')[0]
  126. #print r['$'](report,'par')[1]
  127. #print r['$'](report,'par')[2]
  128. return [b1,b2,rT2]
  129. else:
  130. if not Error:
  131. rT2 = r['$'](report,'par')[1]
  132. b2 = r['$'](report,'par')[0]
  133. else:
  134. print("ERROR DETECTED, regressed values set to default")
  135. return [b2, rT2]
  136. #################################################
  137. # Regress for T2 using rpy2 interface
  138. def regressCurve2(peaks,times,sigma2=[None],intercept=True):
  139. if sigma2[0] != None:
  140. my_weights = robjects.FloatVector( sigma2 )
  141. # TODO, if regression fails, it might be because there is no exponential
  142. # term, maybe do a second regression then on a linear model.
  143. b1 = 0 # Bias
  144. b2 = 0 # Linear
  145. bb2 = 0 # Linear
  146. rT2 = 0.3 # T2 regressed
  147. rrT2 = 1.3 # T2 regressed
  148. r = robjects.r
  149. # Variable shared between R and Python
  150. robjects.globalenv['b1'] = b1
  151. robjects.globalenv['b2'] = b2
  152. robjects.globalenv['rT2'] = rT2
  153. robjects.globalenv['bb2'] = b2
  154. robjects.globalenv['rrT2'] = rT2
  155. #robjects.globalenv['sigma2'] = sigma2
  156. value = robjects.FloatVector(peaks)
  157. times = robjects.FloatVector(numpy.array(times))
  158. if (intercept):
  159. my_list = robjects.r('list(b1=.50, b2=1e2, rT2=0.03, bb2=1e1, rrT2=1.3)')
  160. my_lower = robjects.r('list(b1=0, b2=0, rT2=.005, bb2=0, rrT2=.005 )')
  161. my_upper = robjects.r('list(b1=2000, b2=2000, rT2=.700, bb2=2000, rrT2=1.3 )')
  162. else:
  163. my_list = robjects.r('list(b2=.5, rT2=0.3, bb2=.5, rrT2=1.3)')
  164. my_lower = robjects.r('list(b2=0, rT2=.005, bb2=0, rrT2=.005)')
  165. my_upper = robjects.r('list(b2=1, rT2=2.6, bb2=1, rrT2=2.6)')
  166. my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
  167. if (intercept):
  168. #fmla = robjects.RFormula('value ~ b1 + exp(-times/rT2)')
  169. fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2) + bb2*exp(-times/rrT2)')
  170. #fmla = robjects.RFormula('value ~ b1 + b2*times + exp(-times/rT2)')
  171. else:
  172. fmla = robjects.Formula('value ~ b2*exp(-times/rT2) + bb2*exp(-times/rrT2)')
  173. env = fmla.getenvironment()
  174. env['value'] = value
  175. env['times'] = times
  176. # ugly, but I get errors with everything else I've tried
  177. Error = False
  178. #fit = robjects.r.nls(fmla,start=my_list,control=my_cont,weights=my_weights)
  179. if (sigma2[0] != None):
  180. #print("SIGMA 2")
  181. #fit = robjects.r.tryCatch(robjects.r.suppressWarnings(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port", \
  182. # weights=my_weights)), 'silent=TRUE')
  183. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm='port',weights=my_weights,lower=my_lower,upper=my_upper))#, \
  184. # weights=my_weights))
  185. else:
  186. try:
  187. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port"))#,lower=my_lower,upper=my_upper))
  188. except:
  189. print("regression issue pass")
  190. Error = True
  191. # If failure fall back on zero regression values
  192. if not Error:
  193. #Error = fit[3][0]
  194. report = r.summary(fit)
  195. b1 = 0
  196. b2 = 0
  197. rT2 = 1
  198. if (intercept):
  199. if not Error:
  200. b1 = r['$'](report,'par')[0]
  201. b2 = r['$'](report,'par')[1]
  202. rT2 = r['$'](report,'par')[2]
  203. #print report
  204. #print r['$'](report,'convergence')
  205. #print r['convergence'] #(report,'convergence')
  206. #print r['$'](report,'par')[13]
  207. #print r['$'](report,'par')[14]
  208. else:
  209. print("ERROR DETECTED, regressed values set to default")
  210. b1 = 1e1
  211. b2 = 1e-2
  212. rT2 = 1e-2
  213. #print r['$'](report,'par')[0]
  214. #print r['$'](report,'par')[1]
  215. #print r['$'](report,'par')[2]
  216. return [b1,b2,rT2, bb2, rrT2]
  217. else:
  218. if not Error:
  219. rT2 = r['$'](report,'par')[1]
  220. b2 = r['$'](report,'par')[0]
  221. rrT2 = r['$'](report,'par')[3]
  222. bb2 = r['$'](report,'par')[2]
  223. else:
  224. print("ERROR DETECTED, regressed values set to default")
  225. return [b2, rT2, bb2, rrT2]
  226. def fun(x, t, y):
  227. """ Cost function for regression, single exponential, no DC term
  228. x[0] = A0
  229. x[1] = zeta
  230. x[2] = df
  231. x[3] = T2
  232. """
  233. # concatenated real and imaginary parts
  234. pre = np.concatenate((-x[0]*np.sin(2.*np.pi*x[2]*t + x[1])*np.exp(-t/x[3]), \
  235. +x[0]*np.cos(2.*np.pi*x[2]*t + x[1])*np.exp(-t/x[3])))
  236. return y-pre
  237. def fun2(x, t, y):
  238. """ Cost function for regression, single exponential, no DC term
  239. x[0] = A0
  240. x[1] = zeta
  241. x[2] = T2
  242. """
  243. # concatenated real and imaginary parts
  244. pre = np.concatenate((x[0]*np.cos(x[1])*np.exp(-t/x[2]), \
  245. -1.*x[0]*np.sin(x[1])*np.exp(-t/x[2])))
  246. return y-pre
  247. def quadratureDetect2(X, Y, tt, x0="None"):
  248. """ Pure python quadrature detection using Scipy.
  249. X = real part of NMR signal
  250. Y = imaginary component of NMR signal
  251. tt = time
  252. """
  253. print("Pure Python Quad Det")
  254. # df
  255. if x0=="None":
  256. x0 = np.array( [1., 0., 0., .2] )
  257. res_lsq = least_squares(fun, x0, args=(tt, np.concatenate((X, Y))), loss='cauchy', f_scale=1.0,\
  258. bounds=( [1., 0, -13, .005] , [1000., 2*np.pi, 13, .800] ))
  259. else:
  260. res_lsq = least_squares(fun, x0, args=(tt, np.concatenate((X, Y))), loss='cauchy', f_scale=1.0,\
  261. bounds=( [1., 0, -13, .005] , [1000., 2*np.pi, 13, .800] ))
  262. #bounds=( [0., 0, -20, .0] , [1., np.pi, 20, .6] ))
  263. x = res_lsq.x
  264. return res_lsq.success, x[0], x[2], x[1], x[3]
  265. # no df
  266. #x = np.array( [1., 0., 0.2] )
  267. #res_lsq = least_squares(fun2, x, args=(tt, np.concatenate((X, Y))), loss='soft_l1', f_scale=0.1)
  268. #x = res_lsq.x
  269. #return conv, E0,df,phi,T2
  270. #return res_lsq.success, x[0], 0, x[1], x[2]
  271. def quadratureDetect(X, Y, tt, CorrectFreq=False, BiExp=False, CorrectDC=False):
  272. r = robjects.r
  273. if CorrectDC:
  274. robjects.r('''
  275. Xc1 <- function(E01, df, tt, phi, T2_1, DC) {
  276. DC + E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1)
  277. }
  278. Yc1 <- function(E01, df, tt, phi, T2_1, DC) {
  279. DC - E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1)
  280. }
  281. ''')
  282. else:
  283. robjects.r('''
  284. Xc1 <- function(E01, df, tt, phi, T2_1) {
  285. E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1)
  286. }
  287. Yc1 <- function(E01, df, tt, phi, T2_1) {
  288. -E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1)
  289. }
  290. ''')
  291. # bi-exponential
  292. if CorrectDC:
  293. robjects.r('''
  294. Xc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2, DC) {
  295. DC + E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1) +
  296. DC + E02*cos(2*pi*df*tt + phi) * exp(-tt/T2_2)
  297. }
  298. Yc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2, DC) {
  299. DC - E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1) +
  300. DC - E02*sin(2*pi*df*tt + phi) * exp(-tt/T2_2)
  301. }
  302. ''')
  303. else:
  304. robjects.r('''
  305. Xc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2) {
  306. E01*cos(2*pi*df*tt + phi) * exp(-tt/T2_1) +
  307. E02*cos(2*pi*df*tt + phi) * exp(-tt/T2_2)
  308. }
  309. Yc2 <- function(E01, E02, df, tt, phi, T2_1, T2_2) {
  310. -E01*sin(2*pi*df*tt + phi) * exp(-tt/T2_1) +
  311. -E02*sin(2*pi*df*tt + phi) * exp(-tt/T2_2)
  312. }
  313. ''')
  314. # Make 0 vector
  315. Zero = robjects.FloatVector(numpy.zeros(len(X)))
  316. # Fitted Parameters
  317. E01 = 0.
  318. E02 = 0.
  319. df = 0.
  320. phi = 0.
  321. T2_1 = 0.
  322. T2_2 = 0.
  323. DC = 0.
  324. robjects.globalenv['DC'] = DC
  325. robjects.globalenv['E01'] = E01
  326. robjects.globalenv['E02'] = E02
  327. robjects.globalenv['df'] = df
  328. robjects.globalenv['phi'] = phi
  329. robjects.globalenv['T2_1'] = T2_1
  330. robjects.globalenv['T2_2'] = T2_2
  331. XY = robjects.FloatVector(numpy.concatenate((X,Y)))
  332. # Arrays
  333. tt = robjects.FloatVector(numpy.array(tt))
  334. X = robjects.FloatVector(numpy.array(X))
  335. Y = robjects.FloatVector(numpy.array(Y))
  336. Zero = robjects.FloatVector(numpy.array(Zero))
  337. if BiExp:
  338. if CorrectDC:
  339. fmla = robjects.Formula('XY ~ c(Xc2( E01, E02, df, tt, phi, T2_1, T2_2, DC ), Yc2( E01, E02, df, tt, phi, T2_1, T2_2, DC ))')
  340. if CorrectFreq:
  341. start = robjects.r('list(E01=.100, E02=.01, df=0, phi=0. , T2_1=.100, T2_2=.01, DC=0.0)')
  342. lower = robjects.r('list(E01=1e-6, E02=1e-6, df=-50, phi=-3.14 , T2_1=.001, T2_2=.001, DC=0.0)')
  343. upper = robjects.r('list(E01=1.00, E02=1.0, df=50, phi=3.14 , T2_1=.800, T2_2=.8, DC=0.5)')
  344. else:
  345. start = robjects.r('list(E01=.100, E02=.01, phi=0.9 , T2_1=.100, T2_2=.01, DC=0.0)')
  346. lower = robjects.r('list(E01=1e-6, E02=1e-6, phi=-3.14 , T2_1=.001, T2_2=.001, DC=0.0)')
  347. upper = robjects.r('list(E01=1.00, E02=1.0, phi=3.14 , T2_1=.800, T2_2=.8, DC=0.5)')
  348. else:
  349. fmla = robjects.Formula('XY ~ c(Xc2( E01, E02, df, tt, phi, T2_1, T2_2 ), Yc2( E01, E02, df, tt, phi, T2_1, T2_2))')
  350. if CorrectFreq:
  351. start = robjects.r('list(E01=.100, E02=.01, df=0, phi=0. , T2_1=.100, T2_2=.01)')
  352. lower = robjects.r('list(E01=1e-6, E02=1e-6, df=-50, phi=-3.14 , T2_1=.001, T2_2=.001)')
  353. upper = robjects.r('list(E01=1.00, E02=1.0, df=50, phi=3.14 , T2_1=.800, T2_2=.8)')
  354. else:
  355. start = robjects.r('list(E01=.100, E02=.01, phi=0.9 , T2_1=.100, T2_2=.01)')
  356. lower = robjects.r('list(E01=1e-6, E02=1e-6, phi=-3.14 , T2_1=.001, T2_2=.001)')
  357. upper = robjects.r('list(E01=1.00, E02=1.0, phi=3.14 , T2_1=.800, T2_2=.8)')
  358. else:
  359. if CorrectDC:
  360. fmla = robjects.Formula('XY ~ c(Xc1( E01, df, tt, phi, T2_1, DC), Yc1( E01, df, tt, phi, T2_1,DC))')
  361. if CorrectFreq:
  362. start = robjects.r('list(E01=.100, df=0 , phi=0. , T2_1=.100, DC=0.0)')
  363. lower = robjects.r('list(E01=1e-6, df=-50., phi=-3.14, T2_1=.001, DC=0.0)')
  364. upper = robjects.r('list(E01=1.00, df=50. , phi=3.14 , T2_1=.800, DC=0.5)')
  365. else:
  366. start = robjects.r('list(E01=.100, phi= 0. , T2_1=.100, DC=0.0)')
  367. lower = robjects.r('list(E01=1e-6, phi=-3.13, T2_1=.001, DC=0.0)')
  368. upper = robjects.r('list(E01=1.00, phi= 3.13, T2_1=.800, DC=0.5)')
  369. else:
  370. fmla = robjects.Formula('XY ~ c(Xc1( E01, df, tt, phi, T2_1), Yc1( E01, df, tt, phi, T2_1))')
  371. if CorrectFreq:
  372. start = robjects.r('list(E01=.100, df=0 , phi=0. , T2_1=.100)')
  373. lower = robjects.r('list(E01=1e-6, df=-50. , phi=-3.14 , T2_1=.001)')
  374. upper = robjects.r('list(E01=1.00, df=50. , phi=3.14 , T2_1=.800)')
  375. else:
  376. start = robjects.r('list(E01=.100, phi= 0. , T2_1=.100)')
  377. lower = robjects.r('list(E01=1e-6, phi=-3.13, T2_1=.001)')
  378. upper = robjects.r('list(E01=1.00, phi= 3.13, T2_1=.800)')
  379. env = fmla.getenvironment()
  380. env['Zero'] = Zero
  381. env['X'] = X
  382. env['Y'] = Y
  383. env['XY'] = XY
  384. env['tt'] = tt
  385. cont = robjects.r('nls.control(maxiter=10000, warnOnly=TRUE, printEval=FALSE)')
  386. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont, lower=lower, upper=upper, algorithm='port')) #, \
  387. #fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont)) #, \
  388. report = r.summary(fit)
  389. conv = r['$'](fit,'convergence')[0]
  390. #if conv:
  391. # print (report)
  392. # print ("conv", conv)
  393. print ("Conv", r['$'](fit,'convergence')) # T2
  394. print (report)
  395. if BiExp:
  396. if CorrectFreq:
  397. E0 = r['$'](report,'par')[0] # E01
  398. E0 += r['$'](report,'par')[1] # E02
  399. df = r['$'](report,'par')[2] # offset
  400. phi = r['$'](report,'par')[3] # phase
  401. T2 = r['$'](report,'par')[4] # T2
  402. else:
  403. E0 = r['$'](report,'par')[0] # E01
  404. E0 += r['$'](report,'par')[1] # E02
  405. phi = r['$'](report,'par')[2] # phase
  406. T2 = r['$'](report,'par')[3] # T2
  407. else:
  408. if CorrectFreq:
  409. E0 = r['$'](report,'par')[0] # E01
  410. df = r['$'](report,'par')[1] # offset
  411. phi = r['$'](report,'par')[2] # phase
  412. T2 = r['$'](report,'par')[3] # T2
  413. else:
  414. E0 = r['$'](report,'par')[0] # E01
  415. phi = r['$'](report,'par')[1] # phase
  416. T2 = r['$'](report,'par')[2] # T2
  417. #phi = 0.907655876627
  418. #phi = 0
  419. #print ("df", df)# = 0
  420. return conv, E0,df,phi,T2
  421. #################################################
  422. # Regress for T2 using rpy2 interface
  423. def regressSpec(w, wL, X): #,sigma2=1,intercept=True):
  424. # compute s
  425. s = -1j*w
  426. # TODO, if regression fails, it might be because there is no exponential
  427. # term, maybe do a second regression then on a linear model.
  428. a = 0 # Linear
  429. rT2 = 0.1 # T2 regressed
  430. r = robjects.r
  431. # Variable shared between R and Python
  432. robjects.globalenv['a'] = a
  433. robjects.globalenv['rT2'] = rT2
  434. robjects.globalenv['wL'] = wL
  435. robjects.globalenv['nb'] = 0
  436. s = robjects.ComplexVector(numpy.array(s))
  437. XX = robjects.ComplexVector(X)
  438. Xr = robjects.FloatVector(numpy.real(X))
  439. Xi = robjects.FloatVector(numpy.imag(X))
  440. Xa = robjects.FloatVector(numpy.abs(X))
  441. Xri = robjects.FloatVector(numpy.concatenate((Xr,Xi)))
  442. #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
  443. my_lower = robjects.r('list(a=.001, rT2=.001)')
  444. #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
  445. my_upper = robjects.r('list(a=1.5, rT2=.300)')
  446. #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
  447. my_list = robjects.r('list(a=.2, rT2=0.03)')
  448. my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
  449. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  450. ##fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  451. #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
  452. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
  453. fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
  454. env = fmla.getenvironment()
  455. env['s'] = s
  456. env['Xr'] = Xr
  457. env['Xa'] = Xa
  458. env['Xi'] = Xi
  459. env['Xri'] = Xri
  460. env['XX'] = XX
  461. #fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
  462. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  463. report = r.summary(fit)
  464. #print report
  465. #print r.warnings()
  466. a = r['$'](report,'par')[0]
  467. rT2 = r['$'](report,'par')[1]
  468. nb = r['$'](report,'par')[2]
  469. return a, rT2, nb
  470. #################################################
  471. # Regress for T2 using rpy2 interface
  472. def regressSpecComplex(w, wL, X): #,sigma2=1,intercept=True):
  473. # compute s
  474. s = -1j*w
  475. # TODO, if regression fails, it might be because there is no exponential
  476. # term, maybe do a second regression then on a linear model.
  477. a = 1 # Linear
  478. rT2 = 0.1 # T2 regressed
  479. r = robjects.r
  480. phi2 = 0 # phase
  481. wL2 = wL
  482. # Variable shared between R and Python
  483. robjects.globalenv['a'] = a
  484. robjects.globalenv['rT2'] = rT2
  485. robjects.globalenv['wL'] = wL
  486. robjects.globalenv['wL2'] = 0
  487. robjects.globalenv['nb'] = 0
  488. robjects.globalenv['phi2'] = phi2
  489. s = robjects.ComplexVector(numpy.array(s))
  490. XX = robjects.ComplexVector(X)
  491. Xr = robjects.FloatVector(numpy.real(X))
  492. Xi = robjects.FloatVector(numpy.imag(X))
  493. Xa = robjects.FloatVector(numpy.abs(X))
  494. Xri = robjects.FloatVector(numpy.concatenate((X.real,X.imag)))
  495. robjects.r('''
  496. source('kernel.r')
  497. ''')
  498. #Kw = robjects.globalenv['Kwri']
  499. #print (numpy.shape(X))
  500. #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
  501. #my_lower = robjects.r('list(a=.001, rT2=.001)') # Working
  502. my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14, wL2=wL-5)')
  503. #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
  504. my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14, wL2=wL+5)')
  505. #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
  506. my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0, wL2=wL)')
  507. my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
  508. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  509. #fmla = robjects.Formula('Xi ~ Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 ))') # envelope
  510. #fmla = robjects.Formula('Xri ~ c(Re(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )), Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
  511. #fmlar = robjects.Formula('Xr ~ (Kwr(a, phi2, s, rT2, wL)) ') # envelope
  512. #fmlai = robjects.Formula('Xi ~ (Kwi(a, phi2, s, rT2, wL)) ') # envelope
  513. fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL2), Kwi(a, phi2, s, rT2, wL2) ) ') # envelope
  514. #fmla = robjects.Formula('Xri ~ (Kwri(a, phi2, s, rT2, wL)) ') # envelope
  515. #fmla = robjects.Formula('Xa ~ (abs(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
  516. #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
  517. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
  518. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  519. # self.Gw[iw, iT2] = ((np.sin(phi2) * (alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
  520. # (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
  521. # self.Gw[iw, iT2] = ds * self.sc*((np.sin(phi2)*( alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
  522. # (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
  523. # Works Amplitude Only!
  524. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
  525. env = fmla.getenvironment()
  526. env['s'] = s
  527. env['Xr'] = Xr
  528. env['Xa'] = Xa
  529. env['Xi'] = Xi
  530. env['Xri'] = Xri
  531. env['XX'] = XX
  532. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
  533. #fitr = robjects.r.tryCatch(robjects.r.nls(fmlar, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  534. #env = fmlai.getenvironment()
  535. #fiti = robjects.r.tryCatch(robjects.r.nls(fmlai, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  536. #reportr = r.summary(fitr)
  537. #reporti = r.summary(fiti)
  538. report = r.summary(fit)
  539. #print( report )
  540. #exit()
  541. #print( reportr )
  542. #print( reporti )
  543. #exit()
  544. #print r.warnings()
  545. #a = (r['$'](reportr,'par')[0] + r['$'](reporti,'par')[0]) / 2.
  546. #rT2 = (r['$'](reportr,'par')[1] + r['$'](reporti,'par')[1]) / 2.
  547. #nb = (r['$'](reportr,'par')[2] + r['$'](reporti,'par')[2]) / 2.
  548. a = r['$'](report,'par')[0]
  549. rT2 = r['$'](report,'par')[1]
  550. nb = r['$'](report,'par')[2] #phi2
  551. print ("Python wL2", r['$'](report,'par')[3] )
  552. print ("Python zeta", r['$'](report,'par')[2] )
  553. return a, rT2, nb
  554. ###################################################################
  555. ###################################################################
  556. ###################################################################
  557. if __name__ == "__main__":
  558. dt = .0001
  559. T2 = .1
  560. omega = 2000.*2*numpy.pi
  561. phi = .0
  562. T = 8.*T2
  563. t = numpy.arange(0, T, dt)
  564. # Synthetic data, simple single decaying sinusoid
  565. # with a single decay parameter and gaussian noise added
  566. data = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) + numpy.random.normal(0,.05,len(t)) \
  567. + numpy.random.randint(-1,2,len(t))*numpy.random.exponential(.2,len(t))
  568. cdata = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) #+ numpy.random.normal(0,.25,len(t))
  569. #data = numpy.random.normal(0,.25,len(t))
  570. sigma2 = numpy.std(data[::-len(data)/4])
  571. #sigma2 = numpy.var(data[::-len(data)/4])
  572. print("sigma2", sigma2)
  573. [peaks,times,indices] = peakPicker(data, omega, dt)
  574. [b1,b2,rT2] = regressCurve(peaks,times)
  575. print("rT2 nonweighted", rT2)
  576. [b1,b2,rT2] = regressCurve(peaks,times,sigma2)
  577. print("rT2 weighted", rT2)
  578. envelope = numpy.exp(-t/T2)
  579. renvelope = numpy.exp(-t/rT2)
  580. #outf = file('regress.txt','w')
  581. #for i in range(len(times)):
  582. # outf.write(str(times[i]) + " " + str(peaks[i]) + "\n")
  583. #outf.close()
  584. plt.plot(t,data, 'b')
  585. plt.plot(t,cdata, 'g', linewidth=1)
  586. plt.plot(t,envelope, color='violet', linewidth=4)
  587. plt.plot(t,renvelope, 'r', linewidth=4)
  588. plt.plot(times, numpy.array(peaks), 'bo', markersize=8, alpha=.25)
  589. plt.legend(['noisy data','clean data','real envelope','regressed env','picks'])
  590. plt.savefig("regression.pdf")
  591. # FFT check
  592. fourier = fft(data)
  593. plt.figure()
  594. freq = fftfreq(len(data), d=dt)
  595. plt.plot(freq, (fourier.real))
  596. plt.show()
  597. # TODO do a bunch in batch mode to see if T2 estimate is better with or without
  598. # weighting and which model is best.
  599. # TODO try with real data
  600. # TODO test filters (median, FFT, notch)
  601. # It looks like weighting is good for relatively low sigma, but for noisy data
  602. # it hurts us. Check