Surface NMR processing and inversion GUI
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

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., -np.pi, -5, .005] , [1000., np.pi, 5, .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., -np.pi, -5, .005] , [1000., np.pi, 5, .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