new phase table fit (70ca4f69) · Commits · Zolnierczuk, Piotr / pysen

pysen/echo/init.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -4,5 +4,5 @@ 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_classic # NOQA
		from .phase_table import make_phase_table # NOQA
		from .simulator import SimpleSimulator # NOQA

pysen/echo/phase_table.py

+120 −33

Original line number	Diff line number	Diff line
		@@ -5,7 +5,9 @@ Create and manipulate NSE phase tables
		import logging
		import numpy as np

		from numpy.polynomial import polynomial as Poly
		from scipy.interpolate import RectBivariateSpline
		from scipy.optimize import curve_fit

		from ..mathutil import simple_cluster

		@@ -42,21 +44,20 @@ def get_phifit(y, cur, acur, phis, deg=2):
		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)
		p = [ Poly.polyval(cur, acur[_p]) for _p in phis]
		a = Poly.polyfit(phis, p, deg=deg)
		return Poly.polyval(y, a)

		def check_accuracy(inptab, phasetab, xi00, yphi, func):
		"check accuracy"
		def check_accuracy(inptab, phasetab, xi00, yphi, func, *args):
		"""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]
		for maincur, theta0, ph_fit, _ in np.reshape(inptab, (-1,4)):
		if args:
		ph_tab = func((maincur,theta0), *args)
		else:
		ph_tab = func(maincur, theta0)
		ph_int = interp_bf(phasetab, yphi, xi00, theta0, maincur)
		ph_dif = ph_int-ph_fit
		@@ -68,18 +69,32 @@ def check_accuracy(inptab, phasetab, xi00, yphi, func):
		log.info("max(diff)=(%.4fA,%.3f%%)", max_abs, max_rel)
		return max_abs, max_rel

		def make_phase_table_classic(inptab, polyfit=0, extent=None, **kwargs):
		"""Create phase table

		def sort_dataset(i00, phi, phase):
		"sort phase in i00 and phi"
		idx_i00 = np.argsort(i00)
		idx_phi = np.argsort(phi)
		i00 = i00[idx_i00]
		phi = phi[idx_phi]
		phase = phase[idx_phi, :]
		phase = phase[:, idx_i00]
		return i00, phi, phase

		def make_phase_table_classic(inptab, polyfit=0, extent=None, **kwargs):
		"""Create phase table (old style)
		inptab - input table with three columns i00 [Amp], phi [deg], phase current [Amp]
		polyfit -
		ncur, nphi -
		pos -
		threshold -
		"""
		threshold = kwargs.pop('threshold', 0.2)
		threshold = kwargs.pop('threshold', 0.25)
		ncurrents = kwargs.pop('ncurrents', 10)
		nphiangles = kwargs.pop('nphiangles', 11)
		#degcur = kwargs.pop('polycurrent', 2) # polynomial coefficient in i00
		#degphi = polyfit or 2 # polynomial coefficient in phi


		# 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
		@@ -99,12 +114,9 @@ def make_phase_table_classic(inptab, polyfit=0, extent=None, **kwargs):
		xi00 = np.linspace(mini00, maxi00, ncurrents)
		yphi = np.linspace(minphi, maxphi, nphiangles)

		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
		res = { 'pha': phase, 'phi': phi_set , 'i00': i00_set , 'accuracy': (0,0)}
		#
		if polyfit>0: # fit/extrapolate
		x, y = np.meshgrid(xi00, yphi)
		z = np.empty_like(x)

		@@ -112,25 +124,100 @@ def make_phase_table_classic(inptab, polyfit=0, extent=None, **kwargs):
		# 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)
		afit_i00[p] = Poly.polyfit(i00[j,:],phase[j,:], deg=2)

		# fit/extrapolate phase to
		for i, c in enumerate(xi00):
		z[:,i] = get_phifit(yphi, c, afit_i00, phi_set, deg=polyfit)

		#
		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)
		else: # interpolate using measured data
		i00_set, phi_set, phase = sort_dataset(i00_set, phi_set, phase)
		f = interp_bivspline(i00_set, phi_set, phase)
		x, y = np.meshgrid(i00_set, phi_set)
		z = f(i00_set, phi_set)

		res['cont'] = (x,y,z)
		phase_table = f(xi00,yphi)

		res['accuracy'] = 0, 0
		res['cont'] = (x,y,z)
		if polyfit>0:
		res['accuracy'] = check_accuracy(inptab, phase_table, xi00, yphi, f)
		res['xi00'] = xi00
		res['yphi'] = yphi
		res['classic'] = True

		return phase_table, res


		def make_phase_table_modern(inptab, polyfit=0, extent=None, **kwargs):
		"""Create phase table (new style)
		inptab - input table with three columns i00 [Amp], phi [deg], phase current [Amp]
		polyfit -
		ncur, nphi -
		pos -
		"""
		threshold = kwargs.pop('threshold', 0.25)
		ncurrents = kwargs.pop('ncurrents', 10)
		nphiangles = kwargs.pop('nphiangles', 11)

		degcur = kwargs.pop('polycurrent', 2) # polynomial coefficient in i00
		degphi = polyfit or 2 # polynimial coefficient in phi
		degshape = degcur+1, degphi+1

		#
		i00 = inptab[:,0]
		phi = inptab[:,1]
		phase = inptab[:,2]

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

		#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, ncurrents)
		yphi = np.linspace(minphi, maxphi, nphiangles)

		res = { 'pha': phase, 'phi': phi , 'i00': i00_set , 'accuracy': (0,0)}

		def func(x, *c):
		coeff = np.asarray(c)
		return Poly.polyval2d(x[0], x[1], coeff.reshape(degshape))

		#pylint: disable=unbalanced-tuple-unpacking
		coeff, _ = curve_fit(func, np.vstack((i00,phi)), phase, p0=np.ones(degshape).flatten())
		#pylint: enable=unbalanced-tuple-unpacking
		x, y = np.meshgrid(xi00, yphi)
		phase_table = func( (x,y), *coeff)
		res['exp'] = {}
		for c in i00_set:
		v = func((cnp.ones_like(yphi),yphi), coeff)
		res['exp'][c] = (yphi, v)
		res['cont'] = (x,y, phase_table)
		res['accuracy'] = check_accuracy(inptab, phase_table, xi00, yphi, func, coeff)
		res['xi00'] = xi00
		res['yphi'] = yphi
		#res['inptab'] = inptab.copy()
		res['classic'] = False
		return phase_table, res


		def make_phase_table(inptab, polyfit=0, extent=None, **kwargs):
		"""Create phase table
		inptab - input table with three columns i00 [Amp], phi [deg], phase current [Amp]
		polyfit -
		ncur, nphi -
		pos -
		"""
		if kwargs.pop('classic', False):
		return make_phase_table_classic(inptab, polyfit=polyfit, extent=extent, **kwargs)
		return make_phase_table_modern(inptab, polyfit=polyfit, extent=extent, **kwargs)
		#EOF

