Commit e2ec5235 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

echo/fit fixes and pylint police

parent c317a0b1

pysen/config.py

+2 −0
Original line number Diff line number Diff line
@@ -234,6 +234,7 @@ def _coverage_check_theta(qmin, lmax, detpos): 
    return theta



def get_theta_pix(x,y,z, sa, ca):

    "get theta angle of a pixel"

    xrot =  x*ca + z*sa

    yrot =  y

    zrot = -x*sa + z*ca
@@ -251,6 +252,7 @@ def _gen_plot_data(out, q, tau, lmax, l1, l2): 


def coverage(lmax, qmin, **kwargs):

    "calculate q-tau coverage"

    #pylint: disable=too-many-locals,too-many-statements

    results = dict()

    results.update(kwargs)

pysen/echo/fit.py

+40 −112
Original line number Diff line number Diff line
@@ -20,8 +20,8 @@ from __future__ import print_function 
import numpy as np



from numpy import sum as npsum

from numpy import (outer, sqrt, pi, exp, log, sinc, cos)

from scipy.optimize import fminbound, fmin

from numpy import (outer, sqrt, pi, exp, log, sinc, cos, sign)

from scipy.optimize import minimize_scalar

from pysen import OMEGA_N

SIGMA2FWHM = 2.0*sqrt(2.0*log(2.0))

#GAUSHAPEFAC = 2.0*np.sqrt(np.log(2.0))
@@ -109,106 +109,45 @@ def flux_weighted_eshape(dj, avlam, delam, flux): 
    omega_avlam = outer(dj, avlam)*OMEGA_N

    return npsum(flux * sinc(omega_delam) * cos(omega_avlam), axis=-1)/npsum(flux)



# MINIMIZERS

# def phase_fit_global2(curr, counts, ecounts, omega, phase_zero, shape_func, **kwargs):

#     """extract symmetry phase (global fit)

#     """

#     niter  = kwargs.get('nsteps', 200 ) # max iterations

#     phtol  = kwargs.get('phtol' , 1e-5) # phase fit tolerance

#     phglb  = kwargs.get('phglb' , 100*phtol)

#     chi_min = np.finfo(float).max

#     ph0_min = None

#     for pha_cur in curr:

#         ph0, chi = phase_fit_local(curr, counts, ecounts, omega,

#                                    pha_cur, shape_func, phtol=phglb)

#         if chi < chi_min:

#             chi_min = chi

#             ph0_min = ph0[0]



#     # refine

#     ph0, chi = phase_fit_local(curr, counts, ecounts, omega,

#                                ph0_min, shape_func,

#                                nsteps=1,

#                                phtol=phtol)

#     return ph0, chi





# def phase_fit_global3(curr, counts, ecounts, omega, shape_func, **kwargs):

#     """extract symmetry phase (global fit)

#     """

#     nsteps = kwargs.get('nsteps', 3)    # how many steps from initial guess

#     phtol  = kwargs.get('phtol' , 1e-5) # phase fit tolerance



#     def func(phase0):

#         """fit function - wrapper around linectract

#         returns chi^2 only """

#         return linear_fit(shape_func((curr-phase0)*omega),

#                           counts, ecounts, chionly=True)



#     pstep = nsteps*(curr[-1]-curr[0])/len(curr)

#     rranges = (slice(curr[0], curr[-1], pstep), )

#     resbrute = brute(func, rranges, args=(), full_output=False, finish=None)



#     print "RES-BRUTE"

#     print resbrute



#     # refine with fmin bound

#     phase_zero = resbrute

#     ph0, chi0, _ierr, _niter  = fminbound(func,  phase_zero-pstep, phase_zero+pstep,

#                                           full_output=True, disp=0, xtol=phtol)



#     return (ph0, phtol), chi0



# def phase_fit_global_basinhopping(curr, counts, ecounts, omega, phase_zero, shape_func, **kwargs):

#     """extract symmetry phase (global fit)

#     """

#     niter  = kwargs.get('nsteps', 200 ) # max iterations

#     phtol  = kwargs.get('phtol' , 1e-5) # phase fit tolerance

#     def func(phase0):

