convert.py @ 47fafe9

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests

Last change on this file since 47fafe9 was 47fafe9, checked in by krzywon, 7 years ago
#795: Old style structure factors are now loaded in properly.
Property mode set to `100644`
File size: 22.5 KB

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

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

SasView

source: sasmodels/sasmodels/convert.py @ 47fafe9

Download in other formats: