pr.py @ 2f85af7

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 2f85af7 was 7432acb, checked in by andyfaff, 8 years ago
MAINT: search+replace '!= None' by 'is not None'
Property mode set to `100644`
File size: 45.8 KB

Rev	Line
[959eb01]	1	"""
	2	P(r) perspective for SasView
	3	"""
	4	################################################################################
	5	#This software was developed by the University of Tennessee as part of the
	6	#Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
	7	#project funded by the US National Science Foundation.
	8	#
	9	#See the license text in license.txt
	10	#
	11	#copyright 2009, University of Tennessee
	12	################################################################################
	13
	14
	15	# Make sure the option of saving each curve is available
	16	# Use the I(q) curve as input and compare the output to P(r)
	17
	18	import sys
	19	import wx
	20	import logging
	21	import time
	22	import math
	23	import numpy as np
	24	import pylab
	25	from sas.sasgui.guiframe.gui_manager import MDIFrame
	26	from sas.sasgui.guiframe.dataFitting import Data1D
	27	from sas.sasgui.guiframe.events import NewPlotEvent
	28	from sas.sasgui.guiframe.events import StatusEvent
	29	from sas.sasgui.guiframe.gui_style import GUIFRAME_ID
	30	from sas.sascalc.pr.invertor import Invertor
	31	from sas.sascalc.dataloader.loader import Loader
	32	import sas.sascalc.dataloader
	33
	34	from pr_widgets import load_error
	35	from sas.sasgui.guiframe.plugin_base import PluginBase
	36
	37	logger = logging.getLogger(__name__)
	38
	39
	40	PR_FIT_LABEL = r"$P_{fit}(r)$"
	41	PR_LOADED_LABEL = r"$P_{loaded}(r)$"
	42	IQ_DATA_LABEL = r"$I_{obs}(q)$"
	43	IQ_FIT_LABEL = r"$I_{fit}(q)$"
	44	IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
	45	GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
	46	GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
	47
	48
	49
	50	class Plugin(PluginBase):
	51	"""
	52	P(r) inversion perspective
	53	"""
	54	DEFAULT_ALPHA = 0.0001
	55	DEFAULT_NFUNC = 10
	56	DEFAULT_DMAX = 140.0
	57
	58	def __init__(self):
	59	PluginBase.__init__(self, name="Pr Inversion")
	60	## Simulation window manager
	61	self.simview = None
	62
	63	## State data
	64	self.alpha = self.DEFAULT_ALPHA
	65	self.nfunc = self.DEFAULT_NFUNC
	66	self.max_length = self.DEFAULT_DMAX
	67	self.q_min = None
	68	self.q_max = None
	69	self.has_bck = False
	70	self.slit_height = 0
	71	self.slit_width = 0
	72	## Remember last plottable processed
	73	self.last_data = ""
	74	self._current_file_data = None
	75	## Time elapsed for last computation [sec]
	76	# Start with a good default
	77	self.elapsed = 0.022
	78	self.iq_data_shown = False
	79
	80	## Current invertor
	81	self.invertor = None
	82	self.pr = None
	83	self.data_id = IQ_DATA_LABEL
	84	# Copy of the last result in case we need to display it.
	85	self._last_pr = None
	86	self._last_out = None
	87	self._last_cov = None
	88	## Calculation thread
	89	self.calc_thread = None
	90	## Estimation thread
	91	self.estimation_thread = None
	92	## Result panel
	93	self.control_panel = None
	94	## Currently views plottable
	95	self.current_plottable = None
	96	## Number of P(r) points to display on the output plot
	97	self._pr_npts = 51
	98	self._normalize_output = False
	99	self._scale_output_unity = False
	100
	101	## List of added P(r) plots
	102	self._added_plots = {}
	103	self._default_Iq = {}
	104	self.list_plot_id = []
	105
	106	# Associate the inversion state reader with .prv files
	107	from inversion_state import Reader
	108
	109	# Create a CanSAS/Pr reader
	110	self.state_reader = Reader(self.set_state)
	111	self._extensions = '.prv'
	112	l = Loader()
	113	l.associate_file_reader('.prv', self.state_reader)
	114	#l.associate_file_reader(".svs", self.state_reader)
	115
	116	# Log startup
	117	logger.info("Pr(r) plug-in started")
	118
	119	def delete_data(self, data_id):
	120	"""
	121	delete the data association with prview
	122	"""
	123	self.control_panel.clear_panel()
	124
	125	def get_data(self):
	126	"""
	127	Returns the current data
	128	"""
	129	return self.current_plottable
	130
	131	def set_state(self, state=None, datainfo=None):
	132	"""
	133	Call-back method for the inversion state reader.
	134	This method is called when a .prv file is loaded.
	135
	136	:param state: InversionState object
	137	:param datainfo: Data1D object [optional]
	138
	139	"""
	140	try:
	141	if datainfo.__class__.__name__ == 'list':
	142	if len(datainfo) >= 1:
	143	data = datainfo[0]
	144	else:
	145	data = None
	146	else:
	147	data = datainfo
	148	if data is None:
	149	msg = "Pr.set_state: datainfo parameter cannot "
	150	msg += "be None in standalone mode"
	151	raise RuntimeError, msg
	152
	153	# Ensuring that plots are coordinated correctly
	154	t = time.localtime(data.meta_data['prstate'].timestamp)
	155	time_str = time.strftime("%b %d %H:%M", t)
	156
	157	# Check that no time stamp is already appended
	158	max_char = data.meta_data['prstate'].file.find("[")
	159	if max_char < 0:
	160	max_char = len(data.meta_data['prstate'].file)
	161
	162	datainfo.meta_data['prstate'].file = \
	163	data.meta_data['prstate'].file[0:max_char]\
	164	+ ' [' + time_str + ']'
	165
	166	data.filename = data.meta_data['prstate'].file
	167	# TODO:
	168	#remove this call when state save all information about the gui data
	169	# such as ID , Group_ID, etc...
	170	#make self.current_plottable = datainfo directly
	171	self.current_plottable = self.parent.create_gui_data(data, None)
	172	self.current_plottable.group_id = data.meta_data['prstate'].file
	173
	174	# Make sure the user sees the P(r) panel after loading
	175	#self.parent.set_perspective(self.perspective)
	176	self.on_perspective(event=None)
	177	# Load the P(r) results
	178	#state = self.state_reader.get_state()
	179	data_dict = {self.current_plottable.id:self.current_plottable}
	180	self.parent.add_data(data_list=data_dict)
	181	wx.PostEvent(self.parent, NewPlotEvent(plot=self.current_plottable,
	182	title=self.current_plottable.title))
	183	self.control_panel.set_state(state)
	184	except:
	185	logger.error("prview.set_state: %s" % sys.exc_value)
	186
	187
	188	def help(self, evt):
	189	"""
	190	Show a general help dialog.
	191
	192	:TODO: replace the text with a nice image
	193
	194	"""
	195	from inversion_panel import HelpDialog
	196	dialog = HelpDialog(None, -1)
	197	if dialog.ShowModal() == wx.ID_OK:
	198	dialog.Destroy()
	199	else:
	200	dialog.Destroy()
	201
	202	def _fit_pr(self, evt):
	203	"""
	204	"""
	205	# Generate P(r) for sphere
	206	radius = 60.0
	207	d_max = 2 * radius
	208
	209	r = pylab.arange(0.01, d_max, d_max / 51.0)
	210	M = len(r)
	211	y = np.zeros(M)
	212	pr_err = np.zeros(M)
	213
	214	total = 0.0
	215	for j in range(M):
	216	value = self.pr_theory(r[j], radius)
	217	total += value
	218	y[j] = value
	219	pr_err[j] = math.sqrt(y[j])
	220
	221	y = y / total * d_max / len(r)
	222
	223	# Perform fit
	224	pr = Invertor()
	225	pr.d_max = d_max
	226	pr.alpha = 0
	227	pr.x = r
	228	pr.y = y
	229	pr.err = pr_err
	230	out, cov = pr.pr_fit()
	231	for i in range(len(out)):
	232	print "%g +- %g" % (out[i], math.sqrt(cov[i][i]))
	233
	234	# Show input P(r)
	235	title = "Pr"
	236	new_plot = Data1D(pr.x, pr.y, dy=pr.err)
	237	new_plot.name = "P_{obs}(r)"
	238	new_plot.xaxis("\\rm{r}", 'A')
	239	new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
	240	new_plot.group_id = "P_{obs}(r)"
	241	new_plot.id = "P_{obs}(r)"
	242	new_plot.title = title
	243	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	244	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
	245
	246	# Show P(r) fit
	247	self.show_pr(out, pr)
	248
	249	# Show I(q) fit
	250	q = pylab.arange(0.001, 0.1, 0.01 / 51.0)
	251	self.show_iq(out, pr, q)
	252
	253	def show_shpere(self, x, radius=70.0, x_range=70.0):
	254	"""
	255	"""
	256	# Show P(r)
	257	y_true = np.zeros(len(x))
	258
	259	sum_true = 0.0
	260	for i in range(len(x)):
	261	y_true[i] = self.pr_theory(x[i], radius)
	262	sum_true += y_true[i]
	263
	264	y_true = y_true / sum_true * x_range / len(x)
	265
	266	# Show the theory P(r)
	267	new_plot = Data1D(x, y_true)
	268	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	269	new_plot.name = "P_{true}(r)"
	270	new_plot.xaxis("\\rm{r}", 'A')
	271	new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
	272	new_plot.id = "P_{true}(r)"
	273	new_plot.group_id = "P_{true}(r)"
	274	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	275	#Put this call in plottables/guitools
	276	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot,
	277	title="Sphere P(r)"))
	278
	279	def get_npts(self):
	280	"""
	281	Returns the number of points in the I(q) data
	282	"""
	283	try:
	284	return len(self.pr.x)
	285	except:
	286	return 0
	287
	288	def show_iq(self, out, pr, q=None):
	289	"""
	290	Display computed I(q)
	291	"""
	292	qtemp = pr.x