#         """fit function - wrapper around linectract

#         returns chi^2 only """

#         return linear_fit(shape_func((curr-phase0)*omega),

#                           counts, ecounts, chionly=True)

#     min_kw = dict(method="L-BFGS-B", # see [1]

#                   tol=phtol,         # phase tolerance

#                   bounds=((curr[0], curr[-1]),)) # search limits



#     ret = basinhopping(func, x0=phase_zero, minimizer_kwargs=min_kw, niter=niter)

#     return (ret.x[0], phtol), ret.fun





def phase_fit_local(shape_func, dj, counts, ecounts, **kwargs):

    """extract symmetry phase (local fit)

    """

    djlim  = kwargs.get('djlim', (min(dj), max(dj)))

    #nsteps = kwargs.get('nsteps', 3)         # how many steps from initial guess

    tol    = kwargs.get('tolerance' , 1e-10) # phase fit tolerance

    def func(xdj):

        """fit function - wrapper around linectract

        returns chi^2 only """

        return linear_fit(shape_func((dj-xdj)),

                          counts, ecounts, chionly=True)

    dj0, chi0, _, _  = fminbound(func,  djlim[0], djlim[1],

                                 full_output=True, disp=0, xtol=tol)

    return dj0, chi0

    res = minimize_scalar(func,  bounds=(djlim[0], djlim[1]),

                          method='Bounded', options=dict(xatol=tol))

    if not res.success:

        return None, np.inf

    return res.x, func(res.x)



def phase_fit_global(shape_func, dj, counts, ecounts, **kwargs):

    """extract symmetry phase (local fit)

    """

    dj0  = kwargs.get('dj0', 0.0)

    #nsteps = kwargs.get('nsteps', 3)         # how many steps from initial guess

    tol    = kwargs.get('tolerance' , 1e-10) # phase fit tolerance

    def func(xdj):

        """fit function - wrapper around linectract

        returns chi^2 only """

        return linear_fit(shape_func((dj-xdj)),

                          counts, ecounts, chionly=True)

    djmin, chi0, _ , _ , _ = fmin(func, dj0 , full_output=True, disp=0, xtol=tol)

    return djmin[0], chi0







    """extract symmetry phase (global fit)



    """

    djlim    = kwargs.pop('djlim', (min(dj), max(dj)))

    djstep   = kwargs.pop('djstep', abs(djlim[1]-djlim[0]))

    signum   = kwargs.pop('signum', 1) # +1=want positive, -1=want negative, 0=does not matter

    chi0   = np.inf

    dj0    = 0.5*(min(dj)+max(dj))

    nrange = int(np.ceil(0.5*(max(djlim)-min(djlim))/djstep))

    djmin  = min(dj)

    djmax  = max(dj)

    for i in range(-nrange, nrange+1):

        d1 = max(djlim[0]+i*djstep, djmin)

        d2 = min(djlim[1]+i*djstep, djmax)

        xdj, xchi = phase_fit_local(shape_func, dj, counts, ecounts, djlim=(d1,d2), **kwargs)

        if xchi<chi0:

            _, res = linear_fit(shape_func((dj-xdj)), counts, ecounts)

            if sign(signum)!=0 and sign(signum)!=sign(res['slope'][0]):

                continue

            dj0  = xdj

            chi0 = xchi

    return dj0, chi0



