Loading pysen/echo/fit.py +1 −1 Original line number Diff line number Diff line Loading @@ -273,7 +273,7 @@ def fit_echo_current(current, counts, spectrum, lam0, phase0, **kwargs): phase_pf,_ = fit_poly_current(current, counts, maximum=not incoherent) # max count dj = np.radians(phasesens*(current - phase0))/OMEGA_N/lam0 # convert currents to dJ ddj = dphase/OMEGA_N/lam0 # fit range sfun = partial(spectrum_weighted_eshape_spectrum, spectrum=spectrum) sfun = partial(spectrum_weighted_eshape, spectrum=spectrum) dj0 = np.radians(phasesens*(phase_pf - phase0))/OMEGA_N/lam0 res = echo_fit(sfun, dj, counts[0], counts[1], dj0=dj0, djlim=(dj0-ddj,dj0+ddj)) maxchi = res['chi2'] Loading pysen/echo/simulator.py +24 −25 Original line number Diff line number Diff line Loading @@ -23,22 +23,26 @@ class SimpleSimulator(): #pylint: disable=too-many-instance-attributes def __init__(self, lmax=8.0, ntbins=42, pos='p2'): """init""" self.stats = {} # nsteps = int(2*NFINE*MAX_PHI/DEL_PHI)+1 self.dj = np.linspace(-1, 1, nsteps)*MAX_PHI/(lmax*ANGSTROM*OMEGA_N) self.echo = np.zeros_like(self.dj) self.create_spectrum(lmax, ntbins, pos) self.set_intensities() self.stats = {} def create_spectrum(self, lmax=8.0, ntbins=42, pos='p2'): # "create 'simple' spectrum" # fixme - lmax is set here and in the constructor lmin = lmax - wavelength_bandwidth(pos=pos) lbin = np.linspace(lmin,lmax,ntbins+1)*ANGSTROM # neutron wavelength bins lave = (lbin[1:]+lbin[:-1])/2 # average wavelength dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda flux = (lave[0]/lave)**2 self.spectrum = Spectrum(lave, dlam, flux) # self.set_intensities() def get(self, key, default=None): "get attribute" if key in self.stats: return self.stats[key] return default def set_intensities(self, coherent=18.0, incoherent=6.0, background=0.0, tau_coh=1.0, tau_inc=1.0): "set intensities" Loading @@ -48,9 +52,14 @@ class SimpleSimulator(): self.tau_coh = tau_coh self.tau_inc = tau_inc up = coherent+1/3*incoherent dn = 2/3*incoherent fr = up/dn self.stats.update(up=up, dn=dn, fr=fr, ave=(up+dn)/2, amp=(up-dn)/2) def f_qt(self, tau): "return F(Q,t)" updn = self.coherent + 1/3*self.incoherent - 2/3*self.incoherent updn = self.coherent - 1/3*self.incoherent return fqt(tau, self.coherent, self.tau_coh, self.incoherent, self.tau_inc)/updn def f_coh(self, tau): Loading @@ -61,22 +70,12 @@ class SimpleSimulator(): "return F_incoherent(Q,t)" return f_exp(tau, self.tau_inc) def get_echo(self, tau): def echo(self, tau): "get echo" coh = self.coherent inc = self.incoherent up = coh+1/3*inc dn = 2/3*inc fr = up/dn ave = (up+dn)/2 amp = (up-dn)/2 sqt = self.f_qt(tau)/(up-dn) self.echo = echo_shape(self.dj, self.spectrum)*amp + ave self.stats = {'amp': amp, 'ave': ave, 'up': up, 'dn': dn, 'fr': fr, 'sqt': sqt} return (self.dj, self.echo) amp = self.get('amp',-1) ave = self.get('ave',-1) # yecho = abs(self.f_qt(tau))*echo_shape(self.dj, self.spectrum)*amp + ave return (self.dj, yecho) #EOF pysen/revision.py +2 −2 Original line number Diff line number Diff line Loading @@ -3,8 +3,8 @@ PySEN revision module """ import sys __version__ = "2.1" __release__ = "b3" __date__ = "Mar 31, 2025" __release__ = "b4" __date__ = "Apr 1, 2025" def version(full=False): "get pysen version number" Loading pysen/ui/EchoSimWidget.py +44 −62 Original line number Diff line number Diff line Loading @@ -5,26 +5,16 @@ import numpy as np # from qtpy import QtGui # from numpy import pi, exp, log10 from numpy import exp, log10 # from matplotlib.figure import Figure from matplotlib.backends.backend_qtagg import FigureCanvas #pylint: disable=no-name-in-module from ..config import wavelength_bandwidth from ..constants import ANGSTROM, OMEGA_N from ..echo.fit import flux_weighted_eshape as echo_shape from ..constants import MICRO from ..echo.simulator import SimpleSimulator, NFINE # from .ui_EchoSimWidget import Ui_EchoSimWidget MICRO=1e-6 def f_exp(t, t0=1): "stretched exponential" return exp(-(t/t0)) def sqt(t, coh, tcoh, inc, tinc, bgr=0): "S(Q,t)/S(Q,0)" return coh*f_exp(t,tcoh)-1/3*(inc*f_exp(t,tinc)+bgr) class EchoSimWidget(Ui_EchoSimWidget): """Main for viewing and deleting available groups.""" Loading @@ -32,18 +22,12 @@ class EchoSimWidget(Ui_EchoSimWidget): def __init__(self, parent=None, **kwargs): """init MainWindow""" super().__init__() self.pos = kwargs.pop('pos', 'p2') self.lmax = kwargs.pop('lmax', 8.000) self.lmin = self.lmax - wavelength_bandwidth(pos=self.pos) self.ntbin = kwargs.pop('ntbin', 42) # number of TOF bins pos = kwargs.pop('pos', 'p2') lmax = kwargs.pop('lmax', 8.000) ntbins = kwargs.pop('ntbins', 42) # number of TOF bins self.esim = SimpleSimulator(lmax=lmax, ntbins=ntbins, pos=pos) self.tmin, self.tmax = 1e-3, 1e+3 # self.lave = None self.dlam = None self.flux = None self.dj = None self.dj2 = None self.xt = None self.taus = np.logspace(log10(self.tmin), log10(self.tmax), 201) # for S(Q,t) self.lines = {} # self.initUi(parent) Loading Loading @@ -137,15 +121,16 @@ class EchoSimWidget(Ui_EchoSimWidget): def _init_canvas(self): "initialize canvas" lbin = np.linspace(self.lmin,self.lmax,self.ntbin+1)*ANGSTROM # neutron wavelength bins xt = np.logspace(log10(self.tmin), log10(self.tmax), 201) # for S(Q,t) djmax = 3*pi/(self.lmax*ANGSTROM)/OMEGA_N # field integral range (-3pi,+3pi) self.lave = (lbin[1:]+lbin[:-1])/2 # average wavelength self.dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda self.flux = (self.lave[0]/self.lave)**2 self.dj = np.linspace(-djmax, djmax, 201) # for plotting dj points self.dj2 = np.linspace(-djmax, djmax, 11) # for plotting dj points self.xt = xt #lbin = np.linspace(self.lmin,self.lmax,self.ntbin+1)*ANGSTROM # neutron wavelength bins #djmax = 3*pi/(self.lmax*ANGSTROM)/OMEGA_N # field integral range (-3pi,+3pi) #self.lave = (lbin[1:]+lbin[:-1])/2 # average wavelength #self.dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda #self.flux = (self.lave[0]/self.lave)**2 #self.dj = np.linspace(-djmax, djmax, 201) # for plotting dj points #self.dj2 = np.linspace(-djmax, djmax, 11) # for plotting dj points #self.xt = xt xt = self.taus ax = self._ax[0] self.lines = {} Loading @@ -162,11 +147,11 @@ class EchoSimWidget(Ui_EchoSimWidget): ax.grid(True) #, which='both') ax = self._ax[1] self.lines['echo'], = ax.plot(self.dj/MICRO, np.zeros_like(self.dj), 'k-', lw=2) self.lines['echo2'], = ax.plot(self.dj2/MICRO, np.zeros_like(self.dj2),'ro', lw=2) self.lines['echo'], = ax.plot(self.esim.dj/MICRO, np.zeros_like(self.esim.dj),'g-', lw=3) self.lines['epts'], = ax.plot(self.esim.dj/MICRO, np.zeros_like(self.esim.dj),'ro', ms=5) self.lines['up'] = ax.axhline(2.0, lw=1, ls='--', color='red') self.lines['dn'] = ax.axhline(1.0, lw=1, ls='--', color='blue') self.lines['ave'] = ax.axhline(1.5, lw=2, ls='-', color='k') self.lines['ave'] = ax.axhline(1.5, lw=1, ls='--', color='k') ax.set_xlabel(r'$\Delta J$ [$\mu Tm$]') ax.yaxis.tick_right() ax.grid(True) Loading @@ -184,16 +169,22 @@ class EchoSimWidget(Ui_EchoSimWidget): techo = 10**(self.hSliderTauEcho.value()/1000.0) ax = self._ax[1] try: self.esim.set_intensities(coherent=coh, incoherent=inc, background=bgr, tau_coh=tcoh, tau_inc=tinc) except ArithmeticError as exc: print(exc) ax.set_title(f'!!! {exc} !!!') ax.figure.canvas.draw() return # check inputs if inc<0 or coh<0: ax.set_title('!!! Coherent or Incoherent cannot be negative !!!') ax.figure.canvas.draw() return dn = 2/3*(inc+bgr) up = coh+1/3*(inc+bgr) if up>0 and dn>0: fratio = up/dn else: dn = self.esim.get('up', -1) # 2/3*(inc+bgr) up = self.esim.get('dn', -1) # coh+1/3*(inc+bgr) if up<0 or dn<0: ax.set_title('!!! Up and Down counts must be positive') ax.figure.canvas.draw() return Loading @@ -201,41 +192,32 @@ class EchoSimWidget(Ui_EchoSimWidget): ax.set_title('!!! Invalid input') ax.figure.canvas.draw() return fr = self.esim.get('fr', -1) # up/dn ave= self.esim.get('ave', -1) # (up+dn)/2 ax.set_title(rf'FR={fr:.3g} (U={up:.0f}/D={dn:.0f})') ax.set_title(rf'FR={fratio:.3g} (U={up:.0f}/D={dn:.0f})') ave = (up+dn)/2 amp = (up-dn)/2*sqt(techo,coh,tcoh,inc,tinc,bgr)/abs(up-dn) dj = self.dj lave = self.lave dlam = self.dlam flux = self.flux xt = self.xt yecho = echo_shape(dj, lave, dlam,flux)*amp + ave dj, yecho = self.esim.echo(techo) self.lines['echo'].set_data(dj/MICRO, yecho) yecho2 = echo_shape(self.dj2, lave, dlam,flux)*amp + ave self.lines['echo2'].set_data(self.dj2/MICRO, yecho2) self.lines['epts'].set_data(dj[::NFINE]/MICRO, yecho[::NFINE]) self.lines['up'].set_data(dj/MICRO,up*np.ones_like(dj)) self.lines['dn'].set_data(dj/MICRO,dn*np.ones_like(dj)) self.lines['ave'].set_data(dj/MICRO,ave*np.ones_like(dj)) #bot = min(np.amin(yecho), up, dn)*0.67 #top = max(np.amax(yecho), up, dn)*1.2 ax.set_ylim(bottom=None, top=None) bot = min(np.amin(yecho), up, dn)*0.67 top = max(np.amax(yecho), up, dn)*1.2 ax.set_ylim(bottom=bot, top=top) ysqt = sqt(xt, coh, tcoh, inc, tinc)/(up-dn) ycoh = f_exp(xt, tcoh) yinc = f_exp(xt, tinc) xt = self.taus ysqt = self.esim.f_qt(xt) ycoh = self.esim.f_coh(xt) yinc = self.esim.f_inc(xt) self.lines['sqt'].set_data(xt, ysqt) self.lines['coh'].set_data(xt, ycoh) self.lines['inc'].set_data(xt, yinc) ax = self._ax[0] ymin, ymax = -abs(np.amin([ysqt,ycoh,yinc])), np.amax([ysqt,ycoh,yinc]) self.lines['tau'].set_data([techo,techo],[ymin,ymax]) Loading test/test_echo_fit.py +75 −28 Original line number Diff line number Diff line Loading @@ -10,10 +10,11 @@ from functools import partial import numpy as np from pysen import ANGSTROM, MICRO from pysen import ANGSTROM, MICRO, OMEGA_N from pysen.echo.fit import ( echo_fit, echo_curve, echo_fit_4point, linear_fit, flux_weighted_eshape, spectrum_weighted_eshape, phase_fit_local, phase_fit_global, Spectrum ) phase_fit_local, phase_fit_global, Spectrum, fit_echo_current) TestDataDir = os.path.join(os.path.dirname(__file__), 'data') Loading @@ -35,7 +36,8 @@ def _read_echo_data(fname): flux = data[:,27:] return flux, echo def _phase_sens(lam0): return 198*lam0/ANGSTROM class EchoFitTestCase(unittest.TestCase): "Echo Fit Test Cases" Loading Loading @@ -183,9 +185,30 @@ class EchoFitTestCase(unittest.TestCase): #print("%s %g %g %g" % (key, v4, v0, dv)) self.assertTrue(dv<maxdv) def plot_example(fname, axes, dj0=None, sgn=1): def test_fit_echo_current(self): "fit echo current" avlam, delam, flux = self.flux dj, counts, ecounts = self.echo # phase0 = 2.0 lam0 = avlam[-1] phasesens = _phase_sens(lam0) # currents = np.degrees(dj*lam0*OMEGA_N)/phasesens+phase0 spectrum = Spectrum(lam=avlam, dlam=delam, flux=flux) res = fit_echo_current(currents, (counts, ecounts), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) (phase_ef, phase_err), _ , res = res print() print(lam0, phasesens) print(res) self.assertAlmostEqual( 1.994 , phase_ef , 3) self.assertAlmostEqual( 0.006 , phase_err , 3) self.assertAlmostEqual( 8.498 , res['chi2'], 3) self.assertAlmostEqual( 17.0 , res['ndf'] , 3) def plot_example1(fname, axes, dj0=None, sgn=1): "test function" (avlam, delam, flux), (dj, counts, ecounts) = _read_echo_data(fname) sfun = partial(flux_weighted_eshape, avlam=avlam, delam=delam, flux=flux) Loading Loading @@ -213,13 +236,7 @@ def plot_example(fname, axes, dj0=None, sgn=1): amplitude=res['amplitude'][0], average=res['average'][0], npoints=1001) # for measured point chi2 = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in dj ]) # smooth curve chif = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in djf ]) #chi2 = chi2/res['ndf'] #chif = chif/res['ndf'] chi2 = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in djf ]) # ============================================================================================= #plot test results Loading @@ -236,7 +253,7 @@ def plot_example(fname, axes, dj0=None, sgn=1): dj0 = res['dj0'][0] edj0 = res['dj0'][1] ax2a.errorbar(dj/MICRO, counts, yerr=ecounts, fmt='ko', ms=5) ax2a.plot(djf/MICRO, efit, 'r-') ax2a.plot(djf/MICRO, efit, 'r-', label='echo_fit') if init_dj0 is not None: # initial phase ax2a.axvline(init_dj0/MICRO, color='k', lw=1.0, ls=':') # # Loading @@ -249,41 +266,71 @@ def plot_example(fname, axes, dj0=None, sgn=1): ax2a.axhline(res['average'][0]+res['amplitude'][0], color='g', lw=1.0, ls='--') # average + ampl # plot chi^2 ax2b.plot(djf/MICRO, chif, 'b--', lw=1.0) ax2b.plot(dj /MICRO, chi2, 'bo' , ms=2.0) ax2b.plot(djf/MICRO, chi2, 'b--', lw=1.0) # plot 4pt method ax2a.errorbar(dj[idx_pt]/MICRO, counts[idx_pt], fmt='ms', ms=3) # no error bars ax2a.plot(djf4/MICRO, efit4, 'm--', lw=0.5) ax2a.plot(djf4/MICRO, efit4, 'm--', lw=1, label='4pt_fit') ax2a.axvline(dj0_pt/MICRO , color='m' , lw=0.5 ) # phase fit 4pt ax2a.set_ylabel('Counts [arb.]') ax2a.set_xlabel(r'$\Delta J$ [$\mu$ Tm]') #ax2b.set_ylabel(r'$\chi^2$/n.d.f.', color='b', ) ax2b.set_ylabel(r'$\chi^2$', color='b', ) ax2a.grid(True, lw=0.5, color='k', ls='--') #ax2b.grid(True, lw=0.5, color='b', ls=':') ax2a.legend(loc=0) def plot_example2(fname, axis): "test function" (avlam, delam, flux), (dj, counts, ecounts) = _read_echo_data(fname) lam0 = avlam[-1] phase0 = 2.0 phasesens = _phase_sens(lam0) # currents = np.degrees(dj*lam0*OMEGA_N)/phasesens+phase0 spectrum = Spectrum(lam=avlam, dlam=delam, flux=flux) res = fit_echo_current(currents, (counts, ecounts), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) (phase_ef, phase_err), (curf, efit), res = res print(f"current_fit => p0 ({phase_ef:.4f} +/- {phase_err:.4f})") axis.errorbar(currents, counts, yerr=ecounts, fmt='ko') axis.plot(curf, efit, 'r-', lw=2, label='fit_echo_current') # axis.axvline(phase_ef-phase_err, color='r' , lw=0.5, ls='--' ) # axis.axvline(phase_ef , color='r' , lw=1.5, ls='-' ) # phase fit axis.axvline(phase_ef+phase_err, color='r' , lw=0.5, ls='--' ) # axis.set_xlabel('Phase Current [A]') axis.set_ylabel('Counts [arb.]') axis.grid(True, lw=0.5, color='k', ls='--') axis.legend(loc=0) if __name__ == "__main__": unittest.main(exit=False, verbosity=2) # ================================================================== # graphical tests import matplotlib.pyplot as plt #filename = os.path.join(TestDataDir, 'echo001.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo002.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo003.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo004.dat') #- 4.0-7.0 - P2, incoherent filename = os.path.join(TestDataDir, 'echo005.dat') #- 5.6-8.0 - P4 #filename = os.path.join(TestDataDir, 'echo006.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo005.dat') #- 5.6-8.0 - P4 filename = os.path.join(TestDataDir, 'echo006.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo007.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo008.dat') #- 7.0-8.0 - P2, partial spectrum unittest.main(exit=False, verbosity=1) signum = 0.0 if len(sys.argv)>1: filename=sys.argv[1] # graphical tests import matplotlib.pyplot as plt if len(sys.argv)>2: signum = np.sign(float(sys.argv[2])) else: signum = 0.0 _, axs = plt.subplots(2, 1, figsize=(8,8)) plot_example(filename, axs, dj0=0, sgn=signum) if not os.path.exists(filename): sys.exit("File not found: %s" % filename) _, axs = plt.subplots(3, 1, figsize=(6,8)) plot_example1(filename, axs[:-1], dj0=0, sgn=signum) plot_example2(filename, axs[-1]) plt.tight_layout() plt.show() Loading
pysen/echo/fit.py +1 −1 Original line number Diff line number Diff line Loading @@ -273,7 +273,7 @@ def fit_echo_current(current, counts, spectrum, lam0, phase0, **kwargs): phase_pf,_ = fit_poly_current(current, counts, maximum=not incoherent) # max count dj = np.radians(phasesens*(current - phase0))/OMEGA_N/lam0 # convert currents to dJ ddj = dphase/OMEGA_N/lam0 # fit range sfun = partial(spectrum_weighted_eshape_spectrum, spectrum=spectrum) sfun = partial(spectrum_weighted_eshape, spectrum=spectrum) dj0 = np.radians(phasesens*(phase_pf - phase0))/OMEGA_N/lam0 res = echo_fit(sfun, dj, counts[0], counts[1], dj0=dj0, djlim=(dj0-ddj,dj0+ddj)) maxchi = res['chi2'] Loading
pysen/echo/simulator.py +24 −25 Original line number Diff line number Diff line Loading @@ -23,22 +23,26 @@ class SimpleSimulator(): #pylint: disable=too-many-instance-attributes def __init__(self, lmax=8.0, ntbins=42, pos='p2'): """init""" self.stats = {} # nsteps = int(2*NFINE*MAX_PHI/DEL_PHI)+1 self.dj = np.linspace(-1, 1, nsteps)*MAX_PHI/(lmax*ANGSTROM*OMEGA_N) self.echo = np.zeros_like(self.dj) self.create_spectrum(lmax, ntbins, pos) self.set_intensities() self.stats = {} def create_spectrum(self, lmax=8.0, ntbins=42, pos='p2'): # "create 'simple' spectrum" # fixme - lmax is set here and in the constructor lmin = lmax - wavelength_bandwidth(pos=pos) lbin = np.linspace(lmin,lmax,ntbins+1)*ANGSTROM # neutron wavelength bins lave = (lbin[1:]+lbin[:-1])/2 # average wavelength dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda flux = (lave[0]/lave)**2 self.spectrum = Spectrum(lave, dlam, flux) # self.set_intensities() def get(self, key, default=None): "get attribute" if key in self.stats: return self.stats[key] return default def set_intensities(self, coherent=18.0, incoherent=6.0, background=0.0, tau_coh=1.0, tau_inc=1.0): "set intensities" Loading @@ -48,9 +52,14 @@ class SimpleSimulator(): self.tau_coh = tau_coh self.tau_inc = tau_inc up = coherent+1/3*incoherent dn = 2/3*incoherent fr = up/dn self.stats.update(up=up, dn=dn, fr=fr, ave=(up+dn)/2, amp=(up-dn)/2) def f_qt(self, tau): "return F(Q,t)" updn = self.coherent + 1/3*self.incoherent - 2/3*self.incoherent updn = self.coherent - 1/3*self.incoherent return fqt(tau, self.coherent, self.tau_coh, self.incoherent, self.tau_inc)/updn def f_coh(self, tau): Loading @@ -61,22 +70,12 @@ class SimpleSimulator(): "return F_incoherent(Q,t)" return f_exp(tau, self.tau_inc) def get_echo(self, tau): def echo(self, tau): "get echo" coh = self.coherent inc = self.incoherent up = coh+1/3*inc dn = 2/3*inc fr = up/dn ave = (up+dn)/2 amp = (up-dn)/2 sqt = self.f_qt(tau)/(up-dn) self.echo = echo_shape(self.dj, self.spectrum)*amp + ave self.stats = {'amp': amp, 'ave': ave, 'up': up, 'dn': dn, 'fr': fr, 'sqt': sqt} return (self.dj, self.echo) amp = self.get('amp',-1) ave = self.get('ave',-1) # yecho = abs(self.f_qt(tau))*echo_shape(self.dj, self.spectrum)*amp + ave return (self.dj, yecho) #EOF
pysen/revision.py +2 −2 Original line number Diff line number Diff line Loading @@ -3,8 +3,8 @@ PySEN revision module """ import sys __version__ = "2.1" __release__ = "b3" __date__ = "Mar 31, 2025" __release__ = "b4" __date__ = "Apr 1, 2025" def version(full=False): "get pysen version number" Loading
pysen/ui/EchoSimWidget.py +44 −62 Original line number Diff line number Diff line Loading @@ -5,26 +5,16 @@ import numpy as np # from qtpy import QtGui # from numpy import pi, exp, log10 from numpy import exp, log10 # from matplotlib.figure import Figure from matplotlib.backends.backend_qtagg import FigureCanvas #pylint: disable=no-name-in-module from ..config import wavelength_bandwidth from ..constants import ANGSTROM, OMEGA_N from ..echo.fit import flux_weighted_eshape as echo_shape from ..constants import MICRO from ..echo.simulator import SimpleSimulator, NFINE # from .ui_EchoSimWidget import Ui_EchoSimWidget MICRO=1e-6 def f_exp(t, t0=1): "stretched exponential" return exp(-(t/t0)) def sqt(t, coh, tcoh, inc, tinc, bgr=0): "S(Q,t)/S(Q,0)" return coh*f_exp(t,tcoh)-1/3*(inc*f_exp(t,tinc)+bgr) class EchoSimWidget(Ui_EchoSimWidget): """Main for viewing and deleting available groups.""" Loading @@ -32,18 +22,12 @@ class EchoSimWidget(Ui_EchoSimWidget): def __init__(self, parent=None, **kwargs): """init MainWindow""" super().__init__() self.pos = kwargs.pop('pos', 'p2') self.lmax = kwargs.pop('lmax', 8.000) self.lmin = self.lmax - wavelength_bandwidth(pos=self.pos) self.ntbin = kwargs.pop('ntbin', 42) # number of TOF bins pos = kwargs.pop('pos', 'p2') lmax = kwargs.pop('lmax', 8.000) ntbins = kwargs.pop('ntbins', 42) # number of TOF bins self.esim = SimpleSimulator(lmax=lmax, ntbins=ntbins, pos=pos) self.tmin, self.tmax = 1e-3, 1e+3 # self.lave = None self.dlam = None self.flux = None self.dj = None self.dj2 = None self.xt = None self.taus = np.logspace(log10(self.tmin), log10(self.tmax), 201) # for S(Q,t) self.lines = {} # self.initUi(parent) Loading Loading @@ -137,15 +121,16 @@ class EchoSimWidget(Ui_EchoSimWidget): def _init_canvas(self): "initialize canvas" lbin = np.linspace(self.lmin,self.lmax,self.ntbin+1)*ANGSTROM # neutron wavelength bins xt = np.logspace(log10(self.tmin), log10(self.tmax), 201) # for S(Q,t) djmax = 3*pi/(self.lmax*ANGSTROM)/OMEGA_N # field integral range (-3pi,+3pi) self.lave = (lbin[1:]+lbin[:-1])/2 # average wavelength self.dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda self.flux = (self.lave[0]/self.lave)**2 self.dj = np.linspace(-djmax, djmax, 201) # for plotting dj points self.dj2 = np.linspace(-djmax, djmax, 11) # for plotting dj points self.xt = xt #lbin = np.linspace(self.lmin,self.lmax,self.ntbin+1)*ANGSTROM # neutron wavelength bins #djmax = 3*pi/(self.lmax*ANGSTROM)/OMEGA_N # field integral range (-3pi,+3pi) #self.lave = (lbin[1:]+lbin[:-1])/2 # average wavelength #self.dlam = (lbin[1:]-lbin[:-1])/2 # delta-lambda #self.flux = (self.lave[0]/self.lave)**2 #self.dj = np.linspace(-djmax, djmax, 201) # for plotting dj points #self.dj2 = np.linspace(-djmax, djmax, 11) # for plotting dj points #self.xt = xt xt = self.taus ax = self._ax[0] self.lines = {} Loading @@ -162,11 +147,11 @@ class EchoSimWidget(Ui_EchoSimWidget): ax.grid(True) #, which='both') ax = self._ax[1] self.lines['echo'], = ax.plot(self.dj/MICRO, np.zeros_like(self.dj), 'k-', lw=2) self.lines['echo2'], = ax.plot(self.dj2/MICRO, np.zeros_like(self.dj2),'ro', lw=2) self.lines['echo'], = ax.plot(self.esim.dj/MICRO, np.zeros_like(self.esim.dj),'g-', lw=3) self.lines['epts'], = ax.plot(self.esim.dj/MICRO, np.zeros_like(self.esim.dj),'ro', ms=5) self.lines['up'] = ax.axhline(2.0, lw=1, ls='--', color='red') self.lines['dn'] = ax.axhline(1.0, lw=1, ls='--', color='blue') self.lines['ave'] = ax.axhline(1.5, lw=2, ls='-', color='k') self.lines['ave'] = ax.axhline(1.5, lw=1, ls='--', color='k') ax.set_xlabel(r'$\Delta J$ [$\mu Tm$]') ax.yaxis.tick_right() ax.grid(True) Loading @@ -184,16 +169,22 @@ class EchoSimWidget(Ui_EchoSimWidget): techo = 10**(self.hSliderTauEcho.value()/1000.0) ax = self._ax[1] try: self.esim.set_intensities(coherent=coh, incoherent=inc, background=bgr, tau_coh=tcoh, tau_inc=tinc) except ArithmeticError as exc: print(exc) ax.set_title(f'!!! {exc} !!!') ax.figure.canvas.draw() return # check inputs if inc<0 or coh<0: ax.set_title('!!! Coherent or Incoherent cannot be negative !!!') ax.figure.canvas.draw() return dn = 2/3*(inc+bgr) up = coh+1/3*(inc+bgr) if up>0 and dn>0: fratio = up/dn else: dn = self.esim.get('up', -1) # 2/3*(inc+bgr) up = self.esim.get('dn', -1) # coh+1/3*(inc+bgr) if up<0 or dn<0: ax.set_title('!!! Up and Down counts must be positive') ax.figure.canvas.draw() return Loading @@ -201,41 +192,32 @@ class EchoSimWidget(Ui_EchoSimWidget): ax.set_title('!!! Invalid input') ax.figure.canvas.draw() return fr = self.esim.get('fr', -1) # up/dn ave= self.esim.get('ave', -1) # (up+dn)/2 ax.set_title(rf'FR={fr:.3g} (U={up:.0f}/D={dn:.0f})') ax.set_title(rf'FR={fratio:.3g} (U={up:.0f}/D={dn:.0f})') ave = (up+dn)/2 amp = (up-dn)/2*sqt(techo,coh,tcoh,inc,tinc,bgr)/abs(up-dn) dj = self.dj lave = self.lave dlam = self.dlam flux = self.flux xt = self.xt yecho = echo_shape(dj, lave, dlam,flux)*amp + ave dj, yecho = self.esim.echo(techo) self.lines['echo'].set_data(dj/MICRO, yecho) yecho2 = echo_shape(self.dj2, lave, dlam,flux)*amp + ave self.lines['echo2'].set_data(self.dj2/MICRO, yecho2) self.lines['epts'].set_data(dj[::NFINE]/MICRO, yecho[::NFINE]) self.lines['up'].set_data(dj/MICRO,up*np.ones_like(dj)) self.lines['dn'].set_data(dj/MICRO,dn*np.ones_like(dj)) self.lines['ave'].set_data(dj/MICRO,ave*np.ones_like(dj)) #bot = min(np.amin(yecho), up, dn)*0.67 #top = max(np.amax(yecho), up, dn)*1.2 ax.set_ylim(bottom=None, top=None) bot = min(np.amin(yecho), up, dn)*0.67 top = max(np.amax(yecho), up, dn)*1.2 ax.set_ylim(bottom=bot, top=top) ysqt = sqt(xt, coh, tcoh, inc, tinc)/(up-dn) ycoh = f_exp(xt, tcoh) yinc = f_exp(xt, tinc) xt = self.taus ysqt = self.esim.f_qt(xt) ycoh = self.esim.f_coh(xt) yinc = self.esim.f_inc(xt) self.lines['sqt'].set_data(xt, ysqt) self.lines['coh'].set_data(xt, ycoh) self.lines['inc'].set_data(xt, yinc) ax = self._ax[0] ymin, ymax = -abs(np.amin([ysqt,ycoh,yinc])), np.amax([ysqt,ycoh,yinc]) self.lines['tau'].set_data([techo,techo],[ymin,ymax]) Loading
test/test_echo_fit.py +75 −28 Original line number Diff line number Diff line Loading @@ -10,10 +10,11 @@ from functools import partial import numpy as np from pysen import ANGSTROM, MICRO from pysen import ANGSTROM, MICRO, OMEGA_N from pysen.echo.fit import ( echo_fit, echo_curve, echo_fit_4point, linear_fit, flux_weighted_eshape, spectrum_weighted_eshape, phase_fit_local, phase_fit_global, Spectrum ) phase_fit_local, phase_fit_global, Spectrum, fit_echo_current) TestDataDir = os.path.join(os.path.dirname(__file__), 'data') Loading @@ -35,7 +36,8 @@ def _read_echo_data(fname): flux = data[:,27:] return flux, echo def _phase_sens(lam0): return 198*lam0/ANGSTROM class EchoFitTestCase(unittest.TestCase): "Echo Fit Test Cases" Loading Loading @@ -183,9 +185,30 @@ class EchoFitTestCase(unittest.TestCase): #print("%s %g %g %g" % (key, v4, v0, dv)) self.assertTrue(dv<maxdv) def plot_example(fname, axes, dj0=None, sgn=1): def test_fit_echo_current(self): "fit echo current" avlam, delam, flux = self.flux dj, counts, ecounts = self.echo # phase0 = 2.0 lam0 = avlam[-1] phasesens = _phase_sens(lam0) # currents = np.degrees(dj*lam0*OMEGA_N)/phasesens+phase0 spectrum = Spectrum(lam=avlam, dlam=delam, flux=flux) res = fit_echo_current(currents, (counts, ecounts), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) (phase_ef, phase_err), _ , res = res print() print(lam0, phasesens) print(res) self.assertAlmostEqual( 1.994 , phase_ef , 3) self.assertAlmostEqual( 0.006 , phase_err , 3) self.assertAlmostEqual( 8.498 , res['chi2'], 3) self.assertAlmostEqual( 17.0 , res['ndf'] , 3) def plot_example1(fname, axes, dj0=None, sgn=1): "test function" (avlam, delam, flux), (dj, counts, ecounts) = _read_echo_data(fname) sfun = partial(flux_weighted_eshape, avlam=avlam, delam=delam, flux=flux) Loading Loading @@ -213,13 +236,7 @@ def plot_example(fname, axes, dj0=None, sgn=1): amplitude=res['amplitude'][0], average=res['average'][0], npoints=1001) # for measured point chi2 = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in dj ]) # smooth curve chif = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in djf ]) #chi2 = chi2/res['ndf'] #chif = chif/res['ndf'] chi2 = np.asarray([ linear_fit(sfun(dj-_dj0), counts, ecounts, chionly=True) for _dj0 in djf ]) # ============================================================================================= #plot test results Loading @@ -236,7 +253,7 @@ def plot_example(fname, axes, dj0=None, sgn=1): dj0 = res['dj0'][0] edj0 = res['dj0'][1] ax2a.errorbar(dj/MICRO, counts, yerr=ecounts, fmt='ko', ms=5) ax2a.plot(djf/MICRO, efit, 'r-') ax2a.plot(djf/MICRO, efit, 'r-', label='echo_fit') if init_dj0 is not None: # initial phase ax2a.axvline(init_dj0/MICRO, color='k', lw=1.0, ls=':') # # Loading @@ -249,41 +266,71 @@ def plot_example(fname, axes, dj0=None, sgn=1): ax2a.axhline(res['average'][0]+res['amplitude'][0], color='g', lw=1.0, ls='--') # average + ampl # plot chi^2 ax2b.plot(djf/MICRO, chif, 'b--', lw=1.0) ax2b.plot(dj /MICRO, chi2, 'bo' , ms=2.0) ax2b.plot(djf/MICRO, chi2, 'b--', lw=1.0) # plot 4pt method ax2a.errorbar(dj[idx_pt]/MICRO, counts[idx_pt], fmt='ms', ms=3) # no error bars ax2a.plot(djf4/MICRO, efit4, 'm--', lw=0.5) ax2a.plot(djf4/MICRO, efit4, 'm--', lw=1, label='4pt_fit') ax2a.axvline(dj0_pt/MICRO , color='m' , lw=0.5 ) # phase fit 4pt ax2a.set_ylabel('Counts [arb.]') ax2a.set_xlabel(r'$\Delta J$ [$\mu$ Tm]') #ax2b.set_ylabel(r'$\chi^2$/n.d.f.', color='b', ) ax2b.set_ylabel(r'$\chi^2$', color='b', ) ax2a.grid(True, lw=0.5, color='k', ls='--') #ax2b.grid(True, lw=0.5, color='b', ls=':') ax2a.legend(loc=0) def plot_example2(fname, axis): "test function" (avlam, delam, flux), (dj, counts, ecounts) = _read_echo_data(fname) lam0 = avlam[-1] phase0 = 2.0 phasesens = _phase_sens(lam0) # currents = np.degrees(dj*lam0*OMEGA_N)/phasesens+phase0 spectrum = Spectrum(lam=avlam, dlam=delam, flux=flux) res = fit_echo_current(currents, (counts, ecounts), spectrum, lam0=lam0, phase0=phase0, phasesens=phasesens) (phase_ef, phase_err), (curf, efit), res = res print(f"current_fit => p0 ({phase_ef:.4f} +/- {phase_err:.4f})") axis.errorbar(currents, counts, yerr=ecounts, fmt='ko') axis.plot(curf, efit, 'r-', lw=2, label='fit_echo_current') # axis.axvline(phase_ef-phase_err, color='r' , lw=0.5, ls='--' ) # axis.axvline(phase_ef , color='r' , lw=1.5, ls='-' ) # phase fit axis.axvline(phase_ef+phase_err, color='r' , lw=0.5, ls='--' ) # axis.set_xlabel('Phase Current [A]') axis.set_ylabel('Counts [arb.]') axis.grid(True, lw=0.5, color='k', ls='--') axis.legend(loc=0) if __name__ == "__main__": unittest.main(exit=False, verbosity=2) # ================================================================== # graphical tests import matplotlib.pyplot as plt #filename = os.path.join(TestDataDir, 'echo001.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo002.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo003.dat') #- 4.6-7.0 - P4 #filename = os.path.join(TestDataDir, 'echo004.dat') #- 4.0-7.0 - P2, incoherent filename = os.path.join(TestDataDir, 'echo005.dat') #- 5.6-8.0 - P4 #filename = os.path.join(TestDataDir, 'echo006.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo005.dat') #- 5.6-8.0 - P4 filename = os.path.join(TestDataDir, 'echo006.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo007.dat') #- 5.0-8.0 - P2 #filename = os.path.join(TestDataDir, 'echo008.dat') #- 7.0-8.0 - P2, partial spectrum unittest.main(exit=False, verbosity=1) signum = 0.0 if len(sys.argv)>1: filename=sys.argv[1] # graphical tests import matplotlib.pyplot as plt if len(sys.argv)>2: signum = np.sign(float(sys.argv[2])) else: signum = 0.0 _, axs = plt.subplots(2, 1, figsize=(8,8)) plot_example(filename, axs, dj0=0, sgn=signum) if not os.path.exists(filename): sys.exit("File not found: %s" % filename) _, axs = plt.subplots(3, 1, figsize=(6,8)) plot_example1(filename, axs[:-1], dj0=0, sgn=signum) plot_example2(filename, axs[-1]) plt.tight_layout() plt.show()
