convert.py @ 037aa3d

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests

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

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

SasView

source: sasmodels/sasmodels/convert.py @ 037aa3d

Download in other formats: