Commit 626a3106 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

fixed in phase tab 

- remove calls to obsolete scipy.interpolate.interp[12]d
parent 5a58458c

pysen/echo/__init__.py

+1 −1
Original line number Diff line number Diff line
@@ -4,4 +4,4 @@ echo fit/plot module 


from .fit    import fit_echo_current, Spectrum            # NOQA

from .reduce import get_symmetry_phase, get_phase_indices # NOQA

from .phase_table import make_phase_table                 # NOQA
 
from .phase_table import make_phase_table_classic         # NOQA

pysen/echo/phase_table.py

+97 −84
Original line number Diff line number Diff line 
#!/usr/bin/env python

"PHASE TABLE"



"""

Create and manipulate NSE phase tables

"""

import logging



import numpy as np

import matplotlib.pyplot as plt



from scipy.interpolate import interp2d

from scipy.interpolate import RectBivariateSpline



from pysen import config

from pysen.mathutil import simple_cluster

from ..mathutil import simple_cluster



def brute_force_interpolation(vals, xvals, yvals, x, y):

def interp_bf(vals, xvals, yvals, x, y):

    """

        ORIG CODE FROM nse/src/nsemess.c

    Brute force 2d interpolation. 

    Based on the code from nse/src/nsemess.c

    """

    ix = np.amax(np.asarray(np.where(xvals<=x)))

    iy = np.amax(np.asarray(np.where(yvals<=y)))
@@ -26,98 +25,112 @@ def brute_force_interpolation(vals, xvals, yvals, x, y): 
    res = z1+(z2-z1)*px

    return res



def interp_bivspline(x,y,z):

    """

    wrapper arount RectBivariateSpline

    """

    rinter = RectBivariateSpline(x, y, z.T)

    def func(xnew, ynew):

        return rinter(xnew, ynew).T

    return func



def make_phase_table(phtab, polyfit=0, ncur=10, nphi=11, pos='p2', threshold=0.2):

    "make phase table"

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

    #FIXME: refactor this

    log = logging.getLogger()



    # data must come in same format

    i00_set = simple_cluster(phtab[:,0], threshold=threshold)     # find how many taus we measured

    phtab   = phtab.reshape(-1,len(i00_set), phtab.shape[-1])



    phi_set = np.average(phtab[:,:,1], axis=-1)

def get_phifit(y, cur, acur, phis, deg=2):

    """

    y    - phis where to evalate

    cur  - i00 value where to evaluate

    acur - i00 current polynomial coefficients

    phis - list of phis for acur

    deg  - final polynomial degree 

    """

    p = [ np.polyval(acur[_p], cur) for _p in phis]

    a = np.polyfit(phis, p, deg=deg)

    return np.polyval(a,y)



def check_accuracy(inptab, phasetab, xi00, yphi, func):

    "check accuracy"

    log = logging.getLogger()

    log.info('======> checking ')

    max_abs, max_rel = 0,0

    ni00, kphi, _ = inptab.shape

    for k in range(kphi):

        for i in range(ni00):

            maincur = inptab[i,k,0]

            theta0  = inptab[i,k,1]

            ph_fit  = inptab[i,k,2]

            ph_tab  = func(maincur, theta0)

            ph_int  = interp_bf(phasetab, yphi, xi00, theta0, maincur)

            ph_dif  = ph_int-ph_fit

            ph_rel  = 100*ph_dif/ph_fit

            max_abs = max(max_abs, abs(ph_dif))

            max_rel = max(max_rel, abs(ph_rel))

            log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_tab=%8.5f i5_int=%8.5f diff=(%+.4fA, %+.3f%%)",

                        maincur, theta0, ph_fit, ph_tab, ph_int, ph_dif, ph_rel)

    log.info("max(diff)=(%.4fA,%.3f%%)", max_abs, max_rel)

    return max_abs, max_rel



    i00   = phtab[:,:,0]

    phi   = phtab[:,:,1]

    phase = phtab[:,:,2]