[ac07a3a]	293	if q is not None:
[959eb01]	294	qtemp = q
	295
	296	# Make a plot
	297	maxq = -1
	298	for q_i in qtemp:
	299	if q_i > maxq:
	300	maxq = q_i
	301
	302	minq = 0.001
	303
	304	# Check for user min/max
[ac07a3a]	305	if pr.q_min is not None:
[959eb01]	306	minq = pr.q_min
[ac07a3a]	307	if pr.q_max is not None:
[959eb01]	308	maxq = pr.q_max
	309
	310	x = pylab.arange(minq, maxq, maxq / 301.0)
	311	y = np.zeros(len(x))
	312	err = np.zeros(len(x))
	313	for i in range(len(x)):
	314	value = pr.iq(out, x[i])
	315	y[i] = value
	316	try:
	317	err[i] = math.sqrt(math.fabs(value))
	318	except:
	319	err[i] = 1.0
	320	print "Error getting error", value, x[i]
	321
	322	new_plot = Data1D(x, y)
	323	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	324	new_plot.name = IQ_FIT_LABEL
	325	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
	326	new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
	327	title = "I(q)"
	328	new_plot.title = title
	329
	330	# If we have a group ID, use it
	331	if 'plot_group_id' in pr.info:
	332	new_plot.group_id = pr.info["plot_group_id"]
	333	new_plot.id = IQ_FIT_LABEL
	334	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	335
	336	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
	337
	338	# If we have used slit smearing, plot the smeared I(q) too
	339	if pr.slit_width > 0 or pr.slit_height > 0:
	340	x = pylab.arange(minq, maxq, maxq / 301.0)
	341	y = np.zeros(len(x))
	342	err = np.zeros(len(x))
	343	for i in range(len(x)):
	344	value = pr.iq_smeared(out, x[i])
	345	y[i] = value
	346	try:
	347	err[i] = math.sqrt(math.fabs(value))
	348	except:
	349	err[i] = 1.0
	350	print("Error getting error", value, x[i])
	351
	352	new_plot = Data1D(x, y)
	353	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	354	new_plot.name = IQ_SMEARED_LABEL
	355	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
	356	new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
	357	# If we have a group ID, use it
	358	if 'plot_group_id' in pr.info:
	359	new_plot.group_id = pr.info["plot_group_id"]
	360	new_plot.id = IQ_SMEARED_LABEL
	361	new_plot.title = title
	362	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	363	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title=title))
	364
	365	def _on_pr_npts(self, evt):
	366	"""
	367	Redisplay P(r) with a different number of points
	368	"""
	369	from inversion_panel import PrDistDialog
	370	dialog = PrDistDialog(None, -1)
	371	dialog.set_content(self._pr_npts)
	372	if dialog.ShowModal() == wx.ID_OK:
	373	self._pr_npts = dialog.get_content()
	374	dialog.Destroy()
	375	self.show_pr(self._last_out, self._last_pr, self._last_cov)
	376	else:
	377	dialog.Destroy()
	378
	379
	380	def show_pr(self, out, pr, cov=None):
	381	"""
	382	"""
	383	# Show P(r)
	384	x = pylab.arange(0.0, pr.d_max, pr.d_max / self._pr_npts)
	385
	386	y = np.zeros(len(x))
	387	dy = np.zeros(len(x))
	388	y_true = np.zeros(len(x))
	389
	390	total = 0.0
	391	pmax = 0.0
	392	cov2 = np.ascontiguousarray(cov)
	393
	394	for i in range(len(x)):
