<filename>examples/proc_linewidth.py
"""EPR linewidth
================
This example uses symfit to individually fit several Voigt-shaped ESR spectra

This example does work.
"""
from pyspecdata import *
from pyspecProcScripts import *
from pyspecProcScripts import lookup_table
import symfit as s
import pickle,os
from pylab import ndarray
from symfit import Parameter, Variable, parameters, variables, Fit, Model
from symfit.core.minimizers import MINPACK
from symfit.contrib.interactive_guess import InteractiveGuess
import numpy as np
from itertools import cycle
from scipy.optimize import nnls
init_logging(level='debug')
fl=figlist_var()
B_center = Parameter('B_center', value=-0.2)
sigma = Parameter('sigma', value=3)
R = Parameter('R', value=7.0)
A = Parameter('A', value=3e4)
C = Parameter('C', value=0)
B = Variable('B')
y_var = Variable('y')
conc_list = []
R_list = []
sigma_list = []
data = []
A_list = []
B_center_list = []
save_figure=False
for searchstr,exp_type,postproc,thisguess,interactive,concentration in [
        ("210218_1mM_TEMPOL_2",'ESR','ESR_linewidth',
            {A:2.107721e02,
                B_center: 2.391618e-01,
                R: 5.076452e-01,
                sigma: 6.672132e-01,
                },
            True,
            1e-3
            ),
        #("210114_5mM_4AT",'ESR','ESR_linewidth',
        #    {A:2.862171e02,
        #        B_center: 1.316801e-01,
        #        R: 7.333833e-01,
        #        sigma: 5.831063e-01,
        #        },
        #    False,
        #    5e-3
        #    ),
        #("210114_7mM_4AT",'ESR','ESR_linewidth',
        #    {A: 3.527704e02,
        #        B_center: 9.543103e-03,
        #        R: 8.813178e-01,
        #        sigma: 5.54451e-01,
        #        },
        #    False,
        #    7e-3
        #    ),
        #("210114_10mM_4AT",'ESR','ESR_linewidth',
        #    { # here, I entered based on the next, and then copied and pasted the result
        #        A: 4.020849e02,
        #        B_center: 8.16337e-02,
        #        R: 1.009644,
        #        sigma: 5.880406e-01,
        #        },
        #    False,
        #    10e-3
        #    ),
        #("201118_1mM4AT",'ESR','ESR_linewidth',
        #    {A: 4e2,
        ##        B_center:  4.220642e-01,
        #        R:         3.617994e-01,
        #        sigma:     8.702270e-01,
        #        },
        #    False,
        #    1e-3
        #    )
        ]:
    d = find_file(searchstr + '.DSC', exp_type=exp_type, postproc=postproc,
                  lookup=lookup_table)
    plt.figure()
    plt.title('linewidth for 1mM TEMPOL')
    plot(d)
    d = d['$B_0$':(-9, 9)]
    plot(d, '--', alpha=0.5, linewidth=4)
    d.setaxis('$B_0$', lambda x: x+1) # for a positive B_center, b/c the interactive guess doesn't deal well with negative parameters
    s_integral =d.C.run_nopop(np.cumsum, '$B_0$')
    #{{{fitting with voigt
    if not os.path.exists('dVoigt.pickle'):
        with open('dVoigt.pickle','wb') as fp:
            # cache the expression, which takes some time to generate
            print("no pickle file found -- generating")
            z = ((B-B_center) + s.I*R)/sigma/s.sqrt(2)
            faddeeva = s.simplify(s.exp(-z**2) * s.erfc(-s.I*z))
            voigt = A*s.re(faddeeva)/sigma/s.sqrt(2*s.pi)
            voigt *= sigma * R # so adjusting linewidth doesn't change amplitude
            voigt = voigt.simplify()
            # add real below b/c lambdify was giving complex answer
            dVoigt = s.re(s.re(voigt.diff(B)).simplify())
            pickle.dump(dVoigt,fp)
    else:
        with open('dVoigt.pickle','rb') as fp:
            print("reading expression from pickle")
            dVoigt = pickle.load(fp)
    plt.figure()
    plt.title('plot guess')
    logger.info(strm(A.value,"A value"))
    # {{{ need to re-do b/c defaults are stored in pickle
    for k,v in thisguess.items():
        k.value = v
    # }}}
    model_lambda = s.lambdify([B],dVoigt.subs({
        B_center:B_center.value,
        R:R.value,
        A:A.value,
        sigma:sigma.value}),
        modules=[{'ImmutableMatrix': ndarray}, 'numpy', 'scipy'])
    x_finer = r_[d.getaxis('$B_0$')[0]:d.getaxis('$B_0$')[-1]:500j]
    result = model_lambda(x_finer)
    guess_nddata = nddata(result, [-1], ['$B_0$']).setaxis(
            '$B_0$',x_finer).set_units('$B_0$',d.get_units('$B_0$'))
    plot(d, label='data')
    plot(guess_nddata,':', label='guess')
    model = s.Model({y_var:dVoigt})
    if interactive:
        guess = InteractiveGuess(model, y=d.data.real, B=d.getaxis('$B_0$'), n_points=500)
        guess.execute()
        logger.info(strm(guess))
    y_var = s.Variable('y')
    logger.info(strm("about to run fit"))
    fit = s.Fit(model, d.getaxis('$B_0$'), d.data.real) # really want to use minpack here, but gives "not proper array of floats
    fit_result = fit.execute()
    logger.info(strm("fit is done"))
    plt.figure()
    plt.title('data with fit')
    plot(d, '.', label='data')
    fit_nddata = nddata(
            fit.model(B=x_finer, **fit_result.params).y,
            [-1], ['$B_0$']).setaxis('$B_0$', x_finer)
    plot(fit_nddata, label='fit')
    logger.info(strm(fit_result))
    conc_list.append(concentration)
    R_list.append(fit_result.params['R'])
    sigma_list.append(fit_result.params['sigma'])
    A_list.append(fit_result.params['A'])
    B_center_list.append(fit_result.params['B_center'])
    data.append(d)
