convert.py @ 8d06779

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests

Last change on this file since 8d06779 was 8d06779, checked in by krzywon, 7 years ago
#795: Properly check for parameters during conversion.
Property mode set to `100644`
File size: 22.2 KB

Line
1	"""
2	Convert models to and from sasview.
3	"""
4	from __future__ import print_function, division
5
6	import math
7	import warnings
8
9	from .conversion_table import CONVERSION_TABLE
10	from .core import load_model_info
11
12	# List of models which SasView versions don't contain the explicit 'scale' argument.
13	# When converting such a model, please update this list.
14	MODELS_WITHOUT_SCALE = [
15	'teubner_strey',
16	'broad_peak',
17	'two_lorentzian',
18	"two_power_law",
19	'gauss_lorentz_gel',
20	'be_polyelectrolyte',
21	'correlation_length',
22	'fractal_core_shell',
23	'binary_hard_sphere',
24	'raspberry'
25	]
26
27	# List of models which SasView versions don't contain the explicit 'background' argument.
28	# When converting such a model, please update this list.
29	MODELS_WITHOUT_BACKGROUND = [
30	'guinier',
31	]
32
33	MODELS_WITHOUT_VOLFRACTION = [
34	'fractal',
35	'vesicle',
36	'multilayer_vesicle',
37	]
38
39	MAGNETIC_SASVIEW_MODELS = [
40	'core_shell',
41	'core_multi_shell',
42	'cylinder',
43	'parallelepiped',
44	'sphere',
45	]
46
47
48	# Convert new style names for polydispersity info to old style names
49	PD_DOT = [
50	("", ""),
51	("_pd", ".width"),
52	("_pd_n", ".npts"),
53	("_pd_nsigma", ".nsigmas"),
54	("_pd_type", ".type"),
55	]
56
57	CONVERT_OLD =[
58	"fittable",
59	"std",
60	"upper",
61	"lower",
62	"units"
63	]
64
65	def _rescale(par, scale):
66	return [pkscale for pk in par] if isinstance(par, list) else parscale
67
68	def _is_sld(model_info, id):
69	"""
70	Return True if parameter is a magnetic magnitude or SLD parameter.
71	"""
72	if id.startswith('M0:'):
73	return True
74	if '_pd' in id or '.' in id:
75	return False
76	for p in model_info.parameters.call_parameters:
77	if p.id == id:
78	return p.type == 'sld'
79	# check through kernel parameters in case it is a named as a vector
80	for p in model_info.parameters.kernel_parameters:
81	if p.id == id:
82	return p.type == 'sld'
83	raise ValueError("unknown parameter %r in conversion"%id)
84
85	def _rescale_sld(model_info, pars, scale):
86	"""
87	rescale all sld parameters in the new model definition by scale so the
88	numbers are nicer. Relies on the fact that all sld parameters in the
89	new model definition end with sld. For backward conversion use
90	scale=1e-6. For forward conversion use scale=1e6.
91	"""
92	return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v))
93	for id, v in pars.items())
94
95
96	def _get_translation_table(model_info):
97	_, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
98	translation = translation.copy()
99	for p in model_info.parameters.kernel_parameters:
100	if p.length > 1:
101	newid = p.id
102	oldid = translation.get(p.id, p.id)
103	translation.pop(newid, None)
104	for k in range(1, p.length+1):
105	if newid+str(k) not in translation:
106	translation[newid+str(k)] = oldid+str(k)
107	# Remove control parameter from the result
108	if model_info.control:
109	translation[model_info.control] = "CONTROL"
110	return translation
111
112	# ========= FORWARD CONVERSION sasview 3.x => sasmodels ===========
113	def _dot_pd_to_underscore_pd(par):
114	if par.endswith(".width"):
115	return par[:-6]+"_pd"
116	elif par.endswith(".type"):
117	return par[:-5]+"_pd_type"
118	elif par.endswith(".nsigmas"):
119	return par[:-8]+"_pd_nsigma"
120	elif par.endswith(".npts"):
121	return par[:-5]+"_pd_n"
122	else:
123	return par
124
125	def _pd_to_underscores(pars):
126	return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items())
127
128	def _convert_name(conv_dict, pars):
129	"""
130	Renames parameter values (upper, lower, etc) to v4.0 names
131	:param conv_dict: conversion dictionary mapping new name : old name
132	:param pars: parameters to convert
133	:return:
134	"""
135	import re
136	new_pars = {}
137	for key_par, value_par in pars.iteritems():
138	for key_conv, value_conv in conv_dict.iteritems():
139	if re.search(value_conv, key_par):
140	new_pars[key_par.replace(value_conv, key_conv)] = value_par
141	break
142	elif re.search("background", key_par) or re.search("scale", key_par):
143	new_pars[key_par] = value_par
144	break
145	return new_pars
146
147	def _convert_pars(pars, mapping):
148	"""
149	Rename the parameters and any associated polydispersity attributes.
150	"""
151	newpars = pars.copy()
152	for new, old in mapping.items():
153	if old == new: continue
154	if old is None: continue
155	for underscore, dot in PD_DOT:
156	source = old+dot
157	if source in newpars:
158	if new is not None:
159	target = new+dot
160	else:
161	target = None
162	if source != target:
163	if target:
164	newpars[target] = pars[old+dot]
165	del newpars[source]
166	return newpars
167
168
169	def _conversion_target(model_name):
170	"""
171	Find the sasmodel name which translates into the sasview name.
172
173	Note: CoreShellEllipsoidModel translates into core_shell_ellipsoid:1.
174	This is necessary since there is only one variant in sasmodels for the
175	two variants in sasview.
176	"""
177	for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items():
178	if sasview_name == model_name:
179	return sasmodels_name
180	return None
181
182
183	def _hand_convert(name, oldpars):
184	if name == 'core_shell_parallelepiped':
185	# Make sure pd on rim parameters defaults to zero
186	# ... probably not necessary.
187	oldpars['rimA.width'] = 0.0
188	oldpars['rimB.width'] = 0.0
189	oldpars['rimC.width'] = 0.0
190	elif name == 'core_shell_ellipsoid:1':
191	# Reverse translation (from new to old), from core_shell_ellipsoid.c
192	# equat_shell = equat_core + thick_shell
193	# polar_core = equat_core * x_core
194	# polar_shell = equat_core * x_core + thick_shell*x_polar_shell
195	# Forward translation (from old to new), inverting reverse translation:
196	# thick_shell = equat_shell - equat_core
197	# x_core = polar_core / equat_core
198	# x_polar_shell = (polar_shell - polar_core)/(equat_shell - equat_core)
199	# Auto translation (old <=> new) happens after hand_convert
200	# equat_shell <=> thick_shell
201	# polar_core <=> x_core
202	# polar_shell <=> x_polar_shell
203	# So...
204	equat_core, equat_shell = oldpars['equat_core'], oldpars['equat_shell']
205	polar_core, polar_shell = oldpars['polar_core'], oldpars['polar_shell']
206	oldpars['equat_shell'] = equat_shell - equat_core
207	oldpars['polar_core'] = polar_core / equat_core
208	oldpars['polar_shell'] = (polar_shell-polar_core)/(equat_shell-equat_core)
209	elif name == 'hollow_cylinder':
210	# now uses radius and thickness
211	thickness = oldpars['radius'] - oldpars['core_radius']
212	pd = oldpars['radius.width']*oldpars['radius']/thickness
213	oldpars['radius'] = thickness
214	oldpars['radius.width'] = pd
215	elif name == 'pearl_necklace':
216	pass
217	#_remove_pd(oldpars, 'num_pearls', name)
218	#_remove_pd(oldpars, 'thick_string', name)
219	elif name == 'polymer_micelle':
220	if 'ndensity' in oldpars:
221	oldpars['ndensity'] /= 1e15
222	if 'ndensity.lower' in oldpars:
223	oldpars['ndensity.lower'] /= 1e15
224	if 'ndensity.upper' in oldpars:
225	oldpars['ndensity.upper'] /= 1e15
226	elif name == 'rpa':
227	# convert scattering lengths from femtometers to centimeters
228	for p in "L1", "L2", "L3", "L4":
229	if p in oldpars:
230	oldpars[p] /= 1e-13
231	if p + ".lower" in oldpars:
232	oldpars[p + ".lower"] /= 1e-13
233	if p + ".upper" in oldpars:
234	oldpars[p + ".upper"] /= 1e-13
235	elif name == 'spherical_sld':
236	oldpars["CONTROL"] += 1
237	elif name == 'teubner_strey':
238	# basically undoing the entire Teubner-Strey calculations here.
239	# drho = (sld_a - sld_b)
240	# k = 2.0math.pixi/d
241	# a2 = (1.0 + k2)2
242	# c1 = 2.0 * xi*2 (1.0 - k**2)
243	# c2 = xi**4
244	# prefactor = 8.0math.piphi(1.0-phi)drho*2c2/xi
245	# scale = 1e-4*prefactor
246	# oldpars['scale'] = a2/scale
247	# oldpars['c1'] = c1/scale
248	# oldpars['c2'] = c2/scale
249
250	# need xi, d, sld_a, sld_b, phi=volfraction_a
251	# assume contrast is 1.0e-6, scale=1, background=0
252	sld_a, sld_b = 1.0, 0.
253	drho = sld_a - sld_b
254
255	# find xi
256	p_scale = oldpars['scale']
257	p_c1 = oldpars['c1']
258	p_c2= oldpars['c2']
259	xi = math.sqrt(2/(math.sqrt(p_scale/p_c2) + 0.5*p_c1/p_c2))
260
261	# find d from xi
262	k = math.sqrt(1 - 0.5p_c1/p_c2xi**2)
263	d = 2math.pixi/k
264
265	# solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2)
266	# favour volume fraction in [0, 0.5]
267	c = xi / (1e-4 * 8.0 * math.pi * drho*2 p_c2)
268	phi = 0.5 - math.sqrt(0.25 - c)
269
270	# scale sld_a by 1e-6 because the translator will scale it back
271	oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6,
272	sld_b=sld_b, scale=1.0)
273	oldpars.pop('c1')
274	oldpars.pop('c2')
275
276	return oldpars
277
278	def convert_model(name, pars, use_underscore=False):
279	"""
280	Convert model from old style parameter names to new style.
281	"""
282	newname = _conversion_target(name)
283	if newname is None:
284	return name, pars
285	if ':' in newname: # core_shell_ellipsoid:1
286	model_info = load_model_info(newname[:-2])
287	# Know that the table exists and isn't multiplicity so grab it directly
288	# Can't use _get_translation_table since that will return the 'bare'
289	# version.
290	translation = CONVERSION_TABLE[newname]
291	else:
292	model_info = load_model_info(newname)
293	translation = _get_translation_table(model_info)
294	newpars = _convert_name(translation, pars.copy())
295	newpars = _hand_convert(newname, newpars)
296	newpars = _convert_pars(newpars, translation)
297	newpars = _rescale_sld(model_info, newpars, 1e6)
298	newpars.setdefault('scale', 1.0)
299	newpars.setdefault('background', 0.0)
300	if use_underscore:
301	newpars = _pd_to_underscores(newpars)
302	return newname, newpars
303
304
305	# ========= BACKWARD CONVERSION sasmodels => sasview 3.x ===========
306
307	def _revert_pars(pars, mapping):
308	"""
309	Rename the parameters and any associated polydispersity attributes.
310	"""
311	newpars = pars.copy()
312
313	for new, old in mapping.items():
314	for underscore, dot in PD_DOT:
315	if old and old+underscore == new+dot:
316	continue
317	if new+underscore in newpars:
318	if old is not None:
319	newpars[old+dot] = pars[new+underscore]
320	del newpars[new+underscore]
321	for k in list(newpars.keys()):
322	for underscore, dot in PD_DOT[1:]: # skip "" => ""
323	if k.endswith(underscore):
324	newpars[k[:-len(underscore)]+dot] = newpars[k]
325	del newpars[k]
326	return newpars
327
328	def revert_name(model_info):
329	oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])
330	return oldname
331
332	def _remove_pd(pars, key, name):
333	"""
334	Remove polydispersity from the parameter list.
335
336	Note: operates in place
337	"""
338	# Bumps style parameter names
339	width = pars.pop(key+".width", 0.0)
340	n_points = pars.pop(key+".npts", 0)
341	if width != 0.0 and n_points != 0:
342	warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name))
343	pars.pop(key+".nsigmas", None)
344	pars.pop(key+".type", None)
345	return pars
346
347	def _trim_vectors(model_info, pars, oldpars):
348	_, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
349	for p in model_info.parameters.kernel_parameters:
350	if p.length_control is not None:
351	n = int(pars[p.length_control])
352	oldname = translation.get(p.id, p.id)
353	for k in range(n+1, p.length+1):
354	for _, old in PD_DOT:
355	oldpars.pop(oldname+str(k)+old, None)
356	return oldpars
357
358	def revert_pars(model_info, pars):
359	"""
360	Convert model from new style parameter names to old style.
361	"""
362	if model_info.composition is not None:
363	composition_type, parts = model_info.composition
364	if composition_type == 'product':
365	translation = _get_translation_table(parts[0])
366	# structure factor models include scale:scale_factor mapping
367	translation.update(_get_translation_table(parts[1]))
368	else:
369	raise NotImplementedError("cannot convert to sasview sum")
370	else:
371	translation = _get_translation_table(model_info)
372	oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation)
373	oldpars = _trim_vectors(model_info, pars, oldpars)
374
375	# Make sure the control parameter is an integer
376	if "CONTROL" in oldpars:
377	oldpars["CONTROL"] = int(oldpars["CONTROL"])
378
379	# Note: update compare.constrain_pars to match
380	name = model_info.id
381	if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
382	if oldpars.pop('scale', 1.0) != 1.0:
383	warnings.warn("parameter scale not used in sasview %s"%name)
384	if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
385	if oldpars.pop('background', 0.0) != 0.0:
386	warnings.warn("parameter background not used in sasview %s"%name)
387
388	# Remove magnetic parameters from non-magnetic sasview models
389	if name not in MAGNETIC_SASVIEW_MODELS:
390	oldpars = dict((k, v) for k, v in oldpars.items() if ':' not in k)
391
392	# If it is a product model P*S, then check the individual forms for special
393	# cases. Note: despite the structure factor alone not having scale or
394	# background, the product model does, so this is below the test for
395	# models without scale or background.
396	namelist = name.split('') if '' in name else [name]
397	for name in namelist:
398	if name in MODELS_WITHOUT_VOLFRACTION:
399	del oldpars['volfraction']
400	elif name == 'core_multi_shell':
401	# kill extra shells
402	for k in range(5, 11):
403	oldpars.pop('sld_shell'+str(k), 0)
404	oldpars.pop('thick_shell'+str(k), 0)
405	oldpars.pop('mtheta:sld'+str(k), 0)
406	oldpars.pop('mphi:sld'+str(k), 0)
407	oldpars.pop('M0:sld'+str(k), 0)
408	_remove_pd(oldpars, 'sld_shell'+str(k), 'sld')
409	_remove_pd(oldpars, 'thick_shell'+str(k), 'thickness')
410	elif name == 'core_shell_parallelepiped':
411	_remove_pd(oldpars, 'rimA', name)
412	_remove_pd(oldpars, 'rimB', name)
413	_remove_pd(oldpars, 'rimC', name)
414	elif name == 'hollow_cylinder':
415	# now uses radius and thickness
416	thickness = oldpars['core_radius']
417	oldpars['radius'] += thickness
418	oldpars['radius.width'] *= thickness/oldpars['radius']
419	#elif name in ['mono_gauss_coil', 'poly_gauss_coil']:
420	# del oldpars['i_zero']
421	elif name == 'onion':
422	oldpars.pop('n_shells', None)
423	elif name == 'pearl_necklace':
424	_remove_pd(oldpars, 'num_pearls', name)
425	_remove_pd(oldpars, 'thick_string', name)
426	elif name == 'polymer_micelle':
427	if 'ndensity' in oldpars:
428	oldpars['ndensity'] *= 1e15
429	elif name == 'rpa':
430	# convert scattering lengths from femtometers to centimeters
431	for p in "L1", "L2", "L3", "L4":
432	if p in oldpars: oldpars[p] *= 1e-13
433	if pars['case_num'] < 2:
434	for k in ("a", "b"):
435	for p in ("L", "N", "Phi", "b", "v"):
436	oldpars.pop(p+k, None)
437	for k in "Kab,Kac,Kad,Kbc,Kbd".split(','):
438	oldpars.pop(k, None)
439	elif pars['case_num'] < 5:
440	for k in ("a",):
441	for p in ("L", "N", "Phi", "b", "v"):
442	oldpars.pop(p+k, None)
443	for k in "Kab,Kac,Kad".split(','):
444	oldpars.pop(k, None)
445	elif name == 'spherical_sld':
446	oldpars["CONTROL"] -= 1
447	# remove polydispersity from shells
448	for k in range(1, 11):
449	_remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')
450	_remove_pd(oldpars, 'thick_inter'+str(k), 'interface')
451	# remove extra shells
452	for k in range(int(pars['n_shells']), 11):
453	oldpars.pop('sld_flat'+str(k), 0)
454	oldpars.pop('thick_flat'+str(k), 0)
455	oldpars.pop('thick_inter'+str(k), 0)
456	oldpars.pop('func_inter'+str(k), 0)
457	oldpars.pop('nu_inter'+str(k), 0)
458	elif name == 'stacked_disks':
459	_remove_pd(oldpars, 'n_stacking', name)
460	elif name == 'teubner_strey':
461	# basically redoing the entire Teubner-Strey calculations here.
462	volfraction = oldpars.pop('volfraction_a')
463	xi = oldpars.pop('xi')
464	d = oldpars.pop('d')
465	sld_a = oldpars.pop('sld_a')
466	sld_b = oldpars.pop('sld_b')
467	drho = 1e6*(sld_a - sld_b) # conversion autoscaled these
468	k = 2.0math.pixi/d
469	a2 = (1.0 + k2)2
470	c1 = 2.0 * xi*2 (1.0 - k**2)
471	c2 = xi**4
472	prefactor = 8.0math.pivolfraction(1.0-volfraction)drho*2c2/xi
473	scale = 1e-4*prefactor
474	oldpars['scale'] = a2/scale
475	oldpars['c1'] = c1/scale
476	oldpars['c2'] = c2/scale
477
478	#print("convert from",list(sorted(pars)))
479	#print("convert to",list(sorted(oldpars.items())))
480	return oldpars
481
482	def constrain_new_to_old(model_info, pars):
483	"""
484	Restrict parameter values to those that will match sasview.
485	"""
486	name = model_info.id
487	# Note: update convert.revert_model to match
488	if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
489	pars['scale'] = 1
490	if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
491	pars['background'] = 0
492	# sasview multiplies background by structure factor
493	if '*' in name:
494	pars['background'] = 0
495
496	# Shut off magnetism when comparing non-magnetic sasview models
497	if name not in MAGNETIC_SASVIEW_MODELS:
498	suppress_magnetism = False
499	for key in pars.keys():
500	if key.startswith("M0:"):
501	suppress_magnetism = suppress_magnetism or (pars[key] != 0)
502	pars[key] = 0
503	if suppress_magnetism:
504	warnings.warn("suppressing magnetism for comparison with sasview")
505
506	# Shut off theta polydispersity since algorithm has changed
507	if 'theta_pd_n' in pars:
508	if pars['theta_pd_n'] != 0:
509	warnings.warn("suppressing theta polydispersity for comparison with sasview")
510	pars['theta_pd_n'] = 0
511
512	# If it is a product model P*S, then check the individual forms for special
513	# cases. Note: despite the structure factor alone not having scale or
514	# background, the product model does, so this is below the test for
515	# models without scale or background.
516	namelist = name.split('') if '' in name else [name]
517	for name in namelist:
518	if name in MODELS_WITHOUT_VOLFRACTION:
519	pars['volfraction'] = 1
520	if name == 'core_multi_shell':
521	pars['n'] = min(math.ceil(pars['n']), 4)
522	elif name == 'gel_fit':
523	pars['scale'] = 1
524	elif name == 'line':
525	pars['scale'] = 1
526	pars['background'] = 0
527	elif name == 'mono_gauss_coil':
528	pars['scale'] = 1
529	elif name == 'onion':
530	pars['n_shells'] = math.ceil(pars['n_shells'])
531	elif name == 'pearl_necklace':
532	pars['string_thickness_pd_n'] = 0
533	pars['number_of_pearls_pd_n'] = 0
534	elif name == 'poly_gauss_coil':
535	pars['scale'] = 1
536	elif name == 'rpa':
537	pars['case_num'] = int(pars['case_num'])
538	elif name == 'spherical_sld':
539	pars['n_shells'] = math.ceil(pars['n_shells'])
540	pars['n_steps'] = math.ceil(pars['n_steps'])
541	for k in range(1, 11):
542	pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5)
543	for k in range(2, 11):
544	pars['thickness%d_pd_n'%k] = 0
545	pars['interface%d_pd_n'%k] = 0
546	elif name == 'teubner_strey':
547	pars['scale'] = 1
548	if pars['volfraction_a'] > 0.5:
549	pars['volfraction_a'] = 1.0 - pars['volfraction_a']
550	elif name == 'unified_power_Rg':
551	pars['level'] = int(pars['level'])
552
553	def _check_one(name, seed=None):
554	"""
555	Generate a random set of parameters for name, and check that they can
556	be converted back to SasView 3.x and forward again to sasmodels. Raises
557	an error if the parameters are changed.
558	"""
559	from . import compare
560
561	model_info = load_model_info(name)
562
563	old_name = revert_name(model_info)
564	if old_name is None:
565	return
566
567	pars = compare.get_pars(model_info, use_demo=False)
568	pars = compare.randomize_pars(model_info, pars, seed=seed)
569	if name == "teubner_strey":
570	# T-S model is underconstrained, so fix the assumptions.
571	pars['sld_a'], pars['sld_b'] = 1.0, 0.0
572	compare.constrain_pars(model_info, pars)
573	constrain_new_to_old(model_info, pars)
574	old_pars = revert_pars(model_info, pars)
575	new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True)
576	if 1:
577	print("==== %s in ====="%name)
578	print(str(compare.parlist(model_info, pars, True)))
579	print("==== %s ====="%old_name)
580	for k, v in sorted(old_pars.items()):
581	print(k, v)
582	print("==== %s out ====="%new_name)
583	print(str(compare.parlist(model_info, new_pars, True)))
584	assert name==new_name, "%r != %r"%(name, new_name)
585	for k, v in new_pars.items():
586	assert k in pars, "%s: %r appeared from conversion"%(name, k)
587	if isinstance(v, float):
588	assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r %s != %s"%(name, k, v, pars[k])
589	else:
590	assert v == pars[k], "%s: %r %s != %s"%(name, k, v, pars[k])
591	for k, v in pars.items():
592	assert k in pars, "%s: %r not converted"%(name, k)
593
594	def test_backward_forward():
595	from .core import list_models
596	for name in list_models('all'):
597	L = lambda: _check_one(name, seed=1)
598	L.description = name
599	yield L

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

SasView

source: sasmodels/sasmodels/convert.py @ 8d06779

Download in other formats: