before big refactoring (dfade60d) · Commits · Zolnierczuk, Piotr / pysen

.pylint.rc

+2 −2

Original line number	Diff line number	Diff line
		@@ -59,8 +59,8 @@ confidence=
		# --enable=similarities". If you want to run only the classes checker, but have
		# no Warning level messages displayed, use"--disable=all --enable=classes
		# --disable=W"
		disable=execfile-builtin,input-builtin,setslice-method,getslice-method,old-raise-syntax,useless-suppression,oct-method,map-builtin-not-iterating,parameter-unpacking,print-statement,buffer-builtin,next-method-called,intern-builtin,backtick,cmp-method,import-star-module-level,coerce-method,using-cmp-argument,raw_input-builtin,file-builtin,suppressed-message,delslice-method,dict-iter-method,reduce-builtin,old-division,raising-string,xrange-builtin,long-builtin,no-absolute-import,range-builtin-not-iterating,zip-builtin-not-iterating,hex-method,round-builtin,indexing-exception,filter-builtin-not-iterating,standarderror-builtin,nonzero-method,apply-builtin,old-ne-operator,old-octal-literal,coerce-builtin,metaclass-assignment,cmp-builtin,unpacking-in-except,dict-view-method,basestring-builtin,unicode-builtin,reload-builtin,long-suffix,unichr-builtin,bad-whitespace,invalid-name,locally-disabled,fixme
		#,invalid-name,no-name-in-module,no-member
		disable=fixme,bad-whitespace,invalid-name
		#,too-many-statements,too-many-locals
		#disable=all
		enable=similarities

pysen/config.py

+59 −51