pysen/plot/miscplotlib.py

+38 −27

Original line number	Diff line number	Diff line
		@@ -13,37 +13,47 @@ import h5py

		from ..config import phi_limits
		from ..inout import convert_files
		from ..echo import get_symmetry_phase, make_phase_table_classic
		from ..echo import get_symmetry_phase, make_phase_table

		def plot_phase_table(res):
		def plot_phase_table(fig, res):
		""" plot phase table"""
		plt.subplot(1,2,1)
		#
		i00 = res['i00']
		phi = res['phi']
		pha = res['pha']
		x, y, z = res['cont']
		classic = res['classic']

		ax = fig.subplots(1,2)

		lw, ms= 1.5, 3

		for i, c in enumerate(res['i00']):
		if classic:
		fig.suptitle("Phase Table Plot (classic)")
		for i, c in enumerate(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}')
		ax[0].plot(phi,pha[:,i],'s--', lw=lw, ms=ms, 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}')
		ax[0].plot(phi,pha[:,i],'ks', ms=ms)
		ax[0].plot(yphi, vals, '--', lw=lw, label=rf'i$_{{00}}$={c:.1f}')
		else:
		fig.suptitle("Phase Table Plot ")
		ax[0].plot(phi, pha , 'ks', ms=ms, label='data')
		for i, c in enumerate(i00):
		yphi, vals = res['exp'][c]
		ax[0].plot(yphi, vals, '--', lw=lw, label=rf'i$_{{00}}$={c:.1f}')

		plt.xlabel(r'$\phi$ [deg]')
		plt.ylabel(r'i$_5$ [A]')
		plt.legend(loc='best')
		plt.grid()
		ax[0].set_xlabel(r'$\phi$ [deg]')
		ax[0].set_ylabel(r'i$_5$ [A]')
		ax[0].legend(loc='best')
		ax[0].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')
		x, y, z = res['cont']
		im = ax[1].contourf(y,x,z,21)
		fig.colorbar(im, ax=ax[1])
		ax[1].set_xlabel(r'$\phi$ [deg]')
		ax[1].set_ylabel(r'i$_{00}$ [A]')

		def print_phase_table(phase, res, pos='p2', filename=None):
		"print phasetable"
		@@ -52,7 +62,7 @@ def print_phase_table(phase, res, pos='p2', filename=None):
		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(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 ")
		@@ -81,6 +91,7 @@ def action_phase_table(filenames, **kwargs):
		polyfit = kwargs.pop('polyfit')
		pos = kwargs.pop('pos')
		threshold = kwargs.pop('threshold')
		classic = kwargs.pop('classic')
		#
		log = logging.getLogger()
		inptab = None
		@@ -100,11 +111,11 @@ def action_phase_table(filenames, **kwargs):
		# FIXME hardcoded constants
		extent = (0.0, 90.0, -2.0, phi_limits(pos)[1]) if polyfit>0 else None

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

		phase, res = make_phase_table_classic(inptab, polyfit, extent=extent, threshold=threshold)
		phase, res = make_phase_table(inptab, polyfit, extent=extent, threshold=threshold, classic=classic)
		print_phase_table(phase, res, pos=pos, filename=kwargs.get('savefile'))
		plot_phase_table(res)
		plot_phase_table(fig, res)
		if savefig:
		plt.savefig(savefig)
		fig.savefig(savefig)
		#EOF

pysen/ui/nseplot.py

+10 −6

Original line number	Diff line number	Diff line
		@@ -235,10 +235,12 @@ def add_det_options(subparser, parents=None):
		help='detimage plot (nxs.h5)',
		description='show detector image')
		pars.add_argument('file', metavar='filename', help='file to process', nargs='+')
		pars.set_defaults(minchan=0, maxchan=TDC_MAXCH,
		npix=TDC_MAXPIX, log=False, edge=False, threshold=None)
		pars.add_argument('--num-pix', '-N', dest='npix', metavar='N', type=int, help='set number of pixels, (default %(default)s)')
		pars.add_argument('--edge' , '-E', dest='edge', action='store_true', help='edge detection')
		pars.set_defaults(minchan=0, maxchan=TDC_MAXCH, npix=TDC_MAXPIX,
		log=False, edge=False, threshold=None)
		pars.add_argument('--num-pix' , '-N', dest='npix', metavar='N', type=int,
		help='set number of pixels, (default %(default)s)')
		pars.add_argument('--edge' , '-E', dest='edge', action='store_true',
		help='edge detection')
		pars.add_argument('--threshold', '-T', dest='threshold', metavar='THR', type=float,
		help='set threshold (need more explanation here)')

		@@ -333,12 +335,14 @@ def add_ptab_options(subparser, parents=None):
		pars = subparser.add_parser('phase_table', parents=parents,
		aliases=['pt',],
		help='phase table', description='phase table')
		pars.set_defaults(pos='p2', polyfit=2, threshold=0.3)
		pars.set_defaults(pos='p2', polyfit=2, threshold=0.3, classic=False)
		pars.add_argument('file', metavar='filename', help='file to process', nargs='+')
		pars.add_argument('--pos' , '-p', dest='pos', choices=INST_POSITIONS,
		help='set instrument position: (default %(default)s)')
		pars.add_argument('--poly', '-P', dest='polyfit', metavar='deg', type=int,
		help='use polynomial fit of degree (default %(default)s)')
		pars.add_argument('--classic', dest='classic', action='store_true',
		help='use classic phase table algorithm')
		pars.add_argument('--threshold', dest='threshold', metavar='thres', type=float,
		help='clustering threshold (advanced option)')