[235f514]	395	if cov2 is None:
[959eb01]	396	value = pr.pr(out, x[i])
	397	else:
	398	(value, dy[i]) = pr.pr_err(out, cov2, x[i])
	399	total += value * pr.d_max / len(x)
	400
	401	# keep track of the maximum P(r) value
	402	if value > pmax:
	403	pmax = value
	404
	405	y[i] = value
	406
	407	if self._normalize_output == True:
	408	y = y / total
	409	dy = dy / total
	410	elif self._scale_output_unity == True:
	411	y = y / pmax
	412	dy = dy / pmax
	413
[235f514]	414	if cov2 is None:
[959eb01]	415	new_plot = Data1D(x, y)
	416	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	417	else:
	418	new_plot = Data1D(x, y, dy=dy)
	419	new_plot.name = PR_FIT_LABEL
	420	new_plot.xaxis("\\rm{r}", 'A')
	421	new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
	422	new_plot.title = "P(r) fit"
	423	new_plot.id = PR_FIT_LABEL
	424	# Make sure that the plot is linear
	425	new_plot.xtransform = "x"
	426	new_plot.ytransform = "y"
	427	new_plot.group_id = GROUP_ID_PR_FIT
	428	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	429	wx.PostEvent(self.parent, NewPlotEvent(plot=new_plot, title="P(r) fit"))
	430	return x, pr.d_max
	431
	432	def load(self, data):
	433	"""
	434	Load data. This will eventually be replaced
	435	by our standard DataLoader class.
	436	"""
	437	class FileData(object):
	438	x = None
	439	y = None
	440	err = None
	441	path = None
	442
	443	def __init__(self, path):
	444	self.path = path
	445
	446	self._current_file_data = FileData(data.path)
	447
	448	# Use data loader to load file
	449	dataread = data
	450	# Notify the user if we could not read the file
	451	if dataread is None:
	452	raise RuntimeError, "Invalid data"
	453
	454	x = None
	455	y = None
	456	err = None
	457	if dataread.__class__.__name__ == 'Data1D':
	458	x = dataread.x
	459	y = dataread.y
	460	err = dataread.dy
	461	else:
	462	if isinstance(dataread, list) and len(dataread) > 0:
	463	x = dataread[0].x
	464	y = dataread[0].y
	465	err = dataread[0].dy
	466	msg = "PrView only allows a single data set at a time. "
	467	msg += "Only the first data set was loaded."
	468	wx.PostEvent(self.parent, StatusEvent(status=msg))
	469	else:
	470	if dataread is None:
	471	return x, y, err
	472	raise RuntimeError, "This tool can only read 1D data"
	473
	474	self._current_file_data.x = x
	475	self._current_file_data.y = y
	476	self._current_file_data.err = err
	477	return x, y, err
	478
	479	def load_columns(self, path="sphere_60_q0_2.txt"):
	480	"""
	481	Load 2- or 3- column ascii
	482	"""
	483	# Read the data from the data file
	484	data_x = np.zeros(0)
	485	data_y = np.zeros(0)
	486	data_err = np.zeros(0)
	487	scale = None
	488	min_err = 0.0
[ac07a3a]	489	if path is not None:
[959eb01]	490	input_f = open(path, 'r')
	491	buff = input_f.read()
	492	lines = buff.split('\n')
	493	for line in lines:
	494	try:
	495	toks = line.split()
	496	x = float(toks[0])
	497	y = float(toks[1])
	498	if len(toks) > 2:
	499	err = float(toks[2])
	500	else:
[235f514]	501	if scale is None:
[959eb01]	502	scale = 0.05 * math.sqrt(y)
	503	#scale = 0.05/math.sqrt(y)
	504	min_err = 0.01 * y
	505	err = scale * math.sqrt(y) + min_err
	506	#err = 0
	507
	508	data_x = np.append(data_x, x)
	509	data_y = np.append(data_y, y)
	510	data_err = np.append(data_err, err)
	511	except:
	512	logger.error(sys.exc_value)
	513
[ac07a3a]	514	if scale is not None:
[959eb01]	515	message = "The loaded file had no error bars, statistical errors are assumed."
	516	wx.PostEvent(self.parent, StatusEvent(status=message))
	517	else:
	518	wx.PostEvent(self.parent, StatusEvent(status=''))
	519
	520	return data_x, data_y, data_err
	521
	522	def load_abs(self, path):
	523	"""
	524	Load an IGOR .ABS reduced file
	525
	526	:param path: file path
	527
	528	:return: x, y, err vectors
	529
	530	"""
	531	# Read the data from the data file
	532	data_x = np.zeros(0)
	533	data_y = np.zeros(0)
	534	data_err = np.zeros(0)
	535	scale = None
	536	min_err = 0.0
	537
	538	data_started = False