Original line number	Diff line number	Diff line
		@@ -198,6 +198,7 @@ def get_attenuation(att_pos, wavelength, att_table=None):
		fac = (wavelength-8.0)/3.0 # 11.0 - 8.0
		return att08 + datt*fac

		# ========================for coverage ==============================================

		def _write_ql_min_max(out, qmin, lmax, dlam):
		"write out min max lambda, q"
		@@ -212,10 +213,10 @@ def _write_ql_min_max(out, qmin, lmax, dlam):
		out.write("\n")
		return dict(lmax=lmax,lmin=lmin,lave=lave,qmin=qmin,qmax=qmax,qave=qave)

		def _coverage_check_tau(tau, tlimits, lmax, limits_error=False):
		def _coverage_check_tau(taus, tlimits, lmax, limits_error=False):
		tau_min, tau_max = tlimits
		if tau_min<=tau or tau<=tau_max:
		return
		for tau in taus:
		if tau<tau_min or tau_max<tau:
		msg = ("tau set outside the bounds %s (%s,%s) for lambda=%s A"
		% (tau, tau_min, tau_max, lmax))
		if limits_error:
		@@ -223,6 +224,30 @@ def _coverage_check_tau(tau, tlimits, lmax, limits_error=False):
		else:
		warnings.warn(msg, RuntimeWarning)

		def _coverage_check_theta(qmin, lmax, detpos):
		theta = get_theta(qmin, lmax)
		mintheta, maxtheta = [ radians(_x) for _x in phi_limits(detpos) ]
		if isnan(theta) or theta < mintheta or maxtheta < theta:
		msg = ( "Scattering angle %.5g deg outside instrument limits (%s,%s) for position %s"
		% (degrees(theta), degrees(mintheta), degrees(maxtheta), detpos))
		raise RuntimeError(msg)
		return theta

		def get_theta_pix(x,y,z, sa, ca):
		xrot = xca + zsa
		yrot = y
		zrot = -xsa + zca
		r = sqrt(xrot2+yrot2+zrot**2)
		return arccos(zrot/r)

		def _gen_plot_data(out, q, tau, lmax, l1, l2):
		lave = 0.5*(l1+l2)
		out.write("\ttau set=%7.3f" % tau)
		xtau = tau(lave/lmax)*3
		dt1 = abs(tau(l1/lmax)*3-xtau)
		dt2 = abs(tau(l2/lmax)*3-xtau)
		out.write("\tact=%7.3f\n" % xtau)
		return (xtau, q, dt1, dt2)

		def coverage(lmax, qmin, **kwargs):
		"calculate q-tau coverage"
		@@ -239,39 +264,27 @@ def coverage(lmax, qmin, **kwargs):

		out = StringIO()


		if taus is None:
		taus = 5
		taus = taus or 5
		if tbins is None:
		tbins = 5
		elif len(tbins)==1:
		tbins = tbins[0]

		tau_min, tau_max = tau_limits(lmax, mode)
		if np.isscalar(taus):
		taus = logspace(log10(tau_min), log10(tau_max), taus)
		tau_min, tau_max = taus[0], taus[-1] # to avoid rounding errors
		for tau in taus:
		if tau<tau_min or tau_max<tau:
		msg = ( "Fourier time %s ns outside instrument limits [%.3f,%.3f] for wavelength %s A"
		% ( tau, tau_min, tau_max, lmax))
		raise RuntimeError(msg)
		#
		_coverage_check_tau(taus, (tau_min, tau_max), lmax, limits_error=True)

		if np.isscalar(tbins):
		tbins = linspace(0, ntbin, tbins+1)
		tbin1 = tbins[:-1]
		tbin2 = tbins[1:]

		theta = get_theta(qmin, lmax)
		mintheta, maxtheta = [ radians(_x) for _x in phi_limits(detpos) ]
		if isnan(theta) or theta < mintheta or maxtheta < theta:
		msg = ( "Scattering angle %.5g deg outside instrument limits (%s,%s) for position %s"
		% (degrees(theta), degrees(mintheta), degrees(maxtheta), detpos))
		raise RuntimeError(msg)
		theta0 = _coverage_check_theta(qmin, lmax, detpos)

		lq = lmax*qmin
		out.write("phi=%4.1f tau_range=(%.3f,%.3f)\n" % (degrees(theta), tau_min, tau_max))
		results['theta'] = theta
		out.write("phi=%4.1f tau_range=(%.3f,%.3f)\n" % (degrees(theta0), tau_min, tau_max))
		results['theta'] = theta0
		results['taus'] = taus
		results['tbins'] = tbins
		results['tau_min'] = tau_min
		@@ -283,42 +296,37 @@ def coverage(lmax, qmin, **kwargs):
		l1 = l_binning(lmax, it1, dlam=dlam, nt=ntbin)
		l2 = l_binning(lmax, it2, dlam=dlam, nt=ntbin)
		lave = 0.5*(l1+l2)

		q, dq = q_binning(lmax, qmin, it1, it2, dlam=dlam, nt=ntbin)
		out.write("q=(%.3f +/-%.3f)\t" % (q, dq))
		out.write("lave=%.3f [%d %d]\n" % (lave, it1, it2))
		_coverage_check_tau(taus, (tau_min, tau_max), lmax, limits_error=limits)
		for tau in taus:
		_coverage_check_tau(tau, (tau_min, tau_max), lmax, limits_error=limits)
		out.write("\ttau set=%7.3f" % tau)
		xtau = tau(lave/lmax)*3
		dt1 = abs(tau(l1/lmax)*3-xtau)
		dt2 = abs(tau(l2/lmax)*3-xtau)
		out.write("\tact=%7.3f\n" % xtau)
		plot_data.append((xtau, q, dt1, dt2, dq, dq, tau))
		res = _gen_plot_data(out, q, tau, lmax, l1,l2)
		plot_data.append(res+(dq,dq,tau))
		out.write("\n")
		results['plot_data'] = np.asarray(plot_data).T

		# limits
		lq = lmax*qmin
		q = linspace(qmin, lq/(lmax-dlam), len(tbins)+1)
		t1, t2 = tau_limits(lq/q, mode)
		limits = [(q, t1, t2)]

		# add coverage for left and right
		ca = cos(theta)
		sa = sin(theta)
		z = L2
		ca = cos(theta0)
		sa = sin(theta0)
		# TODO this is way overkill
		xlims = linspace(-XDET/2.0, XDET/2.0, NXCHAN//2)
		ylims = linspace(-YDET/2.0, YDET/2.0, NYCHAN//2)
		for x in xlims:
		for y in ylims:
		xrot = xca + zsa
		yrot = y
		zrot = -xsa + zca
		r = sqrt(xrot2+yrot2+zrot**2)
		_theta = arccos(zrot/r)
		lq = get_q(lmax, _theta)*lmax
		q = linspace(get_q(lmax, _theta), get_q(lmax-dlam, _theta), len(tbins)+1)
		xlims = linspace(-XDET/2.0, XDET/2.0, NXCHAN)
		ylims = linspace( YDET/2.0,-YDET/2.0, NYCHAN)

		for i,x in enumerate(xlims):
		for j,y in enumerate(ylims):
		theta1 = get_theta_pix(x,y,L2,sa,ca)
		theta2, _ = pixel_angle(i,j,theta0)
		print(theta1, theta2)
		theta = theta2
		lq = get_q(lmax, theta)*lmax
		q = linspace(get_q(lmax, theta), get_q(lmax-dlam, theta), len(tbins)+1)
		t1, t2 = tau_limits(lq/q, mode)
		limits.append((q, t1, t2))

pysen/echo/plot.py

+0 −1

Original line number	Diff line number	Diff line
		@@ -80,4 +80,3 @@ def display_echo(hdfile, iecho=None, show=False):
		plt.grid(True)
		if show:
		plt.show()
		return

pysen/echo/reduce.py

+80 −0

Original line number	Diff line number	Diff line
		@@ -14,6 +14,86 @@ from .utils import ( create_mask, calc_sqt, findlevels)

		DEFAULT_NRINGS = 5


		def reduce_data(hdfile, **kwargs):
		"echo plot"
		#
		center_only = kwargs.pop('center_only', False)
		#
		npix = kwargs.pop('npix', 2) # pix
		#
		tbin1 = kwargs.pop('tbin1',0) # TOF bins
		tbin2 = kwargs.pop('tbin2',None)
		#
		min_counts = kwargs.pop('min_counts', 10.0)
		max_chi2 = kwargs.pop('max_chi2' ,500.0)
		min_amp = kwargs.pop('min_amp' , 0.0)
		# FIXME: HARDCODED
		xpix1 = kwargs.pop('xpix1', 10)
		xpix2 = kwargs.pop('xpix2', 22)
		ypix1 = kwargs.pop('ypix1', 10)
		ypix2 = kwargs.pop('ypix2', 22)

		base = os.path.splitext(os.path.basename(hdfile.filename))[0]
		comment = hdfile.attrs['master_comment'].decode("utf-8")

		ntaus = hdfile.attrs['scan_length']
		n_idx = { 'dn': hdfile['/'].attrs['point_to_down'],
		'up': hdfile['/'].attrs['point_to_up'] ,
		'nphases': hdfile['/'].attrs['no_of_phases'] }
		# sample
		sample_name = hdfile['/sample'].attrs['sample'][0].decode("utf-8")
		print(sample_name)
		# detector
		nt = hdfile['/detector'].attrs['no_t_channels']
		ny = hdfile['/detector'].attrs['no_y_channels']
		nx = hdfile['/detector'].attrs['no_x_channels']
		nph = n_idx['dn']

		with h5py.File(base+'_results.h5', 'w') as results:
		for attr in hdfile.attrs:
		results.attrs[attr] = hdfile.attrs[attr]
		results.attrs['reduction'] = time.ctime()


		for echo in list(hdfile['/data'].values()):
		print("processing", echo.name)
		pass

		phase = echo['phase']
		physics = echo['phys']
		params = echo['params']
		tech = echo['tech']

		tau0 = physics.attrs['fouriertime']/NANOSECOND
		q0 = physics.attrs['hkl'][0]
		lam0 = physics.attrs['lambda']
		phase0 = float(tech.attrs['i5'][0])

		lmax = np.array(params['lambdaTable']).flatten()
		dlam = np.ones_like(lmax)*(lmax[1]-lmax[0])

		phasesens = float(params['phaseangle'].attrs['sensitivities'][0])

		det = phase['detector'][...]
		pcha = phase['proton_charge'][...]
		cur = phase['phase_current'][:, 0] # actual value
		pcha0 = np.average(pcha)
		pcha = np.tile(pcha, (ny, nx, nt)).reshape(n_idx['nphases'],ny,nx,nt)
		pcha = pcha/float(pcha0)
		det = det/pcha
		wlen = np.sum(det, axis=(0,1,2))

		spectrum = Spectrum(lam=lmax, dlam=dlam, flux=wlen)

		y = np.sum(det[:,:,:,tbin1:tbin2], axis=(1,2,3))
		ey = np.sqrt(y)
		dn = y[nph:n_idx['up']]
		up = y[n_idx['up']:]
		phase_ef0, (xef,yef), res = fit_echo_current(cur[:nph], (y[:nph], ey[:nph]), spectrum,
		lam0=lam0, phase0=phase0, phasesens=phasesens)


		def process_data(hdfile, dj0=None, qbins=None, tbins=None):
		"""
		process echo data

pysen/io/hdf.py

+3 −3

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


		from ..revision import version
		from ..revision import version as version_info
		from ..iostrings import encode
		from .reader import read_echo, read_magnetic

		@@ -140,7 +140,7 @@ def echo_to_hdf(filename, outdir, **kwargs):

		outfile, _ = os.path.splitext(os.path.basename(filename))
		outfile = os.path.join(outdir, outfile + ".h5")
		converter = HdfConverter(outfile, mode=mode, version=version())
		converter = HdfConverter(outfile, mode=mode, version=version_info())
		data = read_echo(filename)
		converter.write(data, compression=compression, schema=schema)
		return outfile
		@@ -154,7 +154,7 @@ def magnetic_to_hdf(filename, outdir, **kwargs):

		outfile, _ = os.path.splitext(os.path.basename(filename))
		outfile = os.path.join(outdir, outfile + ".h5")
		converter = HdfConverter(outfile, mode=mode, version=version())
		converter = HdfConverter(outfile, mode=mode, version=version_info())
		data = read_magnetic(filename)
		converter.write(data, compression=compression, schema=schema)
		return outfile