test/test_echo_phase_table.py

+44 −7

Original line number	Diff line number	Diff line
		@@ -9,6 +9,7 @@ import numpy as np

		import pysen.config
		import pysen.echo.phase_table as ept
		from pysen.mathutil import simple_cluster

		TestDataDir = os.path.join(os.path.dirname(__file__), 'data')

		@@ -16,6 +17,10 @@ class PhaseTableTestCase(unittest.TestCase):
		"Phase Table Test Cases"
		def setUp(self):
		"setup method"
		self.pt = np.loadtxt(os.path.join(TestDataDir, 'phase_table_inp.dat'))
		self.extent = 0.0, 90.0, -2.0, pysen.config.phi_limits('p2')[1]
		self.threshold = 0.2
		self.i00_exp = simple_cluster(self.pt[:,0], threshold=0.2)

		def tearDown(self):
		"Tear-down method"
		@@ -27,7 +32,7 @@ class PhaseTableTestCase(unittest.TestCase):
		[0.9494, 1.5486, 2.1330, 2.7025, 3.2572, 3.7969, 4.3219],
		[1.0163, 1.6466, 2.2619, 2.8622, 3.4474, 4.0177, 4.5729],
		[0.9869, 1.6641, 2.3250, 2.9695, 3.5978, 4.2097, 4.8054],])
		i00_exp = np.asarray([0.2154, 19.2372, 38.5685, 57.9000, 77.2317])
		#i00_exp = np.asarray([0.2154, 19.2372, 38.5685, 57.9000, 77.2317])
		phi_exp = np.asarray([3.64, 7.30, 10.96, 14.63, 18.32, 29.51])
		pha_exp = np.asarray([ [0.6692, 1.4257, 2.1292, 2.8161, 3.4959],
		[0.7590, 1.5177, 2.2197, 2.9048, 3.5858],
		@@ -36,18 +41,50 @@ class PhaseTableTestCase(unittest.TestCase):
		[0.9287, 1.7014, 2.4309, 3.1365, 3.8401],
		[1.0169, 1.8248, 2.5891, 3.3329, 4.0781]])
		#
		pt = np.loadtxt(os.path.join(TestDataDir, 'phase_table_inp.dat'))
		extent = 0.0, 90.0, -2.0, pysen.config.phi_limits('p2')[1]
		phase_act, pres = ept.make_phase_table_classic(pt, polyfit=2,
		ncurrents=7, nphiangles=5, extent=extent)
		i00_exp = self.i00_exp
		phase_act, pres = ept.make_phase_table_classic(self.pt, polyfit=2,
		ncurrents=7, nphiangles=5, extent=self.extent)
		#
		self.assertEqual(phase_act.shape, phase_exp.shape, "phase.shape")
		self.assertTrue(np.allclose(phase_act, phase_exp, atol=0.001), "phase_table")
		#
		self.assertEqual(pres['phi'][:,0].shape, phi_exp.shape, "phi")
		self.assertEqual(pres['phi'].shape, phi_exp.shape, "phi")
		self.assertEqual(pres['i00'].shape, i00_exp.shape, "i00")
		self.assertEqual(pres['pha'].shape, pha_exp.shape, "pha")
		self.assertTrue(np.allclose(pres['phi'][:,0], phi_exp, atol=0.01), "phi")
		self.assertTrue(np.allclose(pres['phi'], phi_exp, atol=0.01), "phi")
		self.assertTrue(np.allclose(pres['i00'], i00_exp, atol=0.001), "i00")
		self.assertTrue(np.allclose(pres['pha'], pha_exp, atol=0.001), "pha")
		self.assertAlmostEqual(pres['accuracy'][0], 0.027, 3, "max_abs")
		self.assertAlmostEqual(pres['accuracy'][1], 1.782, 3, "max_rel")

		def test_phase_table_modern(self):
		"test phase table 'modern'"
		# This test is similar to the classic one, but uses the modern function
		phase_exp= np.asarray([ [0.5261, 1.1109, 1.6776, 2.2263, 2.7570, 3.2695, 3.7641],
		[0.7860, 1.3700, 1.9382, 2.4906, 3.0270, 3.5476, 4.0523],
		[0.9494, 1.5486, 2.1330, 2.7025, 3.2572, 3.7969, 4.3219],
		[1.0163, 1.6466, 2.2619, 2.8622, 3.4474, 4.0177, 4.5729],
		[0.9869, 1.6641, 2.3250, 2.9695, 3.5978, 4.2097, 4.8054],])
		pha_exp = np.asarray([ [0.6692, 1.4257, 2.1292, 2.8161, 3.4959],
		[0.7590, 1.5177, 2.2197, 2.9048, 3.5858],
		[0.8282, 1.5861, 2.3007, 2.9871, 3.6841],
		[0.8840, 1.6408, 2.3680, 3.0574, 3.7624],
		[0.9287, 1.7014, 2.4309, 3.1365, 3.8401],
		[1.0169, 1.8248, 2.5891, 3.3329, 4.0781]]).flatten()
		#
		phi_exp = self.pt[:,1]
		i00_exp = self.i00_exp

		phase_act, pres = ept.make_phase_table_modern(self.pt, polyfit=2,
		ncurrents=7, nphiangles=5, extent=self.extent)
		#
		self.assertEqual(phase_act.shape, phase_exp.shape, "phase.shape")
		self.assertTrue(np.allclose(phase_act, phase_exp, atol=0.001), "phase_table")
		#
		self.assertEqual(pres['phi'].shape, phi_exp.shape, "phi")
		self.assertEqual(pres['i00'].shape, i00_exp.shape, "i00")
		self.assertEqual(pres['pha'].shape, pha_exp.shape, "pha")
		self.assertTrue(np.allclose(pres['phi'], phi_exp, atol=0.01), "phi")
		self.assertTrue(np.allclose(pres['i00'], i00_exp, atol=0.001), "i00")
		self.assertTrue(np.allclose(pres['pha'], pha_exp, atol=0.001), "pha")
		self.assertAlmostEqual(pres['accuracy'][0], 0.027, 3, "max_abs")