python cleanup

macro

! ==== set paths

! datapath ./examples

! savepath ./out

set r.min_counts       100.00    ! min counts per pixel, default is   1.00

set r.min_sigma        0.03      ! min sigma per pixel (relative) default is   0.01

set r.ave_min          0.00      ! min average value default is   0.00

set r.ave_maxdev       0.25      ! default is   1.50

set r.fqt_min          0.10      ! default is   0.01

set r.fqt_maxsigma     1.00      ! default is 100.00

set r.phase_tol        1.0E-10   ! default is   1.0E-9

set r.phase_maxdev     asin(1.0) ! symmetry phase search limits, default is 1.27rad

!

set r.max_chisq        300.00    ! default is 1000.0

set r.postcoll_qcat    0.00      ! default is 0.00

set r.max_field_var    2.50E-5   ! default is 1e-5

! ==== set binning

bins pix nbins 8

bins tof custom 2 13 26 39

! ==== set histogramming

set A  1e-10

set ns 1e-9

histo tau  nbins 20 min 0.002*ns max 200.0*ns log

!histo q    nbins 20 min 0.000/A  max 0.5/A

histo q custom 0.0484/A 0.0594/A 0.0674/A 0.0746/A 0.0818/A 0.0894/A 0.0982/A 0.1101/A 0.1282/A 0.1552/A 0.1962/A 0.3381/A  ! uniform binning

c ==== read data

read s5848.echo s5849.echo s5850.echo s5851.echo s5852.echo s5853.echo as resolution

read s5854.echo s5855.echo s5857.echo s5858.echo s5859.echo s5860.echo s5861.echo as sample

c === process

match all

fit   res

fit   sam flag offset

collect

!collect

!report start

!report dir details

!report sqtplot

!report end nopdf

c === plot

!plot sqt

all: build

build:

	@echo '# -- AUTOMATICALLY CREATED DO NOT EDIT'        >  $(VERSION_FILE)

	@echo '"-- AUTOMATICALLY CREATED DO NOT EDIT"'        >  $(VERSION_FILE)

	@echo '#pylint: disable=invalid-name'                 >> $(VERSION_FILE)

	@echo 'version = "$(VERSION_MAJOR).$(VERSION_MINOR)"' >> $(VERSION_FILE)

	@echo 'release = "$(VERSION_RELEASE)"'                >> $(VERSION_FILE)

	@echo "PYTHON version: $(PYTHON)"

#

"""

python module to aid plotting drspine results

"""