plt.figure()
plot(conc_list,R_list,'x',label='R')
plot(conc_list,sigma_list,'x',label=r'$\sigma$')
plt.legend(**dict(bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.))
plt.ylabel('concentration')
if save_figure:
    plt.savefig('RvsC.png',transparent=True)
plt.show()

<filename>examples/proc_linewidth.py
"""EPR linewidth
================
This example uses symfit to individually fit several Voigt-shaped ESR spectra

This example does work.
"""
from pyspecdata import *
from pyspecProcScripts import *
from pyspecProcScripts import lookup_table
import symfit as s
import pickle,os
from pylab import ndarray
from symfit import Parameter, Variable, parameters, variables, Fit, Model
from symfit.core.minimizers import MINPACK
from symfit.contrib.interactive_guess import InteractiveGuess
import numpy as np
from itertools import cycle
from scipy.optimize import nnls
init_logging(level='debug')
fl=figlist_var()
B_center = Parameter('B_center', value=-0.2)
sigma = Parameter('sigma', value=3)
R = Parameter('R', value=7.0)
A = Parameter('A', value=3e4)
C = Parameter('C', value=0)
B = Variable('B')
y_var = Variable('y')
conc_list = []
R_list = []
sigma_list = []
data = []
A_list = []
B_center_list = []
save_figure=False
for searchstr,exp_type,postproc,thisguess,interactive,concentration in [
        ("210218_1mM_TEMPOL_2",'ESR','ESR_linewidth',
            {A:2.107721e02,
                B_center: 2.391618e-01,
                R: 5.076452e-01,
                sigma: 6.672132e-01,
                },
            True,
            1e-3
            ),
        #("210114_5mM_4AT",'ESR','ESR_linewidth',
        #    {A:2.862171e02,
        #        B_center: 1.316801e-01,
        #        R: 7.333833e-01,
        #        sigma: 5.831063e-01,
        #        },
        #    False,
        #    5e-3
        #    ),
        #("210114_7mM_4AT",'ESR','ESR_linewidth',
        #    {A: 3.527704e02,
        #        B_center: 9.543103e-03,
        #        R: 8.813178e-01,
        #        sigma: 5.54451e-01,
        #        },
        #    False,
        #    7e-3
        #    ),
        #("210114_10mM_4AT",'ESR','ESR_linewidth',
        #    { # here, I entered based on the next, and then copied and pasted the result
        #        A: 4.020849e02,
        #        B_center: 8.16337e-02,
        #        R: 1.009644,
        #        sigma: 5.880406e-01,
        #        },
        #    False,
        #    10e-3
        #    ),
        #("201118_1mM4AT",'ESR','ESR_linewidth',
        #    {A: 4e2,
        ##        B_center:  4.220642e-01,
        #        R:         3.617994e-01,
        #        sigma:     8.702270e-01,
        #        },
        #    False,
        #    1e-3
        #    )
        ]:
    d = find_file(searchstr + '.DSC', exp_type=exp_type, postproc=postproc,
                  lookup=lookup_table)
    plt.figure()
    plt.title('linewidth for 1mM TEMPOL')
    plot(d)
    d = d['$B_0$':(-9, 9)]
    plot(d, '--', alpha=0.5, linewidth=4)
    d.setaxis('$B_0$', lambda x: x+1) # for a positive B_center, b/c the interactive guess doesn't deal well with negative parameters
    s_integral =d.C.run_nopop(np.cumsum, '$B_0$')
    #{{{fitting with voigt
    if not os.path.exists('dVoigt.pickle'):
        with open('dVoigt.pickle','wb') as fp:
            # cache the expression, which takes some time to generate
            print("no pickle file found -- generating")
            z = ((B-B_center) + s.I*R)/sigma/s.sqrt(2)
            faddeeva = s.simplify(s.exp(-z**2) * s.erfc(-s.I*z))
            voigt = A*s.re(faddeeva)/sigma/s.sqrt(2*s.pi)
            voigt *= sigma * R # so adjusting linewidth doesn't change amplitude
            voigt = voigt.simplify()
            # add real below b/c lambdify was giving complex answer
            dVoigt = s.re(s.re(voigt.diff(B)).simplify())
            pickle.dump(dVoigt,fp)
    else:
        with open('dVoigt.pickle','rb') as fp:
            print("reading expression from pickle")
            dVoigt = pickle.load(fp)
    plt.figure()
    plt.title('plot guess')
    logger.info(strm(A.value,"A value"))
    # {{{ need to re-do b/c defaults are stored in pickle
    for k,v in thisguess.items():
        k.value = v
    # }}}
    model_lambda = s.lambdify([B],dVoigt.subs({
        B_center:B_center.value,
        R:R.value,
        A:A.value,
        sigma:sigma.value}),
        modules=[{'ImmutableMatrix': ndarray}, 'numpy', 'scipy'])
    x_finer = r_[d.getaxis('$B_0$')[0]:d.getaxis('$B_0$')[-1]:500j]
    result = model_lambda(x_finer)
    guess_nddata = nddata(result, [-1], ['$B_0$']).setaxis(
            '$B_0$',x_finer).set_units('$B_0$',d.get_units('$B_0$'))
    plot(d, label='data')
    plot(guess_nddata,':', label='guess')
    model = s.Model({y_var:dVoigt})
    if interactive:
        guess = InteractiveGuess(model, y=d.data.real, B=d.getaxis('$B_0$'), n_points=500)
        guess.execute()
        logger.info(strm(guess))
    y_var = s.Variable('y')
    logger.info(strm("about to run fit"))
    fit = s.Fit(model, d.getaxis('$B_0$'), d.data.real) # really want to use minpack here, but gives "not proper array of floats
    fit_result = fit.execute()
    logger.info(strm("fit is done"))
    plt.figure()
    plt.title('data with fit')
    plot(d, '.', label='data')
    fit_nddata = nddata(
            fit.model(B=x_finer, **fit_result.params).y,
            [-1], ['$B_0$']).setaxis('$B_0$', x_finer)
    plot(fit_nddata, label='fit')
    logger.info(strm(fit_result))
    conc_list.append(concentration)
    R_list.append(fit_result.params['R'])
    sigma_list.append(fit_result.params['sigma'])
    A_list.append(fit_result.params['A'])
    B_center_list.append(fit_result.params['B_center'])
    data.append(d)