[ac07a3a]	539	if path is not None:
[959eb01]	540	input_f = open(path, 'r')
	541	buff = input_f.read()
	542	lines = buff.split('\n')
	543	for line in lines:
	544	if data_started == True:
	545	try:
	546	toks = line.split()
	547	x = float(toks[0])
	548	y = float(toks[1])
	549	if len(toks) > 2:
	550	err = float(toks[2])
	551	else:
[235f514]	552	if scale is None:
[959eb01]	553	scale = 0.05 * math.sqrt(y)
	554	#scale = 0.05/math.sqrt(y)
	555	min_err = 0.01 * y
	556	err = scale * math.sqrt(y) + min_err
	557	#err = 0
	558
	559	data_x = np.append(data_x, x)
	560	data_y = np.append(data_y, y)
	561	data_err = np.append(data_err, err)
	562	except:
	563	logger.error(sys.exc_value)
	564	elif line.find("The 6 columns") >= 0:
	565	data_started = True
	566
[ac07a3a]	567	if scale is not None:
[959eb01]	568	message = "The loaded file had no error bars, statistical errors are assumed."
	569	wx.PostEvent(self.parent, StatusEvent(status=message))
	570	else:
	571	wx.PostEvent(self.parent, StatusEvent(status=''))
	572
	573	return data_x, data_y, data_err
	574
	575	def pr_theory(self, r, R):
	576	"""
	577	Return P(r) of a sphere for a given R
	578	For test purposes
	579	"""
	580	if r <= 2 * R:
	581	return 12.0 * ((0.5 * r / R) ** 2) * ((1.0 - 0.5 * r / R) ** 2) * (2.0 + 0.5 * r / R)
	582	else:
	583	return 0.0
	584
	585	def get_context_menu(self, plotpanel=None):
	586	"""
	587	Get the context menu items available for P(r)
	588
	589	:param graph: the Graph object to which we attach the context menu
	590
	591	:return: a list of menu items with call-back function
	592
	593	"""
	594	graph = plotpanel.graph
	595	# Look whether this Graph contains P(r) data
	596	if graph.selected_plottable not in plotpanel.plots:
	597	return []
	598	item = plotpanel.plots[graph.selected_plottable]
	599	if item.id == PR_FIT_LABEL:
	600	#add_data_hint = "Load a data file and display it on this plot"
	601	#["Add P(r) data",add_data_hint , self._on_add_data],
	602	change_n_hint = "Change the number of"
	603	change_n_hint += " points on the P(r) output"
	604	change_n_label = "Change number of P(r) points"
	605	m_list = [[change_n_label, change_n_hint, self._on_pr_npts]]
	606
	607	if self._scale_output_unity or self._normalize_output:
	608	hint = "Let the output P(r) keep the scale of the data"
	609	m_list.append(["Disable P(r) scaling", hint,
	610	self._on_disable_scaling])
	611	if not self._scale_output_unity:
	612	m_list.append(["Scale P_max(r) to unity",
	613	"Scale P(r) so that its maximum is 1",
	614	self._on_scale_unity])
	615	if not self._normalize_output:
	616	m_list.append(["Normalize P(r) to unity",
	617	"Normalize the integral of P(r) to 1",
	618	self._on_normalize])
	619
	620	return m_list
	621
	622	elif item.id in [PR_LOADED_LABEL, IQ_DATA_LABEL, IQ_FIT_LABEL, IQ_SMEARED_LABEL]:
	623	return []
	624	elif item.id == graph.selected_plottable:
	625	if issubclass(item.__class__, Data1D):
	626	return [["Compute P(r)",
	627	"Compute P(r) from distribution",
	628	self._on_context_inversion]]
	629	return []
	630
	631	def _on_disable_scaling(self, evt):
	632	"""
	633	Disable P(r) scaling
	634
	635	:param evt: Menu event
	636
	637	"""
	638	self._normalize_output = False
	639	self._scale_output_unity = False
	640	self.show_pr(self._last_out, self._last_pr, self._last_cov)
	641
	642	# Now replot the original added data
	643	for plot in self._added_plots:
	644	self._added_plots[plot].y = np.copy(self._default_Iq[plot])
	645	wx.PostEvent(self.parent,
	646	NewPlotEvent(plot=self._added_plots[plot],
	647	title=self._added_plots[plot].name,
	648	update=True))
	649
	650	# Need the update flag in the NewPlotEvent to protect against
	651	# the plot no longer being there...
	652
	653	def _on_normalize(self, evt):
	654	"""
	655	Normalize the area under the P(r) curve to 1.
	656	This operation is done for all displayed plots.
	657
	658	:param evt: Menu event
	659
	660	"""
	661	self._normalize_output = True
	662	self._scale_output_unity = False
	663
	664	self.show_pr(self._last_out, self._last_pr, self._last_cov)
	665
	666	# Now scale the added plots too
	667	for plot in self._added_plots:
	668	total = np.sum(self._added_plots[plot].y)
	669	npts = len(self._added_plots[plot].x)
	670	total *= self._added_plots[plot].x[npts - 1] / npts
	671	y = self._added_plots[plot].y / total
	672
	673	new_plot = Data1D(self._added_plots[plot].x, y)
	674	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	675	new_plot.group_id = self._added_plots[plot].group_id
	676	new_plot.id = self._added_plots[plot].id
	677	new_plot.title = self._added_plots[plot].title
	678	new_plot.name = self._added_plots[plot].name
	679	new_plot.xaxis("\\rm{r}", 'A')
	680	new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
	681	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	682	wx.PostEvent(self.parent,
	683	NewPlotEvent(plot=new_plot, update=True,
	684	title=self._added_plots[plot].name))
	685
	686	def _on_scale_unity(self, evt):
	687	"""
	688	Scale the maximum P(r) value on each displayed plot to 1.
	689
	690	:param evt: Menu event
	691
	692	"""
	693	self._scale_output_unity = True
	694	self._normalize_output = False
	695
	696	self.show_pr(self._last_out, self._last_pr, self._last_cov)
	697
	698	# Now scale the added plots too
	699	for plot in self._added_plots:
	700	_max = 0
	701	for y in self._added_plots[plot].y:
	702	if y > _max:
	703	_max = y
	704	y = self._added_plots[plot].y / _max
	705
	706	new_plot = Data1D(self._added_plots[plot].x, y)
	707	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	708	new_plot.name = self._added_plots[plot].name
	709	new_plot.xaxis("\\rm{r}", 'A')
	710	new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
	711	self.parent.update_theory(data_id=self.data_id, theory=new_plot)
	712	wx.PostEvent(self.parent,
	713	NewPlotEvent(plot=new_plot, update=True,
	714	title=self._added_plots[plot].name))
	715
	716	def start_thread(self):
	717	"""
	718	Start a calculation thread
	719	"""
	720	from pr_thread import CalcPr
	721
	722	# If a thread is already started, stop it