def echo_fit(shape_func, dj, counts, ecounts, **kwargs):

    """echo fit
@@ -244,29 +183,18 @@ def echo_fit(shape_func, dj, counts, ecounts, **kwargs): 
def echo_fit_4point(dj, counts, ecounts, **kwargs):

    """IN15 4 point method"""

    #pylint:  disable=unused-argument

    #_ = kwargs.pop('phase_step', 45.0)

    wavelength = kwargs.pop('wavelength', 0)

    n  = counts.shape[0]



    #al = []

    #bl = []

    #for k in range(1,12):

    #    i = counts[k:k+8:2]

    #    ni = len(i)

    #    b  = np.sum(i)/ni

    #    a  = np.sqrt(0.5*np.sum((i-b)**2))

    #    print a, b, k #, len(i), i

    #    al.append(a)

    #    bl.append(b)

    #print "----------------------------"

    #print np.average(al), np.average(bl)





    # FIXME assuming 45-deg step

    i  = counts[n//2-2:n//2+5:2]

    ni = len(i)

    b  = np.sum(i)/ni

    a  = np.sqrt(0.5*np.sum((i-b)**2))

    #print a, b, i

    return dict(amplitude=(a,0.0), average=(b,0.0))

    dj0 = 0

    if wavelength>0:

        p0 = np.arctan2(i[2]-i[0],i[1]-i[3])

        dj0 = p0/(OMEGA_N*wavelength)

    return dict(amplitude=(a,0.0), average=(b,0.0), dj0=(dj0,0))





def echo_curve(shape_func, dj, dj0, amplitude=1.0, average=0.0, npoints=0):

pysen/echo/reduce.py

+140 −138
Original line number Diff line number Diff line 
#!/usr/bin/env python

"reduce echo data"



from functools import partial

import os.path

import time

#from functools import partial



import numpy as np

import h5py



from pysen      import ANGSTROM, NANOSECOND, OMEGA_N

from ..config   import pixel_angle

from ..mathutil import average, histogram1d, normalize_counts

from .fit   import echo_fit, flux_weighted_eshape, echo_curve

from .utils import ( create_mask, calc_sqt, findlevels)

#from pysen      import ANGSTROM, NANOSECOND, OMEGA_N

#from pysen.atarilib import Spectrum, fit_echo_current

#from ..config   import pixel_angle

#from ..mathutil import average, histogram1d, normalize_counts

#from .fit   import echo_fit, flux_weighted_eshape, echo_curve

#from .utils import ( create_mask, calc_sqt, findlevels)





DEFAULT_NRINGS = 5

#DEFAULT_NRINGS = 5





def reduce_data(hdfile, **kwargs):

    "echo plot"

    #

    center_only = kwargs.pop('center_only', False)

    #center_only = kwargs.pop('center_only', False)

    #

    npix    = kwargs.pop('npix', 2)      # pix

    #

    tbin1   = kwargs.pop('tbin1',0)      # TOF bins

    tbin2   = kwargs.pop('tbin2',None)

    #

    min_counts = kwargs.pop('min_counts', 10.0)

    max_chi2   = kwargs.pop('max_chi2'  ,500.0)

    min_amp    = kwargs.pop('min_amp'   ,  0.0)

    #min_counts = kwargs.pop('min_counts', 10.0)

    #max_chi2   = kwargs.pop('max_chi2'  ,500.0)

    #min_amp    = kwargs.pop('min_amp'   ,  0.0)

    # FIXME: HARDCODED

    xpix1   = kwargs.pop('xpix1', 10)

    xpix2   = kwargs.pop('xpix2', 22)

    ypix1   = kwargs.pop('ypix1', 10)

    ypix2   = kwargs.pop('ypix2', 22)

    #xpix1   = kwargs.pop('xpix1', 10)

    #xpix2   = kwargs.pop('xpix2', 22)

    #ypix1   = kwargs.pop('ypix1', 10)

    #ypix2   = kwargs.pop('ypix2', 22)



    base    = os.path.splitext(os.path.basename(hdfile.filename))[0]

    comment = hdfile.attrs['master_comment'].decode("utf-8")

    #comment = hdfile.attrs['master_comment'].decode("utf-8")



    ntaus   = hdfile.attrs['scan_length']

    #ntaus   = hdfile.attrs['scan_length']

    n_idx   = { 'dn': hdfile['/'].attrs['point_to_down'],

                'up': hdfile['/'].attrs['point_to_up'] ,

                'nphases': hdfile['/'].attrs['no_of_phases'] }
@@ -58,15 +61,14 @@ def reduce_data(hdfile, **kwargs): 


        for echo in list(hdfile['/data'].values()):

            print("processing", echo.name)

            pass



            phase   = echo['phase']

            physics = echo['phys']

            params  = echo['params']

            tech    = echo['tech']



            tau0    = physics.attrs['fouriertime']/NANOSECOND

            q0      = physics.attrs['hkl'][0]

            #tau0    = physics.attrs['fouriertime']/NANOSECOND

            #q0      = physics.attrs['hkl'][0]

            lam0    = physics.attrs['lambda']

            phase0  = float(tech.attrs['i5'][0])
@@ -84,120 +86,120 @@ def reduce_data(hdfile, **kwargs): 
        det  = det/pcha

        wlen = np.sum(det, axis=(0,1,2))



        spectrum =  Spectrum(lam=lmax, dlam=dlam, flux=wlen)



        y  = np.sum(det[:,:,:,tbin1:tbin2], axis=(1,2,3))

        ey = np.sqrt(y)

        dn = y[nph:n_idx['up']]

        up = y[n_idx['up']:]

        phase_ef0, (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum,

                                        lam0=lam0, phase0=phase0, phasesens=phasesens)





def process_data(hdfile, dj0=None, qbins=None, tbins=None):

    """

    process echo data

    """

    n_idx = { 'dn': hdfile['/scan'].attrs['point_to_down'],

              'up': hdfile['/scan'].attrs['point_to_up'] }

    preset  = hdfile['/scan'].attrs['preset']



    nt      = hdfile['/detector'].attrs['no_t_channels']

    ny      = hdfile['/detector'].attrs['no_y_channels']

    nx      = hdfile['/detector'].attrs['no_x_channels']



    # create 3-d grid to represent one phase point

    yy, xx, tt  = np.mgrid[0:ny, 0:nx, 0:nt]

    mask = create_mask(xx, yy, tt, tbins, shape='ring', r2=12.0)



    # loop over echos (taus)

    for echo in hdfile['/data'].values():

        phase   = echo['phase']

        physics = echo['phys']

        params  = echo['params']





        # remove old results if any

        if echo.get('results') is not None:

            del echo['results']



        results = echo.create_group('results')

        lam0    = physics.attrs['lambda']

        #tau0    = physics.attrs['fouriertime']



        lmax    = np.array(params['lambdaTable']).flatten()/ANGSTROM

        dlam    = np.ones_like(lmax)*(lmax[1]-lmax[0])



        tau     = np.array(params['fouriertimeTable']).flatten()/NANOSECOND

        tau     = np.tile(tau, (ny, nx)).reshape((ny, nx, nt))



        omega   = float(params.attrs['phaseangle.sensitivities'][0])

        omega   = np.radians(omega)/max(lmax)

        #spectrum =  Spectrum(lam=lmax, dlam=dlam, flux=wlen)



        # calc q

        polar_angle, _ = pixel_angle(xx, yy,

                                     np.radians(float(params.attrs['phi'][0])),

                                     nx=nx, ny=ny)

        q   = 4*np.pi/lmax * np.sin(polar_angle/2.0)

        #y  = np.sum(det[:,:,:,tbin1:tbin2], axis=(1,2,3))

        #ey = np.sqrt(y)

        #dn = y[nph:n_idx['up']]

        #up = y[n_idx['up']:]

        #phase_ef0, (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum,

        #                                             lam0=lam0, phase0=phase0, phasesens=phasesens)



        if qbins is None:

            qbins = findlevels(q.min(), q.max(), DEFAULT_NRINGS+1)





        #det  = np.array([ phase[pp]['detector'].value for pp in phase ])

        #moni = np.array([ phase[pp]['counter' ].value for pp in phase ])

        #pcha = np.array([ phase[pp].attrs['proton_charge'] for pp in phase ])

        #cur  = np.array([ phase[pp].attrs['phase_current'] for pp in phase ])[:, 0]

        det  = phase['detector'][...]

        moni = phase['counter' ][...]

        pcha = phase['proton_charge'][...]

        dj   = np.radians(phase['phase'])/(OMEGA_N*lam0) #*180.0/np.pi)

        ddj = 2*np.abs(dj[1]-dj[0])

        djlim=(dj0 - ddj, dj0 + ddj)



        norm = pcha/float(np.average(pcha))/preset

        moni = np.tile(moni[-1, 2], (nx, ny)).reshape(ny, nx, nt)



        #results.create_dataset('hist', data=create_qhistogram(q, det, mask, qbins))

        results.create_dataset('hist',

                               data=histogram1d( q[mask],

                                                 weights=np.sum(det[:, mask], axis=0),

                                                 bins=np.linspace(qbins[0], qbins[-1],

                                                                  len(qbins)*10) ))



        grp = results.create_group('parameters')

        grp.create_dataset('mask',  data=mask , dtype=int)

        grp.create_dataset('qbins', data=qbins)

        grp.attrs['dj0'] = dj0

        grp.attrs['nrings'] = len(qbins)-1



        # loop over q bins

        for iqsel, (q1, q2) in enumerate(zip(qbins[:-1], qbins[1:])):



            sel_index = np.logical_and(np.logical_and(q1<q, q<q2), mask)



            det_sel  = det[:, sel_index]

            weights  = np.sum(det_sel, axis=0)



            esig = normalize_counts(np.sum(det_sel, axis=-1), norm)

            flux = np.sum(np.sum(moni*sel_index, axis=0), axis=0)



            sfun = partial(flux_weighted_eshape, avlam=lmax-dlam/2.0, delam=dlam, flux=flux)





            res  = echo_fit(sfun, dj, esig[0], esig[1], dj0=dj0, djlim=djlim,

                            fit_global=False)

            efit = echo_curve(sfun, dj, dj0=res['dj0'][0],

                              amplitude=res['amplitude'][0], average=res['average'][0],

                              npoints=101)

            #"write results"

            grp = results.create_group('ring_%03d' % iqsel)

            grp.create_dataset('echo', data=esig)

            grp.create_dataset('flux', data=flux)

            grp.create_dataset('fit',  data=efit)

            grp.attrs['sqt'] = calc_sqt(esig, res['amplitude'], n_idx['up'], n_idx['dn'])

            grp.attrs['tau'] = average(tau[sel_index],  weights)

            grp.attrs['q']   = average(q  [sel_index],  weights)



            for key in res:

                grp.attrs[key] = res[key]
 
#def process_data(hdfile, dj0=None, qbins=None, tbins=None):

#    """