plt.figure()
plot(conc_list,R_list,'x',label='R')
plot(conc_list,sigma_list,'x',label=r'$\sigma$')
plt.legend(**dict(bbox_to_anchor=(1.05,1), loc=2, borderaxespad=0.))
plt.ylabel('concentration')
if save_figure:
    plt.savefig('RvsC.png',transparent=True)
plt.show()

EPR linewidth ================ This example uses symfit to individually fit several Voigt-shaped ESR spectra  This example does work. #("210114_5mM_4AT",'ESR','ESR_linewidth', #    {A:2.862171e02, #        B_center: 1.316801e-01, #        R: 7.333833e-01, #        sigma: 5.831063e-01, #        }, #    False, #    5e-3 #    ), #("210114_7mM_4AT",'ESR','ESR_linewidth', #    {A: 3.527704e02, #        B_center: 9.543103e-03, #        R: 8.813178e-01, #        sigma: 5.54451e-01, #        }, #    False, #    7e-3 #    ), #("210114_10mM_4AT",'ESR','ESR_linewidth', #    { # here, I entered based on the next, and then copied and pasted the result #        A: 4.020849e02, #        B_center: 8.16337e-02, #        R: 1.009644, #        sigma: 5.880406e-01, #        }, #    False, #    10e-3 #    ), #("201118_1mM4AT",'ESR','ESR_linewidth', #    {A: 4e2, ##        B_center:  4.220642e-01, #        R:         3.617994e-01, #        sigma:     8.702270e-01, #        }, #    False, #    1e-3 #    ) # for a positive B_center, b/c the interactive guess doesn't deal well with negative parameters #{{{fitting with voigt # cache the expression, which takes some time to generate # so adjusting linewidth doesn't change amplitude # add real below b/c lambdify was giving complex answer # {{{ need to re-do b/c defaults are stored in pickle # }}} # really want to use minpack here, but gives "not proper array of floats