__all__     = ['plot_sqt', 'plot_sqt_bin', 'plot_chi2',

               'plot_q', 'plot_qvec', 'plot_qmap', 'plot_elephant',

               'plot_fits', 'plot_maps', 'plot_fit_pixel', 'plot_fit_map',

Python interface to DrSPINE

"""

import sys, os

import ast

from collections import OrderedDict

import numpy as np

import scipy.stats as stats

import matplotlib.pyplot as plt

from collections import OrderedDict

from matplotlib import cm

from scipy.stats import stats

from scipy.optimize import curve_fit #, OptimizeWarning

#NS = 1e-9

def average(arr, axis=None):

    "calculate the average alog an axis"

    marr = np.ma.masked_array(arr,np.isnan(arr))

    mavg = np.mean(marr,axis=axis)

    return mavg.filled(np.nan)

    return value

def _remove_empty_lines(lines, strip=True):

    "remove empty lines"

    for line in lines:

        stripped_line = line.strip()

        if stripped_line:

            else:

                yield line

def load_txt(filename, delimiter=',', cleanup=[]):

    """ """

def load_txt(filename, delimiter=',', cleanup=None):

    """wrapper around np.loadtxt """

    with open(filename) as fd:

        data = fd.readlines()

    # remove lines matching cleanup

    try:

        for remove in cleanup:

            data = [ line for line in data if not remove in line]

    except TypeError:

        pass

    return np.loadtxt(data, delimiter=delimiter)

    #s = np.loadtxt(filename, delimiter=',')

def fit_sqt(q, data, model='diffusion', **kwargs):

    """

    returns a tuple:

        parameter, parameter_error, fit_function

    """

    A       = kwargs.pop('A',  1.0)

    G       = kwargs.pop('G',  1.0)

    t0      = kwargs.pop('t0', 1.0)

    #pylint: disable=too-many-locals, invalid-name

    A       = kwargs.pop('A',  1.0) #pylint: disable=unused-variable

    G       = kwargs.pop('G',  1.0) #pylint: disable=unused-variable

    t0      = kwargs.pop('t0', 1.0) #pylint: disable=unused-variable

    beta    = kwargs.pop('beta', 2.0/3.0)

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

    p0      = kwargs.pop('initvals', None)

        'kww'            : (2, lambda t, t0, beta:       np.exp(-(t/t0)**beta) ), # b-fitting

        'kww-fix'        : (3, lambda t, t0:             np.exp(-(t/t0)**beta) ), # no b-fitting

        'zilman-granek'  : (1, lambda t, G:              np.exp(-(G*t)**beta)  ), # no b-fitting

        'diffusion'      : (1, lambda t, G:              np.exp(-G*q*q*t)      ), # simple diffusion

        'diffusion-norm' : (1, lambda t, G, A:         A*np.exp(-G*q*q*t)      ), # simple diffusion * norm

        'diffusion+const': (2, lambda t, G, A:     (1-A)*np.exp(-G*q*q*t)+A    ), # simple diffusion + const

        'diffusion'      : (1, lambda t, G:              np.exp(-G*q*q*t)      ), # diffusion

        'diffusion-norm' : (1, lambda t, G, A:         A*np.exp(-G*q*q*t)      ), # diffusion/norm

        'diffusion+const': (2, lambda t, G, A:     (1-A)*np.exp(-G*q*q*t)+A    ), # diffusion+const

        #

        'exp'            : (2, lambda t, t0, A:        A*np.exp(-(t/t0))    ), #

        'exp+const'      : (2, lambda t, t0, A:    (1-A)*np.exp(-(t/t0))+A  ), #

    }

    npar, func = function_library.get(model.lower(), (0, None))

    if func is None: raise RuntimeError("model '%s' is invalid" % model)

    if func is None:

        raise RuntimeError(f"model {model} is invalid")

    default_bounds = {

        'diffusion-norm': ((0.0,0.0,),(np.inf, 2.0)),

    tau  = data[:,0]

    sqt  = data[:,1]

    serr = data[:,2]

    popt, pcov = curve_fit(func, tau[tau<max_tau], sqt[tau<max_tau], sigma=serr[tau<max_tau], absolute_sigma=True,

                           p0=p0, bounds=bounds, maxfev=10000, ftol=1e-5)

    #pylint: disable=unbalanced-tuple-unpacking

    popt, pcov = curve_fit(func, tau[tau<max_tau], sqt[tau<max_tau], sigma=serr[tau<max_tau],

                absolute_sigma=True, p0=p0, bounds=bounds, maxfev=10000, ftol=1e-5)

    par, epar  = popt[0], np.sqrt(pcov[0,0])

    fit_func = lambda x : func(x, *popt)

    chisq = np.sum(((fit_func(tau[tau<max_tau]) - sqt[tau<max_tau])/serr[tau<max_tau])**2)/len(sqt[tau<max_tau] - npar)

    sel = tau<max_tau

    chisq = np.sum(((fit_func(tau[sel])-sqt[sel])/serr[sel])**2)/len(sqt[sel] - npar)

    return (par,epar, fit_func, chisq, popt)

                # -------- data section ------------

                dataheader=False

                for line in lines:

                    if line == '#eod' or line == '#nxt':

                    if line in ('#eod','#nxt'):

                        result.append((meta, np.asarray(data)))

                        data = []

                        meta = OrderedDict()

    filename - datreat type file

    """

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

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

    model   = kwargs.pop('fit_model', '')   # fit model (see fit_sqt)

    summary = kwargs.pop('show_summary', False)

    fitres = []

    if model.startswith('diffusion'):

        print("# q     dq      D          dD         chi2  extra_pars")

        print("# q     dq      D          dD         chi2\textra_pars")

    elif model.startswith('zilman-granek'):

        print("# q     dq      Gamma      dGamma     chi2  extra_pars")

        print("# q     dq      Gamma      dGamma     chi2\textra_pars")

    elif model:

        print("# q     dq      tau0       dtau0      chi2  extra_pars")

        print("# q     dq      tau0       dtau0      chi2\textra_pars")

    for meta, data in result:

        if not len(data): continue

        if not len(data):

            continue

        q  = meta.get('q')

        dq = meta.get('q_var')

        if not q or q=='NaN': continue

        if not dq: continue

        if not q or q=='NaN' or not dq:

            continue

        q  = float(q)  #*ANGSTROEM

        dq = float(dq) #*ANGSTROEM

        label = r"Q=(%.3f $\pm$ %.3f) $\AA^{-1}$" % (q, dq)

        sdat = _cleanup(data)

        if sdat is None: continue

        if sdat is None:

            continue

        #

        tau  = sdat[:,0]

        sqt  = sdat[:,1]

        dsqt = abs(sdat[:,2])

        p = plt.errorbar(tau, sqt, xerr=dtau, yerr=dsqt, fmt=fmt, label=label)

        color=p[0].get_color()

        if (model):

        if model:

            try:

                par0, epar0, fit_func, chi2, popt = fit_sqt(q, sdat, model=model, **mpars)

                ftau = np.logspace(np.log10(min(tau))-0.1, np.log10(max(tau))+0.1)

                fsqt = fit_func(ftau)

                plt.plot(ftau, fsqt, '--', color=color) # fit

                fitres.append((q,dq,par0,epar0))

                print("%.4f, %.4f, %.3e, %.3e, %.3g,\t" % (q, dq, par0, epar0, chi2), end="")

                print(f"{q:.4f}, {dq:.4f}, {par0:.3e}, {epar0:.3e}, chi={chi2:.3g},", end="\t")

                if len(popt)>1:

                    print(", ".join(["%.3g" % _x for _x in popt[1:]]), end=",")

                    print("".join([f"{_:.3g}, " for _ in popt[1:]]), end=" ")

                print("")

            except RuntimeError as _exc:

                print("fitting: %s" % _exc)

                print(f"fitting error: {_exc}")

        else:

            plt.plot(tau, sqt, ':', color=color)

        plt.xlabel(r'$Q \; [\AA^{-1}]$')

        #plt.legend(loc='best')

        plt.grid(True, which='both')

        plt.title("Model: %s" % model)

        plt.title(f"Model: {model}")

# #####################################################################

        plt.savefig(filename+".pdf")

def _plot_chi2_histo(bins, chi2, color='r', fmt='.', label=None):

    db   = (bins[1:]-bins[:-1])     # bin width

    db   =  bins[1:]-bins[:-1]      # bin width

    b    = (bins[1:]+bins[:-1])/2.0 # bin center

    N=float(sum(chi2))*db

    plt.step(    b,  chi2/N, where='mid', color=color, ls='-')

    plt.errorbar(b,  chi2/N, yerr=np.sqrt(chi2)/N, fmt=color+fmt, label=label)

    n    = float(sum(chi2))*db

    plt.step(    b,  chi2/n, where='mid', color=color, ls='-')

    plt.errorbar(b,  chi2/n, yerr=np.sqrt(chi2)/n, fmt=color+fmt, label=label)

def _plot_chi2_theory(bins, ndf, fmt='.', color='k', label=None):

def _plot_chi2_theory(bins, ndf,color='k', label=None):

    b    = (bins[1:]+bins[:-1])/2.0 # bin center

    xchi2 =  stats.chi2.pdf(b, ndf)

    plt.plot(b, xchi2, color=color, ls=':', lw=1, label=label)

def plot_chi2(filename, chimax=100.0, nbins=101, ndf=19, full_bw=True):

    """plot chi2"""

    d = _load_sqt_bin(filename, full_bw)

    label = ''

    if full_bw:

    rchi2, _ = np.histogram(d[:,-3], bins=bins, density=norm)

    schi2, _ = np.histogram(d[:,-2], bins=bins, density=norm)

    _plot_chi2_histo(bins, rchi2, color='r', fmt='o', label='resolution %s' % label)

    _plot_chi2_histo(bins, schi2, color='b', fmt='o', label='signal     %s' % label)

    _plot_chi2_theory(bins, ndf=ndf, color='k', fmt='.', label=r'$\chi^2 (NDF=%d)$' % ndf)

    _plot_chi2_histo(bins, rchi2, color='r', fmt='o', label=f'resolution {label}')

    _plot_chi2_histo(bins, schi2, color='b', fmt='o', label=f'signal     {label}')

    _plot_chi2_theory(bins, ndf=ndf, color='k', label=rf'$\chi^2 (NDF={ndf})$')

    plt.grid(True, which='both')

    plt.xlabel(r'$\tau$ [ns]')

    plt.hist2d(q[:,idx1], q[:,idx2], weights=w, bins=bins, **kwargs)

    plt.colorbar()

    plt.xlabel(r'Q$_%s$ [$\AA$]' % maptype[0])

    plt.ylabel(r'Q$_%s$ [$\AA$]' % maptype[1])

    plt.xlabel(rf'Q$_{maptype[0]}$ [$\AA$]')

    plt.ylabel(rf'Q$_{maptype[1]}$ [$\AA$]')

    plt.grid(True, which='both')

#!/usr/bin/env python

"OBSOLETE PLOT PHASE PORTRAITS"

from __future__ import print_function

#pylint: disable=consider-using-f-string, too-many-locals

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

#pylint: disable=too-many-arguments, too-many-positional-arguments

#pylint: disable=nested-min-max

import os

import sys

import re

import glob

import argparse

import os, sys

import re, glob

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.ticker import MultipleLocator, LinearLocator

from matplotlib.ticker import LinearLocator

_FLOAT    = r'[-+]?\d+\.?\d*(E[-+]?\d+)?'

_RE_FLOAT = re.compile(r'(%s|nan|infinity)' % _FLOAT, re.IGNORECASE)

def read_data(fname):

    meta = dict()

    data = None

    "read fits data"

    meta, data = {}, None

    if not os.path.exists(fname):

        return meta, data

    return meta, data

def save_fig(datadir, prefix="", suffix=""):

    "save fig"

    datadir  = os.path.normpath(datadir)

    dirname  = os.path.dirname(datadir)

    basename = os.path.basename(datadir)

def plot_fit_pixel(itau, pix, axes=None, datadir='.', **kwargs):

    "plot fit pixels"

    #miny=None, maxy=None,

    #xticks=True, yticks=True, mark_pixels=False):

    miny   = kwargs.get('miny', None)

    mfit, fit = read_data(fname)

    if dat is None:

        return

        return {}

    if echo_only:

        dat = dat[:,:nud]

    amp   = mfit.get('amp', 0)

    sym   = mfit.get('sym', 0) or mfit.get('pha', 0)

    try: # matplotlib 2.0+

    bgcolor = axes.set_facecolor

    except: # older matplotlib

        bgcolor = axes.set_axis_bgcolor

    if mark_pixels:

        if fit is None:

            bgcolor('gray')

    return dict(axes=axes, title=title, tau=tau, limits=(miny, maxy))

def plotspect(itau, pix=0, subpl=(211,212), datadir='.'):

    "plot spectrum"

    fname = os.path.join(datadir,"%s.spect.%s" % (itau,pix))

    print("loading %s" % fname, end="")

    dat   = np.loadtxt(fname, skiprows=10, dtype='float').T

    print('done')

    cols, rows = dat.shape

    _cols, rows = dat.shape

    tbins = range(1,rows+1)

    lbins = dat[5]*1e10      #  Angstroems

    npix = kwargs.pop('npix', 8)

    scale = kwargs.pop('scale', 1.0)

    rows, cols = npix, npix

    #rows, cols = npix, npix

    center = (npix+1)*npix//2

    tit=''

        icol = pix // npix

        fname = os.path.join(datadir,"%s.firef.%s" % (itau,pix+1)) # fit

        meta, data = read_data(fname)

        meta, _data = read_data(fname)

        if pix==center and meta and meta.get('title'):

            tit = meta.get('title')

        d[irow,icol] = meta.get(what, None)

    """

    if not os.path.isdir(datadir):

        raise RuntimeError("directory %s does not exist" % datadir)

        return

    print("plotting fits tau=%-s" % tau, end="")

    rows, cols = npix, npix

    koff  = cols - rows + 1

    nplot = rows*cols

    #rows, cols = npix, npix

    #koff  = cols - rows + 1

    #nplot = rows*cols

    center = (npix+1)*npix//2

    miny, maxy = limy

    fig, axes = plt.subplots(npix,npix, figsize=(8,8))

    _fig, axes = plt.subplots(npix,npix, figsize=(8,8))

    plt.subplots_adjust(hspace=0.05, wspace=0.05)

    tit=''

    """

    if not os.path.isdir(datadir):

        raise RuntimeError("directory %s does not exist" % datadir)

        return

    print("plotting maps tau=%-s" % tau, end="")

    minz, maxz = limz

    fig, ax = plt.subplots()

    _fig, ax = plt.subplots()

    t, img = plot_fit_map(tau, ax, datadir=datadir, what=what, scale=scale, npix=npix)

    tit = t.split()[0]

def find_minmax(datadir, tau=None):

    "find min and max"

    if tau:

        datafiles = glob.glob(datadir+'/%s.daref.*' % tau)

    else:

    miny = np.finfo(np.float).max

    maxy = np.finfo(np.float).min

    for filename in datafiles:

        meta, data = read_data(filename)

        _meta, data = read_data(filename)

        miny = min(miny, min(data[1]))

        maxy = max(maxy, max(data[1]))

    return miny, maxy

def get_info(datadir):

    # fixme

    "GET INFO"

    datafiles = [ os.path.basename(_f) for _f in glob.glob(datadir+'/*.daref.*') ]

    taus = [ int(_f.split('.')[0]) for _f in datafiles ]

    pixs = [ int(_f.split('.')[2]) for _f in datafiles ]

def main():

    import argparse

    "the main"

    parser = argparse.ArgumentParser(description='plot fits/maps')

    parser.set_defaults(datadir='.', tau=[],

                        npix=8, miny=0.0, maxy=0.0,

    args = parser.parse_args()

    taus, pixs =  get_info(args.datadir)

    taus, _pixs =  get_info(args.datadir)

    if not taus:

        parser.error("directory '%s' does not have any fit files" % args.datadir)

        else:

            plot_maps(datadir=args.datadir,

                      tau=tau,

                      mark_pixels=args.mark_pixels,

                      #mark_pixels=args.mark_pixels,

                      npix=args.npix,

                      what=args.plot,