[7432acb]	723	if self.calc_thread is not None and self.calc_thread.isrunning():
[959eb01]	724	self.calc_thread.stop()
	725	## stop just raises the flag -- the thread is supposed to
	726	## then kill itself. In August 2014 it was shown that this is
	727	## incorrectly handled by fitting.py and a fix implemented.
	728	## It is not clear whether it is properly used here, but the
	729	## "fix" of waiting for the previous thread to end breaks the
	730	## pr perspective completely as it causes an infinite loop.
	731	## Thus it is likely the threading is bing properly handled.
	732	## While the "fix" is no longer implemented the comment is
	733	## left here till somebody ascertains that in fact the threads
	734	## are being properly handled.
	735	##
	736	## -PDB January 25, 2015
	737
	738	pr = self.pr.clone()
	739	self.calc_thread = CalcPr(pr, self.nfunc,
	740	error_func=self._thread_error,
	741	completefn=self._completed, updatefn=None)
	742	self.calc_thread.queue()
	743	self.calc_thread.ready(2.5)
	744
	745	def _thread_error(self, error):
	746	"""
	747	Call-back method for calculation errors
	748	"""
	749	wx.PostEvent(self.parent, StatusEvent(status=error))
	750
	751	def _estimate_completed(self, alpha, message, elapsed):
	752	"""
	753	Parameter estimation completed,
	754	display the results to the user
	755
	756	:param alpha: estimated best alpha
	757	:param elapsed: computation time
	758
	759	"""
	760	# Save useful info
	761	self.elapsed = elapsed
	762	self.control_panel.alpha_estimate = alpha
[ac07a3a]	763	if message is not None:
[959eb01]	764	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
	765	self.perform_estimateNT()
	766
	767	def _estimateNT_completed(self, nterms, alpha, message, elapsed):
	768	"""
	769	Parameter estimation completed,
	770	display the results to the user
	771
	772	:param alpha: estimated best alpha
	773	:param nterms: estimated number of terms
	774	:param elapsed: computation time
	775
	776	"""
	777	# Save useful info
	778	self.elapsed = elapsed
	779	self.control_panel.nterms_estimate = nterms
	780	self.control_panel.alpha_estimate = alpha
[ac07a3a]	781	if message is not None:
[959eb01]	782	wx.PostEvent(self.parent, StatusEvent(status=str(message)))
	783
	784	def _completed(self, out, cov, pr, elapsed):
	785	"""
	786	wxCallAfter Method called with the results when the inversion
	787	is done
	788
	789	:param out: output coefficient for the base functions
	790	:param cov: covariance matrix
	791	:param pr: Invertor instance
	792	:param elapsed: time spent computing
	793	"""
	794	# Ensure hat you have all inputs are ready at the time call happens:
	795	# Without CallAfter, it will freeze with wx >= 2.9.
	796	wx.CallAfter(self._completed_call, out, cov, pr, elapsed)
	797
	798	def _completed_call(self, out, cov, pr, elapsed):
	799	"""
	800	Method called with the results when the inversion
	801	is done
	802
	803	:param out: output coefficient for the base functions
	804	:param cov: covariance matrix
	805	:param pr: Invertor instance
	806	:param elapsed: time spent computing
	807
	808	"""
	809	# Save useful info
	810	self.elapsed = elapsed
	811	# Keep a copy of the last result
	812	self._last_pr = pr.clone()
	813	self._last_out = out
	814	self._last_cov = cov
	815
	816	# Save Pr invertor
	817	self.pr = pr
	818	cov = np.ascontiguousarray(cov)
	819
	820	# Show result on control panel
	821	self.control_panel.chi2 = pr.chi2
	822	self.control_panel.elapsed = elapsed
	823	self.control_panel.oscillation = pr.oscillations(out)
	824	self.control_panel.positive = pr.get_positive(out)
	825	self.control_panel.pos_err = pr.get_pos_err(out, cov)
	826	self.control_panel.rg = pr.rg(out)
	827	self.control_panel.iq0 = pr.iq0(out)
	828	self.control_panel.bck = pr.background
	829
	830	# Show I(q) fit
	831	self.show_iq(out, self.pr)
	832
	833	# Show P(r) fit
	834	self.show_pr(out, self.pr, cov)
	835
	836	def show_data(self, path=None, data=None, reset=False):
	837	"""
	838	Show data read from a file
	839
	840	:param path: file path
	841	:param reset: if True all other plottables will be cleared
	842
	843	"""
	844	#if path is not None:
	845	if data is not None:
	846	try:
	847	pr = self._create_file_pr(data)
	848	except:
	849	status = "Problem reading data: %s" % sys.exc_value
	850	wx.PostEvent(self.parent, StatusEvent(status=status))
	851	raise RuntimeError, status
	852
	853	# If the file contains nothing, just return
	854	if pr is None:
	855	raise RuntimeError, "Loaded data is invalid"
	856
	857	self.pr = pr
	858
	859	# Make a plot of I(q) data