def make_phase_table_classic(inptab, polyfit=0, extent=None, **kwargs):

    """Create phase table 



    if polyfit>0:

        mini00, maxi00 =  0.0, 90.0 # FIXME: hard coded

        minphi, maxphi = -2.0, config.phi_limits(pos)[1]

    else:

    inptab   - input table with three columns i00 [Amp], phi [deg], phase current [Amp]

    polyfit     -

    ncur, nphi  -

    pos         -

    threshold   -

    """

    threshold  = kwargs.pop('threshold',  0.2)

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

    nphiangles = kwargs.pop('nphiangles', 11)

    # i00_set - 'set' values of i00

    # phi_set - 'set' values of phi (scattering angle)

    i00_set = simple_cluster(inptab[:,0], threshold=threshold)     # find how many taus we measured

    inptab   = inptab.reshape(-1,len(i00_set), inptab.shape[-1])

    phi_set = np.average(inptab[:,:,1], axis=-1)

    i00   = inptab[:,:,0]

    phi   = inptab[:,:,1]

    phase = inptab[:,:,2]



    #get the extent of the interpolation/extrapolation table

    if extent is None:

        mini00, maxi00 = np.amin(i00), np.amax(i00)

        minphi, maxphi = np.amin(phi), np.amax(phi)

    else:

        mini00, maxi00, minphi, maxphi = extent



    xi00  = np.linspace(mini00, maxi00, ncur)

    yphi  = np.linspace(minphi, maxphi, nphi)



    plt.subplot(1,2,1)



    # plot measured data

    for i, c in enumerate(i00_set):

        plt.plot(phi[:,i],phase[:,i],'s', label=r'i$_{00}$=%.1f' % c)

    xi00  = np.linspace(mini00, maxi00, ncurrents)

    yphi  = np.linspace(minphi, maxphi, nphiangles)



    if polyfit>0: # fit/extratpolate

    res = { 'pha': phase, 'phi': phi , 'i00': i00_set }

    if polyfit<=0:  # interpolate using measured data

        f = interp_bivspline( i00_set, phi_set, phase)

        x, y  = np.meshgrid(i00_set, phi_set)

        z = f(i00_set, phi_set)

    else:           # fit/extratpolate

        x, y = np.meshgrid(xi00, yphi)

        z    = np.empty_like(x)



        # fit phase = f(i00) for a given phi angle

        # this is fixed to 2nd degree polynomial

        afit_i00 = {}

        for j, p in enumerate(phi_set):

            afit_i00[p]    = np.polyfit(i00[j,:],phase[j,:], deg=2)



        # fit/extrapolate phase to

        for i, c in enumerate(xi00):

            # values of symmetry phase for given current

            pha_fit = [ np.polyval(afit_i00[_phi], c) for _phi in phi_set]

            # fit phase = f(phi_angle) for a given current

            a = np.polyfit(phi_set, pha_fit, deg=polyfit)

            # fill in data

            z[:,i] = np.polyval(a,yphi)

            plt.plot(yphi, np.polyval(a,yphi), '--', lw=1, label=r'i$_{00}$=%.1f' % c)

        f = interp2d( x, y, z, kind='linear', fill_value=None, bounds_error=False)

    else:

        # interpolate using measured data

        f = interp2d(i00, phi, phase, kind='linear', fill_value=None, bounds_error=False)

        x, y  = np.meshgrid(i00_set, phi_set)

        z = f(i00_set, phi_set)

            z[:,i]  = get_phifit(yphi, c, afit_i00, phi_set, deg=polyfit)



    phase_table =  f(xi00,yphi)

        res['exp'] = {}

        for i,c in enumerate(i00_set):

            v  = get_phifit(yphi, c, afit_i00, phi_set, deg=polyfit)

            res['exp'][c] = (yphi, v)

        f = interp_bivspline( xi00, yphi, z)



    res['cont'] = (x,y,z)

    phase_table =  f(xi00,yphi)



    max_rel = 0

    max_abs = 0

    res['accuracy'] = 0, 0

    if polyfit>0:

        log.info('======> checking ')

        ni00, kphi, _ = phtab.shape

        for k in range(kphi):

            for i in range(ni00):

                maincur = i00[i,k]

                theta0  = phi[i,k]

                ph_fit  = phase[i,k]

                ph_tab  = f(maincur, theta0)

                ph_int  = brute_force_interpolation(phase_table, yphi, xi00, theta0, maincur)

                ph_dif  = ph_int-ph_fit

                ph_rel  = 100*ph_dif/ph_fit

                max_abs = max(max_abs, abs(ph_dif))

                max_rel = max(max_rel, abs(ph_rel))

                log.info("i00=%9.5f phi=%9.5f, i5_fit=%8.5f i5_tab=%8.5f i5_int=%8.5f diff=(%+.4fA, %+.3f%%)",

                        maincur, theta0, ph_fit, ph_tab, ph_int, ph_dif, ph_rel)

        log.info("max(diff)=(%.4fA,%.3f%%)", max_abs, max_rel)



    plt.xlabel(r'$\phi$ [deg]')

    plt.ylabel(r'i$_5$ [A]')

    plt.legend(loc='best')

    plt.grid()



    plt.subplot(1,2,2)

    plt.contourf(y,x,z,21)

    plt.colorbar()

    plt.xlabel(r'$\phi$ [deg]')

    plt.ylabel(r'i$_{00}$ [A]')

    #

    plt.suptitle('Phase Table Plot')

    return xi00, yphi, phase_table, (max_abs, max_rel)
 
        res['accuracy'] = check_accuracy(inptab, phase_table, xi00, yphi, f)

    res['xi00'] = xi00

    res['yphi'] = yphi

    return phase_table, res