#    process echo data

#    """

#    n_idx = { 'dn': hdfile['/scan'].attrs['point_to_down'],

#              'up': hdfile['/scan'].attrs['point_to_up'] }

#    preset  = hdfile['/scan'].attrs['preset']

#

#    nt      = hdfile['/detector'].attrs['no_t_channels']

#    ny      = hdfile['/detector'].attrs['no_y_channels']

#    nx      = hdfile['/detector'].attrs['no_x_channels']

#

#    # create 3-d grid to represent one phase point

#    yy, xx, tt  = np.mgrid[0:ny, 0:nx, 0:nt]

#    mask = create_mask(xx, yy, tt, tbins, shape='ring', r2=12.0)

#

#    # loop over echos (taus)

#    for echo in hdfile['/data'].values():

#        phase   = echo['phase']

#        physics = echo['phys']

#        params  = echo['params']

#

#

#        # remove old results if any

#        if echo.get('results') is not None:

#            del echo['results']

#

#        results = echo.create_group('results')

#        lam0    = physics.attrs['lambda']

#        #tau0    = physics.attrs['fouriertime']

#

#        lmax    = np.array(params['lambdaTable']).flatten()/ANGSTROM

#        dlam    = np.ones_like(lmax)*(lmax[1]-lmax[0])