[235f514]	860	if self.pr.err is None:
[959eb01]	861	new_plot = Data1D(self.pr.x, self.pr.y)
	862	new_plot.symbol = GUIFRAME_ID.CURVE_SYMBOL_NUM
	863	else:
	864	new_plot = Data1D(self.pr.x, self.pr.y, dy=self.pr.err)
	865	new_plot.name = IQ_DATA_LABEL
	866	new_plot.xaxis("\\rm{Q}", 'A^{-1}')
	867	new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
	868	new_plot.interactive = True
	869	new_plot.group_id = GROUP_ID_IQ_DATA
	870	new_plot.id = self.data_id
	871	new_plot.title = "I(q)"
	872	wx.PostEvent(self.parent,
	873	NewPlotEvent(plot=new_plot, title="I(q)", reset=reset))
	874
	875	self.current_plottable = new_plot
	876	# Get Q range
	877	self.control_panel.q_min = min(self.pr.x)
	878	self.control_panel.q_max = max(self.pr.x)
	879
	880	def save_data(self, filepath, prstate=None):
	881	"""
	882	Save data in provided state object.
	883
	884	:TODO: move the state code away from inversion_panel and move it here.
	885	Then remove the "prstate" input and make this method private.
	886
	887	:param filepath: path of file to write to
	888	:param prstate: P(r) inversion state
	889
	890	"""
	891	#TODO: do we need this or can we use DataLoader.loader.save directly?
	892
	893	# Add output data and coefficients to state
	894	prstate.coefficients = self._last_out
	895	prstate.covariance = self._last_cov
	896
	897	# Write the output to file
	898	# First, check that the data is of the right type
	899	if issubclass(self.current_plottable.__class__,
	900	sas.sascalc.dataloader.data_info.Data1D):
	901	self.state_reader.write(filepath, self.current_plottable, prstate)
	902	else:
	903	msg = "pr.save_data: the data being saved is not a"
	904	msg += " sas.data_info.Data1D object"
	905	raise RuntimeError, msg
	906
	907	def setup_plot_inversion(self, alpha, nfunc, d_max, q_min=None, q_max=None,
	908	bck=False, height=0, width=0):
	909	"""
	910	Set up inversion from plotted data
	911	"""
	912	self.alpha = alpha
	913	self.nfunc = nfunc
	914	self.max_length = d_max
	915	self.q_min = q_min
	916	self.q_max = q_max
	917	self.has_bck = bck
	918	self.slit_height = height
	919	self.slit_width = width
	920
	921	try:
	922	pr = self._create_plot_pr()
[ac07a3a]	923	if pr is not None:
[959eb01]	924	self.pr = pr
	925	self.perform_inversion()
	926	except:
	927	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
	928
	929	def estimate_plot_inversion(self, alpha, nfunc, d_max,
	930	q_min=None, q_max=None,
	931	bck=False, height=0, width=0):
	932	"""
	933	Estimate parameters from plotted data
	934	"""
	935	self.alpha = alpha
	936	self.nfunc = nfunc
	937	self.max_length = d_max
	938	self.q_min = q_min
	939	self.q_max = q_max
	940	self.has_bck = bck
	941	self.slit_height = height
	942	self.slit_width = width
	943
	944	try:
	945	pr = self._create_plot_pr()
[ac07a3a]	946	if pr is not None:
[959eb01]	947	self.pr = pr
	948	self.perform_estimate()
	949	except:
	950	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
	951
	952	def _create_plot_pr(self, estimate=False):
	953	"""
	954	Create and prepare invertor instance from
	955	a plottable data set.
	956
	957	:param path: path of the file to read in
	958
	959	"""
	960	# Sanity check
	961	if self.current_plottable is None:
	962	msg = "Please load a valid data set before proceeding."
	963	wx.PostEvent(self.parent, StatusEvent(status=msg))
	964	return None
	965
	966	# Get the data from the chosen data set and perform inversion
	967	pr = Invertor()
	968	pr.d_max = self.max_length
	969	pr.alpha = self.alpha
	970	pr.q_min = self.q_min
	971	pr.q_max = self.q_max
	972	pr.x = self.current_plottable.x
	973	pr.y = self.current_plottable.y
	974	pr.has_bck = self.has_bck
	975	pr.slit_height = self.slit_height
	976	pr.slit_width = self.slit_width
	977
	978	# Keep track of the plot window title to ensure that
	979	# we can overlay the plots
	980	pr.info["plot_group_id"] = self.current_plottable.group_id
	981
	982	# Fill in errors if none were provided
	983	err = self.current_plottable.dy
	984	all_zeros = True
[235f514]	985	if err is None:
[959eb01]	986	err = np.zeros(len(pr.y))
	987	else:
	988	for i in range(len(err)):
	989	if err[i] > 0:
	990	all_zeros = False
	991
	992	if all_zeros:
	993	scale = None
	994	min_err = 0.0
	995	for i in range(len(pr.y)):
	996	# Scale the error so that we can fit over several decades of Q
[235f514]	997	if scale is None:
[959eb01]	998	scale = 0.05 * math.sqrt(pr.y[i])
	999	min_err = 0.01 * pr.y[i]
	1000	err[i] = scale * math.sqrt(math.fabs(pr.y[i])) + min_err
	1001	message = "The loaded file had no error bars, "
	1002	message += "statistical errors are assumed."
	1003	wx.PostEvent(self.parent, StatusEvent(status=message))
	1004
	1005	pr.err = err
	1006
	1007	return pr
	1008
	1009	def setup_file_inversion(self, alpha, nfunc, d_max, data,
	1010	path=None, q_min=None, q_max=None,
	1011	bck=False, height=0, width=0):
	1012	"""
	1013	Set up inversion
	1014	"""
	1015	self.alpha = alpha
	1016	self.nfunc = nfunc
	1017	self.max_length = d_max
	1018	self.q_min = q_min
	1019	self.q_max = q_max
	1020	self.has_bck = bck
	1021	self.slit_height = height
	1022	self.slit_width = width
	1023
	1024	try:
	1025	pr = self._create_file_pr(data)
