# -*- coding: iso-8859-1 -*-

from numpy import *
from scipy.signal import *
from scipy.optimize import *
from os import path, rename

def result():
  
    measurement = MeasurementResult('Magnetization')
  
    suffix = '' 	# output file name's suffix and...
    counter = 1		# counter for arrayed experiments
        
    # loop over the incoming results: 
    for timesignal in results:
        if not isinstance(timesignal,ADC_Result): 
            continue   
        
        # read experiment parameters:
        pars = timesignal.get_description_dictionary()
        
        # rotate timesignal by current receiver's phase:
        timesignal.phase(pars['rec_phase'])

        # provide timesignal to the display tab:
        data['Current scan'] = timesignal

        # accumulate...
        if not locals().get('accu'):
            accu = Accumulation()

        # skip dummy scans, if any:
        if pars['run'] < 0: continue

        # add up:
        accu += timesignal

        # provide accumulation to the display tab:
        data['Accumulation'] = accu

        # check how many scans are done:
        if accu.n == pars['NS']: # accumulation is complete

            # get number of echoes:
            num_echoes = pars['NECH'] 

            # downsize accu to one point per echo:   
            echodecay = accu + 0  		
            echodecay.x = resize(echodecay.x, int(num_echoes))
            echodecay.y[0] = resize(echodecay.y[0], int(num_echoes))
            echodecay.y[1] = resize(echodecay.y[1], int(num_echoes))

            # specify noise level: 
            if not locals().get('noise'):
                echo = accu.get_accu_by_index(0)
                noise = 0.1*max(abs(echo.y[0]))
                samples = abs(echo.y[0]) > noise

            # set echo times and intensities:
            for i in range(num_echoes):
               # get ith echo:
                echo = accu.get_accu_by_index(i)                          
               # set echo timing:
                echodecay.x[i] = i*2*pars['TAU']
               # set echo value:
                echodecay.y[0][i] = sum(echo.y[0][samples])	# the sum of echo points that exeed noise
                echodecay.y[1][i] = sum(echo.y[1][samples]) 
                #echodecay.y[0][i] = sum(echo.y[0])		# the sum of all echo points
                #echodecay.y[1][i] = sum(echo.y[1])
                #echodecay.y[0][i] = echo.y[0][echo.x.size/2]	# a middle echo point
                #echodecay.y[1][i] = echo.y[1][echo.x.size/2]
                
            # compute a signal's phase:
            phi0 = arctan2(echodecay.y[1][0], echodecay.y[0][0]) * 180 / pi
            if not locals().get('ref'): ref = phi0   
            print 'phi0 = ', phi0

            # rotate signal to maximize Re (optional): 
            #echodecay.phase(-phi0)

            # provide echo decay to the display tab:
            data['Echo Decay'] = echodecay

            # fit a mono-exponential function to the echo decay:
            [amplitude, rate] = fitfunc(echodecay.x, echodecay.y[0])
            print '%s%02g' % ('Amplitude = ', amplitude)
            print '%s%02g' % ('T2 [s] = ', 1./rate)

            # provide the fit to the display tab:
            fit = MeasurementResult('Mono-Exponential Fit')
            for i, key in enumerate(echodecay.x):
                fit[key] = echodecay.y[0][i]
            fit.y = func([amplitude, rate], echodecay.x)
            data[fit.get_title()] = fit
                 
            # check whether it is an arrayed experiment: 
            var_key = pars.get('VAR_PAR')
            if var_key:
                # get variable parameter's value:
                var_value = pars.get(var_key)
                
                # provide data recorded with different var_value's to the display tab:
                data['Accumulation'+"/"+var_key+"=%e"%(var_value)] = accu
                data['Echo Decay'+"/"+var_key+"=%e"%(var_value)] = echodecay
                data[fit.get_title()+"/"+var_key+"=%e"%(var_value)] = fit

                # measure a signal parameter vs. var_value:
                measurement[var_value] = amplitude
                #measurement[var_value] = sum(echodecay.y[0][:])
                #measurement[var_value] = 1./rate

                # provide measurement to the display tab:
                data[measurement.get_title()] = measurement

            # save accu if required:   
            outfile = pars.get('OUTFILE')
            if outfile: 
                datadir = pars.get('DATADIR')
                                
                # write data in Simpson format:
                filename = datadir+outfile+suffix+'.dat'
                if path.exists(filename):
                    rename(filename, datadir+'~'+outfile+suffix+'.dat')
                accu.write_to_simpson(filename)

                # write data in Tecmag format:
                # filename = datadir+outfile+'.tnt'
                # accu.write_to_tecmag(filename, nrecords=20)

               # write parameters in a text file:
                filename = datadir+outfile+suffix+'.par' 
                if path.exists(filename):
                    rename(filename, datadir+'~'+outfile+suffix+'.par')

                fileobject = open(filename, 'w')
                for key in sorted(pars.iterkeys()):
                    if key=='run': continue
                    if key=='rec_phase': continue
                    fileobject.write('%s%s%s'%(key,'=', pars[key]))
                    fileobject.write('\n')
                fileobject.close()

            # reset accumulation:
            del accu

# the fitting procedure:
def fitfunc(xdata, ydata): 

    # initialize variable parameters:
    try: 
        # solve Az = b:
        A = array((ones(xdata.size/2), xdata[0:xdata.size/2]))    
        b = log(abs(ydata[0:xdata.size/2]))     
        z = linalg.lstsq(transpose(A), b)
        amplitude = exp(z[0][0])
        rate = -z[0][1]
    except:
        amplitude = abs(ydata[0])
        rate = 1./(xdata[-1] - xdata[0])
    p0 = [amplitude, rate]

    # run least-squares optimization:
    plsq = leastsq(residuals, p0, args=(xdata, ydata))

    return plsq[0] # best-fit parameters

def residuals(p, xdata, ydata):
    return ydata - func(p, xdata)

# here is the function to fit:
def func(p, xdata):
    return p[0]*exp(-p[1]*xdata)


pass