#

#        tau     = np.array(params['fouriertimeTable']).flatten()/NANOSECOND

#        tau     = np.tile(tau, (ny, nx)).reshape((ny, nx, nt))

#

#        omega   = float(params.attrs['phaseangle.sensitivities'][0])

#        omega   = np.radians(omega)/max(lmax)

#

#        # calc q

#        polar_angle, _ = pixel_angle(xx, yy,

#                                     np.radians(float(params.attrs['phi'][0])),

#                                     nx=nx, ny=ny)

#        q   = 4*np.pi/lmax * np.sin(polar_angle/2.0)

#

#        if qbins is None:

#            qbins = findlevels(q.min(), q.max(), DEFAULT_NRINGS+1)

#

#

#        #det  = np.array([ phase[pp]['detector'].value for pp in phase ])

#        #moni = np.array([ phase[pp]['counter' ].value for pp in phase ])

#        #pcha = np.array([ phase[pp].attrs['proton_charge'] for pp in phase ])

#        #cur  = np.array([ phase[pp].attrs['phase_current'] for pp in phase ])[:, 0]

#        det  = phase['detector'][...]

#        moni = phase['counter' ][...]

#        pcha = phase['proton_charge'][...]

#        dj   = np.radians(phase['phase'])/(OMEGA_N*lam0) #*180.0/np.pi)