[ac07a3a]	1026	if pr is not None:
[959eb01]	1027	self.pr = pr
	1028	self.perform_inversion()
	1029	except:
	1030	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
	1031
	1032	def estimate_file_inversion(self, alpha, nfunc, d_max, data,
	1033	path=None, q_min=None, q_max=None,
	1034	bck=False, height=0, width=0):
	1035	"""
	1036	Estimate parameters for inversion
	1037	"""
	1038	self.alpha = alpha
	1039	self.nfunc = nfunc
	1040	self.max_length = d_max
	1041	self.q_min = q_min
	1042	self.q_max = q_max
	1043	self.has_bck = bck
	1044	self.slit_height = height
	1045	self.slit_width = width
	1046
	1047	try:
	1048	pr = self._create_file_pr(data)
	1049	if pr is not None:
	1050	self.pr = pr
	1051	self.perform_estimate()
	1052	except:
	1053	wx.PostEvent(self.parent, StatusEvent(status=sys.exc_value))
	1054
	1055	def _create_file_pr(self, data):
	1056	"""
	1057	Create and prepare invertor instance from
	1058	a file data set.
	1059
	1060	:param path: path of the file to read in
	1061
	1062	"""
	1063	# Reset the status bar so that we don't get mixed up
	1064	# with old messages.
	1065	#TODO: refactor this into a proper status handling
	1066	wx.PostEvent(self.parent, StatusEvent(status=''))
	1067	try:
	1068	class FileData(object):
	1069	x = None
	1070	y = None
	1071	err = None
	1072	path = None
	1073	def __init__(self, path):
	1074	self.path = path
	1075
	1076	self._current_file_data = FileData(data.path)
	1077	self._current_file_data.x = data.x
	1078	self._current_file_data.y = data.y
	1079	self._current_file_data.err = data.dy
	1080	x, y, err = data.x, data.y, data.dy
	1081	except:
	1082	load_error(sys.exc_value)
	1083	return None
	1084
	1085	# If the file contains no data, just return
	1086	if x is None or len(x) == 0:
	1087	load_error("The loaded file contains no data")
	1088	return None
	1089
	1090	# If we have not errors, add statistical errors
	1091	if y is not None:
[235f514]	1092	if err is None or np.all(err) == 0:
[959eb01]	1093	err = np.zeros(len(y))
	1094	scale = None
	1095	min_err = 0.0
	1096	for i in range(len(y)):
	1097	# Scale the error so that we can fit over several decades of Q
[235f514]	1098	if scale is None:
[959eb01]	1099	scale = 0.05 * math.sqrt(y[i])
	1100	min_err = 0.01 * y[i]
	1101	err[i] = scale * math.sqrt(math.fabs(y[i])) + min_err
	1102	message = "The loaded file had no error bars, "
	1103	message += "statistical errors are assumed."
	1104	wx.PostEvent(self.parent, StatusEvent(status=message))
	1105
	1106	try:
	1107	# Get the data from the chosen data set and perform inversion
	1108	pr = Invertor()
	1109	pr.d_max = self.max_length
	1110	pr.alpha = self.alpha
	1111	pr.q_min = self.q_min
	1112	pr.q_max = self.q_max
	1113	pr.x = x
	1114	pr.y = y
	1115	pr.err = err
	1116	pr.has_bck = self.has_bck
	1117	pr.slit_height = self.slit_height
	1118	pr.slit_width = self.slit_width
	1119	return pr
	1120	except:
	1121	load_error(sys.exc_value)
	1122	return None
	1123
	1124	def perform_estimate(self):
	1125	"""
	1126	Perform parameter estimation
	1127	"""
	1128	from pr_thread import EstimatePr
	1129
	1130	# If a thread is already started, stop it
[7432acb]	1131	if self.estimation_thread is not None and \
[959eb01]	1132	self.estimation_thread.isrunning():
	1133	self.estimation_thread.stop()
	1134	## stop just raises the flag -- the thread is supposed to
	1135	## then kill itself. In August 2014 it was shown that this is
	1136	## incorrectly handled by fitting.py and a fix implemented.
	1137	## It is not clear whether it is properly used here, but the
	1138	## "fix" of waiting for the previous thread to end breaks the
	1139	## pr perspective completely as it causes an infinite loop.
	1140	## Thus it is likely the threading is bing properly handled.
	1141	## While the "fix" is no longer implemented the comment is
	1142	## left here till somebody ascertains that in fact the threads
	1143	## are being properly handled.
	1144	##
	1145	## -PDB January 25, 2015
	1146	pr = self.pr.clone()
	1147	self.estimation_thread = EstimatePr(pr, self.nfunc,
	1148	error_func=self._thread_error,
	1149	completefn=self._estimate_completed,
	1150	updatefn=None)
	1151	self.estimation_thread.queue()
	1152	self.estimation_thread.ready(2.5)
	1153
	1154	def perform_estimateNT(self):
	1155	"""
	1156	Perform parameter estimation
	1157	"""
	1158	from pr_thread import EstimateNT
	1159
	1160	# If a thread is already started, stop it