#EOF

pysen/mathutil.py

+3 −1
Original line number Diff line number Diff line
@@ -83,7 +83,9 @@ def linear_fit(t, v): 
    return a, b, chi2, r, maxd



def simple_cluster(data, threshold=0.1):

    "simple cluster"

    """find 'clusters' of data based on a threshold

    returns an array of averages of the clusters

    """

    data = sorted(data)

    v0   = data[0]

    indices = [0,]

pysen/physics/water_model.py

+37 −31
Original line number Diff line number Diff line
@@ -7,15 +7,13 @@ H2O/D2O Model Parameters 
"""



from scipy.special import spherical_jn

from   scipy         import interpolate

from numpy import asarray, arange, ones_like, newaxis, sqrt, exp, log

from numpy import asarray, arange, ones_like, newaxis, sqrt, exp, log, interp

from numpy import sum as npsum





# Model parameters, see Table

ns   = 1000.0       # unit of time is ps in this model

kB   = 0.08617      # meV/K

hbar = 0.6582  # meV/ps

hbar = 0.6582119569 # meV*ps



# Atomic mass

M_H = 0.1044 # meV * ps^2 / A^2
@@ -325,38 +323,46 @@ tau1_0_tx = 0.0485 # ps 
EA_tx     = 1850      # cal/mol

RA        = 1.9872036 # cal/mol/K



# TABLE I in [3]

# TABLE I in [3], sorted in increasing temp

Teixeira_Data = asarray([

    #T[K]	     t0[ps]	L[A]

    273.15+20.0, 1.25,	1.29,

    273.15+12.0, 1.66,	1.25,

    273.15 +5.0, 2.33,	1.32,

    273.15 -5.0, 4.66,	1.54,

    273.15-10.0, 6.47, 	1.65,

    273.15-12.0, 7.63, 	1.70,

    273.15-15.0, 8.90, 	1.73,

    273.15-20.0, 22.7, 	2.39,

    273.15-17.0, 10.8, 	1.80,

    273.15-20.0, 22.7, 	2.39]).reshape((-1,3))

    273.15-15.0, 8.90, 	1.73,

    273.15-12.0, 7.63, 	1.70,

    273.15-10.0, 6.47, 	1.65,

    273.15 -5.0, 4.66,	1.54,

    273.15 +5.0, 2.33,	1.32,

    273.15+12.0, 1.66,	1.25,

    273.15+20.0, 1.25,	1.29,

    ]).reshape((-1,3))



t0_tx = interpolate.interp1d(Teixeira_Data[:,0], Teixeira_Data[:,1], fill_value='extrapolate')

L_tx  = interpolate.interp1d(Teixeira_Data[:,0], Teixeira_Data[:,2], fill_value='extrapolate')

def tau1_Teixeira(T):

    "tau1 [ps]"

    return tau1_0_tx * exp(EA_tx/(RA*T)) # Eq.(6) in [3]



def W_Teixeira(Q):

def U_Teixeira(Q):

    "Debye-Waller component, Eq.(1) in [3]"

    return exp(-(Q*u_Debye)**2/3.0)



def R_Teixeira(t, Q, T):

    "Rotational diffusion component, Eq.(2) in [3]"

    jn   = spherical_jn(2,Q*R_HO)

    tau1 = tau1_0_tx * exp(EA_tx/(RA*T)) # Eq.(6) in [3]

    jn   = spherical_jn([0,1,2],Q*R_HO)

    tau1 = tau1_Teixeira(T)

    return jn[0]**2 + 3*jn[1]**2*exp(-t/(3*tau1)) + 5*jn[2]**2*exp(-t/tau1)



def Gamma_Teixeira(Q, T):

    """Translational diffusion width Eq.(3) in [3]

    in units of 1/ps

    """

    t0 = interp(T, Teixeira_Data[:,0], Teixeira_Data[:,1])

    L  = interp(T, Teixeira_Data[:,0], Teixeira_Data[:,2])

    DQ2 = L**2/(6*t0)*Q**2

    return  DQ2/(1+DQ2*t0)



def T_Teixeira(t, Q, T):

    "Translational diffusion component Eq.(3) in [3]"

    t0 = t0_tx(T)

    D  = L_tx(T)**2/(6*t0)

    DQ2 = D*Q**2

    return exp(-t*DQ2/(1+DQ2*t0))

    return exp(-Gamma_Teixeira(Q,T)*t)





#
@@ -382,11 +388,11 @@ def em_model(t, Q, T, components='all'): 


def tx_model(t, Q, T, components='all'):

    "Teixeira Model"

    w = W_Teixeira(Q)*ones_like(t)

    r = R_Teixeira(t, Q, T)

    t = T_Teixeira(t, Q, T)

    res = { 'all' : w*r*t,

            'w'   : w,

    u = U_Teixeira(Q)*ones_like(t)  # Debye-Waller

    r = R_Teixeira(t, Q, T)         # rotational

    t = T_Teixeira(t, Q, T)         # translational

    res = { 'all' : u*r*t,

            'u'   : u,

            'r'   : r,

            't'   : t}

    return res.get(components)

pysen/plot/miscplotlib.py

+74 −37
Original line number Diff line number Diff line
@@ -5,69 +5,106 @@ import sys 
import time

import logging



from io import StringIO



import numpy as np

import matplotlib.pyplot as plt

import h5py



from ..config import phi_limits

from ..inout import convert_files

from ..echo  import get_symmetry_phase, make_phase_table

from ..echo  import get_symmetry_phase, make_phase_table_classic



def plot_phase_table(res):

    """ plot phase table"""

    plt.subplot(1,2,1)

    #

    phi     = res['phi']

    pha     = res['pha']

    x, y, z = res['cont']



    for i, c in enumerate(res['i00']):

        # show where we measured data

        if res.get('exp') is None:

            plt.plot(phi[:,i],pha[:,i],'s--', lw=1, ms=4, label=rf'i$_{{00}}$={c:.1f}')

        else:

            yphi, vals = res['exp'][c]

            p = plt.plot(phi[:,i],pha[:,i],'s', lw=1, ms=4)

            plt.plot(yphi, vals, '--', lw=1, color=p[0].get_color(), label=rf'i$_{{00}}$={c:.1f}')



    plt.xlabel(r'$\phi$ [deg]')

    plt.ylabel(r'i$_5$ [A]')

    plt.legend(loc='best')

    plt.grid()



    plt.subplot(1,2,2)

    plt.contourf(y,x,z,21)

    plt.colorbar()

    plt.xlabel(r'$\phi$ [deg]')

    plt.ylabel(r'i$_{00}$ [A]')

    #

    plt.suptitle('Phase Table Plot')



def print_phase_table(phase, res, pos='p2', filename=None):

    "print phasetable"

    with StringIO() as strbuf:

        strbuf.write( "input phasetable\n")

        strbuf.write(f"phasetable pos{pos[-1]} rebuild\n")

        strbuf.write(f"c time: {time.asctime()}\n")

        accuracy = res['accuracy']

        strbuf.write(f"c max interpolation error [{accuracy[0]:.3g}A, {accuracy[1]:.3f}%%]\n")

        strbuf.write( "c command used:\n")

        strbuf.write(f"c {' '.join(sys.argv)}\n")

        strbuf.write("*\t angles ")

        i00 = res['xi00']

        phi = res['yphi']

        for p in phi:

            strbuf.write(f"{p:7.4f} ")

        strbuf.write("\n")

        strbuf.write("c amperes(main) --- phase current(i5)\n")

        for i,c in enumerate(i00):

            strbuf.write(f"*\t{c:7.3f} ")

            for j,_ in enumerate(phi):

                strbuf.write(f"{phase[j,i]:7.4f} ")

            strbuf.write("\n")

        strbuf.write("eof\n")

        if filename is None:

            print(strbuf.getvalue())

        else:

            with open(filename, 'w+', encoding='ascii') as fdout:

                fdout.write(strbuf.getvalue())



def action_phase_table(filenames, **kwargs):

    "action for phase table plots"

    phtab = None

    log = logging.getLogger()

    outdir    = kwargs.pop('outdir')

    savefig   = kwargs.pop('savefig')

    savefile  = kwargs.pop('savefile')

    polyfit   = kwargs.pop('polyfit')

    pos       = kwargs.pop('pos')

    threshold = kwargs.pop('threshold')

    #

    log   = logging.getLogger()

    inptab = None

    for filename in convert_files(*filenames, converter='from_taco', ignore_errors=True, outdir=outdir):

        try:

            basename, _ = os.path.splitext(os.path.basename(filename))

            log.info('======> input data %s', basename)

            with h5py.File(filename, 'r') as hdf5file:

                res = get_symmetry_phase(hdf5file, **kwargs)

                if phtab is None:

                    phtab = res

                if inptab is None:

                    inptab = res

                else:

                    phtab = np.vstack((phtab,res))

                    inptab = np.vstack((inptab,res))

        except (OSError, RuntimeError) as exc:

            log.error("%s",exc)

            return

    # FIXME hardcoded constants

    extent = (0.0, 90.0, -2.0, phi_limits(pos)[1]) if polyfit>0 else None



    plt.figure(figsize=(12,5))

    i00, phi, phase, max_diff = make_phase_table(phtab, polyfit, pos=pos, threshold=threshold)



    phase, res = make_phase_table_classic(inptab, polyfit, extent=extent, threshold=threshold)

    print_phase_table(phase, res, pos=pos, filename=kwargs.get('savefile'))

    plot_phase_table(res)

    if savefig:

        plt.savefig(savefig)



    def print_phase_table(stream=None):

        if stream:

            orig_stdout = sys.stdout

            sys.stdout  = stream

        print("input phasetable")

        print("phasetable pos%s rebuild" % pos[-1])

        print("c time: ", time.asctime())

        print("c max interpolation error [%.3gA, %.3f%%]" % (max_diff[0],max_diff[1]))

        print("c command used:")

        print("c %s" % " ".join(sys.argv))

        print("*\t angles ", end='')

        for p in phi:

            print("%7.4f " % p, end='')

        print("\nc amperes(main) --- phase current(i5)")

        for i,c in enumerate(i00):

            print("*\t%7.3f " % c,end='')

            for j,_ in enumerate(phi):

                print("%7.4f " % phase[j,i], end='')

            print()

        print("eof")

        if stream:

            sys.stdout = orig_stdout



    if savefile:

        with open(savefile, 'w+', encoding='ascii') as fd:

            print_phase_table(fd)

    else:

        print_phase_table()

#EOF