#        ddj = 2*np.abs(dj[1]-dj[0])

#        djlim=(dj0 - ddj, dj0 + ddj)

#

#        norm = pcha/float(np.average(pcha))/preset

#        moni = np.tile(moni[-1, 2], (nx, ny)).reshape(ny, nx, nt)

#

#        #results.create_dataset('hist', data=create_qhistogram(q, det, mask, qbins))

#        results.create_dataset('hist',

#                               data=histogram1d( q[mask],

#                                                 weights=np.sum(det[:, mask], axis=0),

#                                                 bins=np.linspace(qbins[0], qbins[-1],

#                                                                  len(qbins)*10) ))

#

#        grp = results.create_group('parameters')

#        grp.create_dataset('mask',  data=mask , dtype=int)

#        grp.create_dataset('qbins', data=qbins)

#        grp.attrs['dj0'] = dj0

#        grp.attrs['nrings'] = len(qbins)-1

#

#        # loop over q bins

#        for iqsel, (q1, q2) in enumerate(zip(qbins[:-1], qbins[1:])):

#

#            sel_index = np.logical_and(np.logical_and(q1<q, q<q2), mask)

#

#            det_sel  = det[:, sel_index]

#            weights  = np.sum(det_sel, axis=0)

#

#            esig = normalize_counts(np.sum(det_sel, axis=-1), norm)

#            flux = np.sum(np.sum(moni*sel_index, axis=0), axis=0)

#

#            sfun = partial(flux_weighted_eshape, avlam=lmax-dlam/2.0, delam=dlam, flux=flux)

#

#

#            res  = echo_fit(sfun, dj, esig[0], esig[1], dj0=dj0, djlim=djlim,

#                            fit_global=False)

#            efit = echo_curve(sfun, dj, dj0=res['dj0'][0],

#                              amplitude=res['amplitude'][0], average=res['average'][0],

#                              npoints=101)

#            #"write results"

#            grp = results.create_group('ring_%03d' % iqsel)

#            grp.create_dataset('echo', data=esig)

#            grp.create_dataset('flux', data=flux)

#            grp.create_dataset('fit',  data=efit)

#            grp.attrs['sqt'] = calc_sqt(esig, res['amplitude'], n_idx['up'], n_idx['dn'])

#            grp.attrs['tau'] = average(tau[sel_index],  weights)

#            grp.attrs['q']   = average(q  [sel_index],  weights)

#

#            for key in res:

#                grp.attrs[key] = res[key]

pysen/io/writer.py

+2 −0
Original line number Diff line number Diff line
@@ -159,6 +159,7 @@ def generate_phase_point(fd, phase, phase_current, proton_charge, det, mon): 


def generate_phase_scan(fd, scan, detector, physics, tech):

    "generate phase scan"

    #pylint: disable=too-many-locals,too-many-statements

    preset     = scan['preset']

    phase_step = scan['step']

    nphases    = scan['nphases']
@@ -239,6 +240,7 @@ def generate_phase_scan(fd, scan, detector, physics, tech): 


def generate_echo(filename, **kwargs):

    "generate .echo file"

    #pylint: disable=too-many-locals,too-many-statements

    fd = open(filename, 'w+')



    mlog = logging.getLogger()

pysen/revision.py

+2 −2
Original line number Diff line number Diff line
@@ -2,9 +2,9 @@ 
PySEN revision module

"""

import sys

__version__  = "0.6.0"

__version__  = "0.6.1"

__release__  = "dev1"

__date__     = "Feb 6, 2019"

__date__     = "Feb 11, 2019"

#

VERSION  = __version__

RELEASE  = __release__