123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550 |
- import numpy, pylab,array
-
- from rpy2.robjects.packages import importr
-
- import rpy2.robjects as robjects
- import rpy2.robjects.numpy2ri
-
-
- from numpy.fft import fft, fftfreq
-
-
-
- def peakPicker(data, omega, dt):
-
-
-
- window = (2*numpy.pi) / (omega*dt)
-
- data = numpy.array(data)
- peaks = []
- troughs = []
- times = []
- times2 = []
- indices = []
- ws = 0
- we = window
- ii = 0
- for i in range((int)(len(data)/window)):
-
-
-
-
- peaks.append(numpy.max(data[ws:we]))
- times.append( (ws + data[ws:we].argmax()) * dt )
- indices.append( ii + data[ws:we].argmax() )
-
- troughs.append(numpy.min(data[ws:we]))
- times2.append( (ws + (data[ws:we]).argmin()) * dt )
- indices.append( ii + data[ws:we].argmin() )
-
- ws += window
- we += window
- ii += (int)(we-ws)
-
-
-
-
- return (numpy.array(peaks)-numpy.array(troughs))/2., \
- (numpy.array(times)+numpy.array(times2))/2., \
- indices
-
-
-
-
- def regressCurve(peaks,times,sigma2=None ,intercept=True):
-
-
-
- b1 = 0
- b2 = 0
- rT2 = 0.3
- r = robjects.r
-
-
- robjects.globalenv['b1'] = b1
- robjects.globalenv['b2'] = b2
- robjects.globalenv['rT2'] = rT2
-
- value = robjects.FloatVector(peaks)
- times = robjects.FloatVector(numpy.array(times))
-
-
-
-
- if sigma2 != None:
-
-
- my_weights = robjects.FloatVector( sigma2 )
-
-
-
- robjects.globalenv['my_weights'] = my_weights
-
-
-
- if (intercept):
- my_list = robjects.r('list(b1=50, b2=1e2, rT2=0.03)')
- my_lower = robjects.r('list(b1=0, b2=0, rT2=.005)')
- my_upper = robjects.r('list(b1=20000, b2=2000, rT2=.700)')
- else:
- my_list = robjects.r('list(b2=1e2, rT2=0.3)')
- my_lower = robjects.r('list(b2=0, rT2=.005)')
- my_upper = robjects.r('list(b2=2000, rT2=.700)')
-
- my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
-
-
- if (intercept):
-
- fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2)')
-
- else:
- fmla = robjects.Formula('value ~ b2*exp(-times/rT2)')
-
- env = fmla.getenvironment()
- env['value'] = value
- env['times'] = times
-
-
-
-
-
-
-
- Error = False
-
- if (sigma2 != None):
-
-
-
- try:
- fit = robjects.r.tryCatch( robjects.r.nls(fmla, start=my_list, control=my_cont, weights=my_weights, algorithm="port" , \
- lower=my_lower,upper=my_upper))
- except:
- print("regression issue pass")
- Error = True
-
- else:
- try:
- fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port",lower=my_lower,upper=my_upper))
- except:
- print("regression issue pass")
- Error = True
-
- if not Error:
-
- report = r.summary(fit)
- b1 = 0
- b2 = 0
- rT2 = 1
- if (intercept):
- if not Error:
- b1 = r['$'](report,'par')[0]
- b2 = r['$'](report,'par')[1]
- rT2 = r['$'](report,'par')[2]
-
-
-
-
-
- else:
- print("ERROR DETECTED, regressed values set to default")
- b1 = 1e1
- b2 = 1e-2
- rT2 = 1e-2
-
-
-
- return [b1,b2,rT2]
- else:
- if not Error:
- rT2 = r['$'](report,'par')[1]
- b2 = r['$'](report,'par')[0]
- else:
- print("ERROR DETECTED, regressed values set to default")
- return [b2, rT2]
-
- def quadratureDetect(X, Y, tt):
-
- r = robjects.r
-
- robjects.r('''
- Xc <- function(E0, df, tt, phi, T2) {
- E0 * -sin(2*pi*df*tt + phi) * exp(-tt/T2)
- }
-
- Yc <- function(E0, df, tt, phi, T2) {
- E0 * cos(2*pi*df*tt + phi) * exp(-tt/T2)
- }
- ''')
-
-
- Zero = robjects.FloatVector(numpy.zeros(len(X)))
-
-
- E0 = 0.
- df = 0.
- phi = 0.
- T2 = 0.
- robjects.globalenv['E0'] = E0
- robjects.globalenv['df'] = df
- robjects.globalenv['phi'] = phi
- robjects.globalenv['T2'] = T2
- XY = robjects.FloatVector(numpy.concatenate((X,Y)))
-
-
- tt = robjects.FloatVector(numpy.array(tt))
- X = robjects.FloatVector(numpy.array(X))
- Y = robjects.FloatVector(numpy.array(Y))
- Zero = robjects.FloatVector(numpy.array(Zero))
-
-
-
-
- fmla = robjects.Formula('XY ~ c(Xc( E0, df, tt, phi, T2 ), Yc( E0, df, tt, phi, T2 ))')
-
- env = fmla.getenvironment()
- env['Zero'] = Zero
- env['X'] = X
- env['Y'] = Y
- env['XY'] = XY
- env['tt'] = tt
-
-
- start = robjects.r('list(E0=100, df= 0 , phi= 0.00, T2=.100)')
- lower = robjects.r('list(E0=1, df=-13.0, phi= -3.14, T2=.005)')
- upper = robjects.r('list(E0=1000, df= 13.0, phi= 3.14, T2=.800)')
-
- cont = robjects.r('nls.control(maxiter=10000, warnOnly=TRUE, printEval=FALSE)')
-
- fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont, lower=lower, upper=upper, algorithm='port'))
- report = r.summary(fit)
-
-
- E0 = r['$'](report,'par')[0]
- df = r['$'](report,'par')[1]
- phi = r['$'](report,'par')[2]
- T2 = r['$'](report,'par')[3]
-
- return E0,df,phi,T2
-
-
-
-
- def regressSpec(w, wL, X):
-
-
- s = -1j*w
-
-
-
- a = 0
- rT2 = 0.1
- r = robjects.r
-
-
- robjects.globalenv['a'] = a
-
- robjects.globalenv['rT2'] = rT2
- robjects.globalenv['wL'] = wL
- robjects.globalenv['nb'] = 0
-
-
- w = robjects.FloatVector(numpy.array(w))
- XX = robjects.FloatVector(X)
-
-
-
-
-
-
- my_lower = robjects.r('list(a=.001, rT2=.001)')
-
- my_upper = robjects.r('list(a=1.5, rT2=.300)')
-
-
- my_list = robjects.r('list(a=.2, rT2=0.03)')
- my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
-
-
-
-
-
-
- fmla = robjects.Formula('XX ~ a*(.5/rT2) / ((1./rT2)^2 + (w-wL)^2 )')
-
-
- env = fmla.getenvironment()
-
- env['w'] = w
-
-
-
-
- env['XX'] = XX
-
-
- fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port'))
- report = r.summary(fit)
-
-
-
- a = r['$'](report,'par')[0]
- rT2 = r['$'](report,'par')[1]
- nb = r['$'](report,'par')[2]
-
- return a, rT2, nb
-
-
-
- def regressModulus(w, wL, X):
-
-
- s = -1j*w
-
-
-
- a = 0
- rT2 = 0.1
- r = robjects.r
-
-
- robjects.globalenv['a'] = a
- robjects.globalenv['rT2'] = rT2
- robjects.globalenv['wL'] = wL
- robjects.globalenv['nb'] = 0
-
- s = robjects.ComplexVector(numpy.array(s))
- XX = robjects.ComplexVector(X)
- Xr = robjects.FloatVector(numpy.real(X))
- Xi = robjects.FloatVector(numpy.imag(X))
- Xa = robjects.FloatVector(numpy.abs(X))
- Xri = robjects.FloatVector(numpy.concatenate((Xr,Xi)))
-
-
- my_lower = robjects.r('list(a=.001, rT2=.001)')
-
- my_upper = robjects.r('list(a=1.5, rT2=.300)')
-
-
- my_list = robjects.r('list(a=.2, rT2=0.03)')
- my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
-
-
-
-
-
- fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))')
-
- env = fmla.getenvironment()
- env['s'] = s
- env['Xr'] = Xr
- env['Xa'] = Xa
- env['Xi'] = Xi
- env['Xri'] = Xri
- env['XX'] = XX
-
-
- fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port'))
- report = r.summary(fit)
-
-
-
- a = r['$'](report,'par')[0]
- rT2 = r['$'](report,'par')[1]
- nb = r['$'](report,'par')[2]
-
- return a, rT2
-
-
-
- def regressSpecComplex(w, wL, X, known=True, win=None):
-
-
- s = -1j*w
-
-
-
- a = 1
- rT2 = 0.1
- r = robjects.r
- phi2 = 0
- wL2 = wL
-
-
- robjects.globalenv['a'] = a
- robjects.globalenv['rT2'] = rT2
- robjects.globalenv['wL'] = wL
- robjects.globalenv['wL2'] = 0
- robjects.globalenv['nb'] = 0
- robjects.globalenv['phi2'] = phi2
-
- s = robjects.ComplexVector(numpy.array(s))
- XX = robjects.ComplexVector(X)
- Xr = robjects.FloatVector(numpy.real(X))
- Xi = robjects.FloatVector(numpy.imag(X))
- Xa = robjects.FloatVector(numpy.abs(X))
- Xri = robjects.FloatVector(numpy.concatenate((X.real,X.imag)))
-
- robjects.r('''
- source('kernel.r')
- ''')
-
-
-
-
-
-
- if known:
-
- my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14)')
- my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14)')
- my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0)')
- else:
-
- my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14, wL2=wL-5)')
- my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14, wL2=wL+5)')
- my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0, wL2=wL)')
-
- my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
-
-
-
-
-
-
-
-
-
-
- if known:
-
-
- fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL), Kwi(a, phi2, s, rT2, wL) ) ')
- else:
-
-
- fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL2), Kwi(a, phi2, s, rT2, wL2) ) ')
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- env = fmla.getenvironment()
- env['s'] = s
- env['Xr'] = Xr
- env['Xa'] = Xa
- env['Xi'] = Xi
- env['Xri'] = Xri
- env['XX'] = XX
-
-
-
- fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port'))
-
-
-
-
-
-
- report = r.summary(fit)
-
-
-
-
-
-
-
-
-
-
- a = r['$'](report,'par')[0]
- rT2 = r['$'](report,'par')[1]
- nb = r['$'](report,'par')[2]
-
-
-
-
- return a, rT2, nb
-
-
-
-
-
-
- if __name__ == "__main__":
-
- dt = .0001
- T2 = .1
- omega = 2000.*2*numpy.pi
- phi = .0
- T = 8.*T2
-
- t = numpy.arange(0, T, dt)
-
-
-
- data = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) + numpy.random.normal(0,.05,len(t)) \
- + numpy.random.randint(-1,2,len(t))*numpy.random.exponential(.2,len(t))
- cdata = numpy.exp(-t/T2) * numpy.sin(omega * t + phi)
-
-
- sigma2 = numpy.std(data[::-len(data)/4])
-
- print("sigma2", sigma2)
-
- [peaks,times,indices] = peakPicker(data, omega, dt)
-
- [b1,b2,rT2] = regressCurve(peaks,times)
- print("rT2 nonweighted", rT2)
-
- [b1,b2,rT2] = regressCurve(peaks,times,sigma2)
- print("rT2 weighted", rT2)
-
- envelope = numpy.exp(-t/T2)
- renvelope = numpy.exp(-t/rT2)
-
- outf = file('regress.txt','w')
- for i in range(len(times)):
- outf.write(str(times[i]) + " " + str(peaks[i]) + "\n")
- outf.close()
-
- pylab.plot(t,data, 'b')
- pylab.plot(t,cdata, 'g', linewidth=1)
- pylab.plot(t,envelope, color='violet', linewidth=4)
- pylab.plot(t,renvelope, 'r', linewidth=4)
- pylab.plot(times, numpy.array(peaks), 'bo', markersize=8, alpha=.25)
- pylab.legend(['noisy data','clean data','real envelope','regressed env','picks'])
- pylab.savefig("regression.pdf")
-
-
-
- fourier = fft(data)
- pylab.figure()
- freq = fftfreq(len(data), d=dt)
- pylab.plot(freq, (fourier.real))
-
- pylab.show()
-
-
-
-
-
-
-
-
-
-
|