[7432acb]	1161	if self.estimation_thread is not None and self.estimation_thread.isrunning():
[959eb01]	1162	self.estimation_thread.stop()
	1163	## stop just raises the flag -- the thread is supposed to
	1164	## then kill itself. In August 2014 it was shown that this is
	1165	## incorrectly handled by fitting.py and a fix implemented.
	1166	## It is not clear whether it is properly used here, but the
	1167	## "fix" of waiting for the previous thread to end breaks the
	1168	## pr perspective completely as it causes an infinite loop.
	1169	## Thus it is likely the threading is bing properly handled.
	1170	## While the "fix" is no longer implemented the comment is
	1171	## left here till somebody ascertains that in fact the threads
	1172	## are being properly handled.
	1173	##
	1174	## -PDB January 25, 2015
	1175	pr = self.pr.clone()
	1176	# Skip the slit settings for the estimation
	1177	# It slows down the application and it doesn't change the estimates
	1178	pr.slit_height = 0.0
	1179	pr.slit_width = 0.0
	1180	self.estimation_thread = EstimateNT(pr, self.nfunc,
	1181	error_func=self._thread_error,
	1182	completefn=self._estimateNT_completed,
	1183	updatefn=None)
	1184	self.estimation_thread.queue()
	1185	self.estimation_thread.ready(2.5)
	1186
	1187	def perform_inversion(self):
	1188	"""
	1189	Perform inversion
	1190	"""
	1191	self.start_thread()
	1192
	1193	def _on_context_inversion(self, event):
	1194	"""
	1195	Call-back method for plot context menu
	1196	"""
	1197	panel = event.GetEventObject()
	1198	Plugin.on_perspective(self, event=event)
	1199
	1200	# If we have more than one displayed plot, make the user choose
	1201	if len(panel.plots) >= 1 and \
	1202	panel.graph.selected_plottable in panel.plots:
	1203	dataset = panel.plots[panel.graph.selected_plottable].name
	1204	else:
	1205	logger.info("Prview Error: No data is available")
	1206	return
	1207
	1208	# Store a reference to the current plottable
	1209	# If we have a suggested value, use it.
	1210	try:
	1211	estimate = float(self.control_panel.alpha_estimate)
	1212	self.control_panel.alpha = estimate
	1213	except:
	1214	self.control_panel.alpha = self.alpha
	1215	logger.info("Prview :Alpha Not estimate yet")
	1216	try:
	1217	estimate = int(self.control_panel.nterms_estimate)
	1218	self.control_panel.nfunc = estimate
	1219	except:
	1220	self.control_panel.nfunc = self.nfunc
	1221	logger.info("Prview : ntemrs Not estimate yet")
	1222
	1223	self.current_plottable = panel.plots[panel.graph.selected_plottable]
	1224	self.set_data([self.current_plottable])
	1225	self.control_panel.plotname = dataset
	1226	#self.control_panel.nfunc = self.nfunc
	1227	self.control_panel.d_max = self.max_length
	1228	#self.parent.set_perspective(self.perspective)
	1229	self.control_panel._on_invert(None)
	1230
	1231	def get_panels(self, parent):
	1232	"""
	1233	Create and return a list of panel objects
	1234	"""
	1235	from inversion_panel import InversionControl
	1236
	1237	self.parent = parent
	1238	self.frame = MDIFrame(self.parent, None, 'None', (100, 200))
	1239	self.control_panel = InversionControl(self.frame, -1,
	1240	style=wx.RAISED_BORDER)
	1241	self.frame.set_panel(self.control_panel)
	1242	self._frame_set_helper()
	1243	self.control_panel.set_manager(self)
	1244	self.control_panel.nfunc = self.nfunc
	1245	self.control_panel.d_max = self.max_length
	1246	self.control_panel.alpha = self.alpha
	1247	self.perspective = []
	1248	self.perspective.append(self.control_panel.window_name)
	1249
	1250	return [self.control_panel]
	1251
	1252	def set_data(self, data_list=None):
	1253	"""
	1254	receive a list of data to compute pr
	1255	"""
	1256	if data_list is None:
	1257	data_list = []
	1258	if len(data_list) >= 1:
	1259	if len(data_list) == 1:
	1260	data = data_list[0]
	1261	else:
	1262	data_1d_list = []
	1263	data_2d_list = []
	1264	error_msg = ""
	1265	# separate data into data1d and data2d list
	1266	for data in data_list:
	1267	if data is not None:
	1268	if issubclass(data.__class__, Data1D):
	1269	data_1d_list.append(data)
	1270	else:
	1271	error_msg += " %s type %s \n" % (str(data.name),
	1272	str(data.__class__.__name__))
	1273	data_2d_list.append(data)
	1274	if len(data_2d_list) > 0:
	1275	msg = "PrView does not support the following data types:\n"
	1276	msg += error_msg
	1277	if len(data_1d_list) == 0:
	1278	wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
	1279	return
	1280	msg = "Prview does not allow multiple data!\n"
	1281	msg += "Please select one.\n"
	1282	if len(data_list) > 1:
	1283	from pr_widgets import DataDialog
	1284	dlg = DataDialog(data_list=data_1d_list, text=msg)
	1285	if dlg.ShowModal() == wx.ID_OK:
	1286	data = dlg.get_data()
	1287	else:
	1288	data = None
	1289	dlg.Destroy()
	1290	if data is None:
	1291	msg += "PrView receives no data. \n"
	1292	wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
	1293	return
	1294	if issubclass(data.__class__, Data1D):
	1295	try:
	1296	wx.PostEvent(self.parent,
	1297	NewPlotEvent(action='remove',
	1298	group_id=GROUP_ID_IQ_DATA,
	1299	id=self.data_id))
	1300	self.data_id = data.id
	1301	self.control_panel._change_file(evt=None, data=data)
	1302	except:
	1303	msg = "Prview Set_data: " + str(sys.exc_value)
	1304	wx.PostEvent(self.parent, StatusEvent(status=msg, info="error"))
	1305	else:
	1306	msg = "Pr cannot be computed for data of "
	1307	msg += "type %s" % (data_list[0].__class__.__name__)
	1308	wx.PostEvent(self.parent, StatusEvent(status=msg, info='error'))
	1309	else:
	1310	msg = "Pr contain no data"
	1311	wx.PostEvent(self.parent, StatusEvent(status=msg, info='warning'))
	1312
	1313	def post_init(self):
	1314	"""
	1315	Post initialization call back to close the loose ends
	1316	[Somehow openGL needs this call]
	1317	"""
	1318	pass

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/src/sas/sasgui/perspectives/pr/pr.py @ 2f85af7

Download in other formats: