forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
fitting.py
1600 lines (1419 loc) · 71.2 KB
/
fitting.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
# NScD Oak Ridge National Laboratory, European Spallation Source,
# Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
from mantid.api import IFunction, AlgorithmManager, mtd
from mantid.kernel import logger
from mantid.simpleapi import CalculateChiSquared, FunctionFactory, plotSpectrum
from .function import PeaksFunction, PhysicalProperties, ResolutionModel, Background, Function
from .energies import energies
from .normalisation import split2range, ionname2Nre
from .CrystalFieldMultiSite import CrystalFieldMultiSite
from scipy.constants import physical_constants
import numpy as np
import re
import scipy.optimize as sp
import warnings
# RegEx pattern matching a composite function parameter name, eg f2.Sigma.
FN_PATTERN = re.compile('f(\\d+)\\.(.+)')
# RegEx pattern matching a composite function parameter name, eg f2.Sigma. Multi-spectrum case.
FN_MS_PATTERN = re.compile('f(\\d+)\\.f(\\d+)\\.(.+)')
CONSTRAINTS_PATTERN = re.compile(r'constraints=\((.*?)\)')
FWHM_PATTERN = re.compile(r'FWHM[X|Y]\d+=\(\),')
PHYSICAL_PROPERTIES_PATTERN = re.compile(r'(name=.*?,)(.*?)(PhysicalProperties=\(.*?\),)')
TEMPERATURES_PATTERN = re.compile(r'(name=.*?,)(.*?)(Temperatures=\(.*?\),)')
TIES_PATTERN = re.compile(r',ties=\((.*?)\)')
def makeWorkspace(xArray, yArray):
"""Create a workspace that doesn't appear in the ADS"""
alg = AlgorithmManager.createUnmanaged('CreateWorkspace')
alg.initialize()
alg.setChild(True)
alg.setProperty('DataX', xArray)
alg.setProperty('DataY', yArray)
alg.setProperty('OutputWorkspace', 'dummy')
alg.execute()
return alg.getProperty('OutputWorkspace').value
def islistlike(arg):
return (not hasattr(arg, "strip")) and (hasattr(arg, "__getitem__") or hasattr(arg, "__iter__")) and hasattr(arg, "__len__")
def cfpstrmaker(x, pref='B'):
return [pref+str(k)+str(x) for k in [2, 4, 6] if x <= k]
def getSymmAllowedParam(sym_str):
if 'T' in sym_str or 'O' in sym_str:
return ['B40', 'B60']
if any([sym_str == val for val in ['C1', 'Ci']]):
return sum([cfpstrmaker(i) for i in range(7)]
+ [cfpstrmaker(i, 'IB') for i in range(1, 7)],[])
retval = cfpstrmaker(0)
if '6' in sym_str or '3' in sym_str:
retval += cfpstrmaker(6)
if any([sym_str == val for val in ['C6', 'C3h', 'C6h']]):
retval += cfpstrmaker(6, 'IB')
if ('3' in sym_str and '3h' not in sym_str) or 'S6' in sym_str:
retval += cfpstrmaker(3)
if any([sym_str == val for val in ['C3', 'S6']]):
retval += cfpstrmaker(3, 'IB') + cfpstrmaker(6, 'IB')
if '4' in sym_str or '2' in sym_str:
retval += cfpstrmaker(4)
if any([sym_str == val for val in ['C4', 'S4', 'C4h']]):
retval += cfpstrmaker(4, 'IB')
if ('2' in sym_str and '2d' not in sym_str) or 'Cs' in sym_str:
retval += cfpstrmaker(2)
if any([sym_str == val for val in ['C2', 'Cs', 'C2h']]):
retval += cfpstrmaker(2, 'IB') + cfpstrmaker(4, 'IB')
return retval
#pylint: disable=too-many-instance-attributes,too-many-public-methods
class CrystalField(object):
"""Calculates the crystal fields for one ion"""
allowed_symmetries = ['C1', 'Ci', 'C2', 'Cs', 'C2h', 'C2v', 'D2', 'D2h', 'C4', 'S4', 'C4h',
'D4', 'C4v', 'D2d', 'D4h', 'C3', 'S6', 'D3', 'C3v', 'D3d', 'C6', 'C3h',
'C6h', 'D6', 'C6v', 'D3h', 'D6h', 'T', 'Td', 'Th', 'O', 'Oh']
lande_g = [6.0 / 7., 4.0 / 5., 8.0 / 11., 3.0 / 5., 2.0 / 7., 0.0, 2.0,
3.0 / 2., 4.0 / 3., 5.0 / 4., 6.0 / 5., 7.0 / 6., 8.0 / 7.]
default_peakShape = 'Gaussian'
default_background = 'FlatBackground'
default_spectrum_size = 200
field_parameter_names = ['BmolX','BmolY','BmolZ','BextX','BextY','BextZ',
'B20','B21','B22','B40','B41','B42','B43','B44','B60','B61','B62','B63','B64','B65','B66',
'IB21','IB22','IB41','IB42','IB43','IB44','IB61','IB62','IB63','IB64','IB65','IB66']
def __init__(self, Ion, Symmetry, **kwargs):
"""
Constructor.
@param Ion: A rare earth ion. Possible values:
Ce, Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb
@param Symmetry: Symmetry of the field. Possible values:
C1, Ci, C2, Cs, C2h, C2v, D2, D2h, C4, S4, C4h, D4, C4v, D2d, D4h, C3,
S6, D3, C3v, D3d, C6, C3h, C6h, D6, C6v, D3h, D6h, T, Td, Th, O, Oh
@param kwargs: Other field parameters and attributes. Acceptable values include:
ToleranceEnergy: energy tolerance,
ToleranceIntensity: intensity tolerance,
ResolutionModel: A resolution model.
FWHMVariation: Absolute value of allowed variation of a peak width during a fit.
FixAllPeaks: A boolean flag that fixes all parameters of the peaks.
Field parameters:
BmolX: The x-component of the molecular field,
BmolY: The y-component of the molecular field,
BmolZ: The z-component of the molecular field,
BextX: The x-component of the external field,
BextY: The y-component of the external field,
BextZ: The z-component of the external field,
B20: Real part of the B20 field parameter,
B21: Real part of the B21 field parameter,
B22: Real part of the B22 field parameter,
B40: Real part of the B40 field parameter,
B41: Real part of the B41 field parameter,
B42: Real part of the B42 field parameter,
B43: Real part of the B43 field parameter,
B44: Real part of the B44 field parameter,
B60: Real part of the B60 field parameter,
B61: Real part of the B61 field parameter,
B62: Real part of the B62 field parameter,
B63: Real part of the B63 field parameter,
B64: Real part of the B64 field parameter,
B65: Real part of the B65 field parameter,
B66: Real part of the B66 field parameter,
IB21: Imaginary part of the B21 field parameter,
IB22: Imaginary part of the B22 field parameter,
IB41: Imaginary part of the B41 field parameter,
IB42: Imaginary part of the B42 field parameter,
IB43: Imaginary part of the B43 field parameter,
IB44: Imaginary part of the B44 field parameter,
IB61: Imaginary part of the B61 field parameter,
IB62: Imaginary part of the B62 field parameter,
IB63: Imaginary part of the B63 field parameter,
IB64: Imaginary part of the B64 field parameter,
IB65: Imaginary part of the B65 field parameter,
IB66: Imaginary part of the B66 field parameter,
Each of the following parameters can be either a single float or an array of floats.
They are either all floats or all arrays of the same size.
IntensityScaling: A scaling factor for the intensity of each spectrum.
FWHM: A default value for the full width at half maximum of the peaks.
Temperature: A temperature "of the spectrum" in Kelvin
PhysicalProperty: A list of PhysicalProperties objects denoting the required data type
Note that physical properties datasets should follow inelastic spectra
See the Crystal Field Python Interface help page for more details.
"""
self._background = None
if 'Temperature' in kwargs:
temperature = kwargs['Temperature']
del kwargs['Temperature']
else:
temperature = -1
# Create self.function attribute
self._makeFunction(Ion, Symmetry, temperature)
self.Temperature = temperature
self.Ion = Ion
self.Symmetry = Symmetry
self._resolutionModel = None
self._physprop = None
free_parameters = {key: kwargs[key] for key in kwargs if key in CrystalField.field_parameter_names}
if 'ResolutionModel' in kwargs and 'FWHM' in kwargs:
msg = 'Both ''ResolutionModel'' and ''FWHM'' specified but can only accept one width option.'
msg += ' Preferring to use ResolutionModel, and ignoring FWHM.'
kwargs.pop('FWHM')
warnings.warn(msg, SyntaxWarning)
for key in kwargs:
if key == 'ToleranceEnergy':
self.ToleranceEnergy = kwargs[key]
elif key == 'ToleranceIntensity':
self.ToleranceIntensity = kwargs[key]
elif key == 'IntensityScaling':
self.IntensityScaling = kwargs[key]
elif key == 'FWHM':
self.FWHM = kwargs[key]
elif key == 'ResolutionModel':
self.ResolutionModel = kwargs[key]
elif key == 'NPeaks':
self.NPeaks = kwargs[key]
elif key == 'FWHMVariation':
self.FWHMVariation = kwargs[key]
elif key == 'FixAllPeaks':
self.FixAllPeaks = kwargs[key]
elif key == 'PhysicalProperty':
self.PhysicalProperty = kwargs[key]
elif key not in free_parameters:
raise RuntimeError('Unknown attribute/parameters %s' % key)
# Cubic is a special case where B44=5*B40, B64=-21*B60
is_cubic = self.Symmetry.startswith('T') or self.Symmetry.startswith('O')
symm_allowed_par = getSymmAllowedParam(self.Symmetry)
for param in CrystalField.field_parameter_names:
if param in free_parameters:
self.function.setParameter(param, free_parameters[param])
if is_cubic and (param == 'B44' or param == 'B64'):
continue
if param not in symm_allowed_par:
self.function.fixParameter(param)
else:
self.function.freeParameter(param)
self._setPeaks()
# Required to build the target function for the first time (which includes applying all ties)
self.crystalFieldFunction.initialize()
# Eigensystem
self._dirty_eigensystem = True
self._eigenvalues = None
self._eigenvectors = None
self._hamiltonian = None
# Peak lists
self._dirty_peaks = True
self._peakList = None
# Spectra
self._plot_window = {}
# self._setDefaultTies()
self.chi2 = None
def _makeFunction(self, ion, symmetry, temperature):
if temperature is not None and islistlike(temperature) and len(temperature) > 1:
self.function = FunctionFactory.createFunction('CrystalFieldMultiSpectrum')
self._isMultiSpectrum = True
tempStr = 'Temperatures'
else:
self.function = FunctionFactory.createFunction('CrystalFieldSpectrum')
self._isMultiSpectrum = False
tempStr = 'Temperature'
self.function.setAttributeValue('Ion', ion)
self.function.setAttributeValue('Symmetry', symmetry)
if temperature:
temperature = [float(val) for val in temperature] if islistlike(temperature) else float(temperature)
self.function.setAttributeValue(tempStr, temperature)
def _remakeFunction(self, temperature):
"""Redefines the internal function, e.g. when `Temperature` (number of datasets) change"""
fieldParams = self._getFieldParameters()
self._makeFunction(self.Ion, self.Symmetry, temperature)
for item in fieldParams.items():
self.function.setParameter(item[0], item[1])
for param in CrystalField.field_parameter_names:
if param not in fieldParams.keys():
self.function.fixParameter(param)
def _setPeaks(self):
if self._isMultiSpectrum:
self._peaks = []
# Even a single spectrum when used in conjunction with physical properties
# is treated as a multi-spectra and the number of spectra includes the number
# of physical properties. However, only spectra have the peaks.
len_physical_properties = 0
if self._physprop is not None:
if isinstance(self._physprop, list):
len_physical_properties = len(self._physprop)
else:
len_physical_properties = 1
for i in range(self.NumberOfSpectra - len_physical_properties):
self._peaks.append(PeaksFunction(self.crystalFieldFunction, 'f%s.' % i, 1))
else:
self._peaks = PeaksFunction(self.crystalFieldFunction, '', 0)
@property
def crystalFieldFunction(self):
if not self._isMultiSpectrum and self.background is not None:
return self.function[1]
else:
return self.function
def makePeaksFunction(self, i):
"""Form a definition string for the CrystalFieldPeaks function
@param i: Index of a spectrum.
"""
temperature = self._getTemperature(i)
out = 'name=CrystalFieldPeaks,Ion=%s,Symmetry=%s,Temperature=%s' % (self.Ion, self.Symmetry, temperature)
out += ',ToleranceEnergy=%s,ToleranceIntensity=%s' % (self.ToleranceEnergy, self.ToleranceIntensity)
out += ',%s' % ','.join(['%s=%s' % item for item in self._getFieldParameters().items()])
return out
def makeSpectrumFunction(self, i=0):
"""Form a definition string for the CrystalFieldSpectrum function
@param i: Index of a spectrum.
"""
if not self._isMultiSpectrum:
return str(self.function)
else:
funs = self.function.createEquivalentFunctions()
return str(funs[i])
def makePhysicalPropertiesFunction(self, i=0):
"""Form a definition string for one of the crystal field physical properties functions
@param i: Index of the dataset (default=0), or a PhysicalProperties object.
"""
if hasattr(i, 'toString'):
out = i.toString()
ppobj = i
else:
if self._physprop is None:
raise RuntimeError('Physical properties environment not defined.')
ppobj = self._physprop[i] if islistlike(self._physprop) else self._physprop
if hasattr(ppobj, 'toString'):
out = ppobj.toString()
else:
return ''
out += ',Ion=%s,Symmetry=%s' % (self.Ion, self.Symmetry)
fieldParams = self._getFieldParameters()
if len(fieldParams) > 0:
out += ',%s' % ','.join(['%s=%s' % item for item in fieldParams.items()])
ties = self._getFieldTies()
if ppobj.TypeID == PhysicalProperties.SUSCEPTIBILITY:
ties += ',' if ties else ''
ties += 'Lambda=0' if ppobj.Lambda == 0. else ''
ties += ',Chi0=0' if ppobj.Chi0 == 0. else ''
if len(ties) > 0:
out += ',ties=(%s)' % ties
constraints = self._getFieldConstraints()
if len(constraints) > 0:
out += ',constraints=(%s)' % constraints
return out
def makeMultiSpectrumFunction(self):
fun = re.sub(FWHM_PATTERN, '', str(self.function))
fun = re.sub(TEMPERATURES_PATTERN, r'\1\3\2', fun)
fun = re.sub(PHYSICAL_PROPERTIES_PATTERN, r'\1\3\2', fun)
return fun
@property
def Ion(self):
"""Get value of Ion attribute. For example:
cf = CrystalField(...)
...
ion = cf.Ion
"""
return self.crystalFieldFunction.getAttributeValue('Ion')
@Ion.setter
def Ion(self, value):
"""Set new value of Ion attribute. For example:
cf = CrystalField(...)
...
cf.Ion = 'Pr'
"""
self._nre = ionname2Nre(value)
self.crystalFieldFunction.setAttributeValue('Ion', value)
self._dirty_eigensystem = True
self._dirty_peaks = True
@property
def Symmetry(self):
"""Get value of Symmetry attribute. For example:
cf = CrystalField(...)
...
symm = cf.Symmetry
"""
return self.crystalFieldFunction.getAttributeValue('Symmetry')
@Symmetry.setter
def Symmetry(self, value):
"""Set new value of Symmetry attribute. For example:
cf = CrystalField(...)
...
cf.Symmetry = 'Td'
"""
if value not in self.allowed_symmetries:
msg = 'Value %s is not allowed for attribute Symmetry.\nList of allowed values: %s' % \
(value, ', '.join(self.allowed_symmetries))
raise RuntimeError(msg)
self.crystalFieldFunction.setAttributeValue('Symmetry', value)
self._dirty_eigensystem = True
self._dirty_peaks = True
self.crystalFieldFunction.initialize()
@property
def ToleranceEnergy(self):
"""Get energy tolerance"""
return self.crystalFieldFunction.getAttributeValue('ToleranceEnergy')
@ToleranceEnergy.setter
def ToleranceEnergy(self, value):
"""Set energy tolerance"""
self.crystalFieldFunction.setAttributeValue('ToleranceEnergy', float(value))
self._dirty_peaks = True
@property
def ToleranceIntensity(self):
"""Get intensity tolerance"""
return self.crystalFieldFunction.getAttributeValue('ToleranceIntensity')
@ToleranceIntensity.setter
def ToleranceIntensity(self, value):
"""Set intensity tolerance"""
self.crystalFieldFunction.setAttributeValue('ToleranceIntensity', float(value))
self._dirty_peaks = True
@property
def IntensityScaling(self):
if not self._isMultiSpectrum:
return self.crystalFieldFunction.getParameterValue('IntensityScaling')
iscaling = []
for i in range(self.NumberOfSpectra):
paramName = 'IntensityScaling%s' % i
iscaling.append(self.crystalFieldFunction.getParameterValue(paramName))
return iscaling
@IntensityScaling.setter
def IntensityScaling(self, value):
if not self._isMultiSpectrum:
if islistlike(value):
if len(value) == 1:
value = value[0]
else:
raise ValueError('IntensityScaling is expected to be a single floating point value')
self.crystalFieldFunction.setParameter('IntensityScaling', value)
else:
n = self.NumberOfSpectra
if not islistlike(value) or len(value) != n:
raise ValueError('IntensityScaling is expected to be a list of %s values' % n)
for i in range(n):
paramName = 'IntensityScaling%s' % i
self.crystalFieldFunction.setParameter(paramName, value[i])
self._dirty_peaks = True
@property
def Temperature(self):
attrName = 'Temperatures' if self._isMultiSpectrum else 'Temperature'
return self.crystalFieldFunction.getAttributeValue(attrName)
@Temperature.setter
def Temperature(self, value):
if islistlike(value) and len(value) == 1:
value = value[0]
if self._isMultiSpectrum:
if not islistlike(value):
# Try to keep current set of field parameters.
self._remakeFunction(float(value))
return
self.crystalFieldFunction.setAttributeValue('Temperatures', value)
else:
if islistlike(value):
self._remakeFunction(value)
return
self.crystalFieldFunction.setAttributeValue('Temperature', float(value))
self._dirty_peaks = True
@property
def FWHM(self):
attrName = 'FWHMs' if self._isMultiSpectrum else 'FWHM'
fwhm = self.crystalFieldFunction.getAttributeValue(attrName)
if self._isMultiSpectrum:
nDatasets = len(self.Temperature)
if len(fwhm) != nDatasets:
return list(fwhm) * nDatasets
return fwhm
@FWHM.setter
def FWHM(self, value):
if islistlike(value) and len(value) == 1:
value = value[0]
if self._isMultiSpectrum:
if not islistlike(value):
value = [value] * self.NumberOfSpectra
if len(value) != len(self.Temperature):
if self.PhysicalProperty is not None and len(value) == len(self.Temperature) - len(self.PhysicalProperty):
value = value + [0] * len(self.PhysicalProperty)
# Cast all types to match the first elem so we don't have mixed lists of int/doubles
value = [float(v) for v in value]
else:
raise RuntimeError('Vector of FWHMs must either have same size as '
'Temperatures (%i) or have size 1.' % (len(self.Temperature)))
self.crystalFieldFunction.setAttributeValue('FWHMs', value)
else:
if islistlike(value):
raise ValueError('FWHM is expected to be a single floating point value')
self.crystalFieldFunction.setAttributeValue('FWHM', float(value))
# If both FWHM and ResolutionModel is set, may cause runtime errors
self._resolutionModel = None
if self._isMultiSpectrum:
for i in range(self.NumberOfSpectra):
if self.crystalFieldFunction.getAttributeValue('FWHMX%s' % i):
self.crystalFieldFunction.setAttributeValue('FWHMX%s' % i, [])
if self.crystalFieldFunction.getAttributeValue('FWHMY%s' % i):
self.crystalFieldFunction.setAttributeValue('FWHMY%s' % i, [])
else:
if self.crystalFieldFunction.getAttributeValue('FWHMX'):
self.crystalFieldFunction.setAttributeValue('FWHMX', [])
if self.crystalFieldFunction.getAttributeValue('FWHMY'):
self.crystalFieldFunction.setAttributeValue('FWHMY', [])
@property
def FWHMVariation(self):
return self.crystalFieldFunction.getAttributeValue('FWHMVariation')
@FWHMVariation.setter
def FWHMVariation(self, value):
self.crystalFieldFunction.setAttributeValue('FWHMVariation', float(value))
def __getitem__(self, item):
return self.crystalFieldFunction.getParameterValue(item)
def __setitem__(self, key, value):
self.crystalFieldFunction.setParameter(key, value)
@property
def ResolutionModel(self):
return self._resolutionModel
@ResolutionModel.setter
def ResolutionModel(self, value):
if hasattr(value, 'model'):
self._resolutionModel = value
else:
self._resolutionModel = ResolutionModel(value)
if self._isMultiSpectrum:
NumberOfPhysProp = len(self._physprop) if islistlike(self._physprop) else (0 if self._physprop is None else 1)
NSpec = self.NumberOfSpectra - NumberOfPhysProp
if not self._resolutionModel.multi or self._resolutionModel.NumberOfSpectra != NSpec:
raise RuntimeError('Resolution model is expected to have %s functions, found %s' %
(NSpec, self._resolutionModel.NumberOfSpectra))
for i in range(self._resolutionModel.NumberOfSpectra):
model = self._resolutionModel.model[i]
self.crystalFieldFunction.setAttributeValue('FWHMX%s' % i, model[0])
self.crystalFieldFunction.setAttributeValue('FWHMY%s' % i, model[1])
else:
model = self._resolutionModel.model
self.crystalFieldFunction.setAttributeValue('FWHMX', model[0])
self.crystalFieldFunction.setAttributeValue('FWHMY', model[1])
# If FWHM is set, it overrides resolution model, so unset it
if self._isMultiSpectrum and any(self.crystalFieldFunction.getAttributeValue('FWHMs')):
self.crystalFieldFunction.setAttributeValue('FWHMs', [0.] * self.NumberOfSpectra)
elif not self._isMultiSpectrum and self.crystalFieldFunction.getAttributeValue('FWHM'):
self.crystalFieldFunction.setAttributeValue('FWHM', 0.)
@property
def FixAllPeaks(self):
return self.crystalFieldFunction.getAttributeValue('FixAllPeaks')
@FixAllPeaks.setter
def FixAllPeaks(self, value):
self.crystalFieldFunction.setAttributeValue('FixAllPeaks', value)
@property
def PeakShape(self):
return self.crystalFieldFunction.getAttributeValue('PeakShape')
@PeakShape.setter
def PeakShape(self, value):
self.crystalFieldFunction.setAttributeValue('PeakShape', value)
@property
def NumberOfSpectra(self):
return self.crystalFieldFunction.getNumberDomains()
@property
def NPeaks(self):
return self.crystalFieldFunction.getAttributeValue('NPeaks')
@NPeaks.setter
def NPeaks(self, value):
self.crystalFieldFunction.setAttributeValue('NPeaks', value)
@property
def peaks(self):
return self._peaks
@property
def background(self):
return self._background
@background.setter
def background(self, value):
"""
Define the background function.
Args:
value: an instance of function.Background class or a list of instances
in a multi-spectrum case
"""
if self._background is not None:
raise ValueError('Background has been set already')
if not hasattr(value, 'toString'):
raise TypeError('Expected a Background object, found %s' % str(value))
if not self._isMultiSpectrum:
fun_str = value.toString() + ';' + str(self.function)
self.function = FunctionFactory.createInitialized(fun_str)
self._background = self._makeBackgroundObject(value)
self._setPeaks()
else:
self.function.setAttributeValue("Background", value.toString())
self._background = []
for ispec in range(self.NumberOfSpectra):
prefix = 'f%s.' % ispec
self._background.append(self._makeBackgroundObject(value, prefix))
def _makeBackgroundObject(self, value, prefix=''):
if len(value.functions) > 1:
prefix += 'f0.'
n_functions = 0
peak, background = None, None
if value.peak is not None:
peak = Function(self.function, prefix=prefix + f'f{n_functions}.')
n_functions += 1
if value.background is not None:
background = Function(self.function, prefix=prefix + f'f{n_functions}.')
n_functions += 1
other_functions = []
for function_index in range(n_functions, len(value.functions)):
other_functions.append(Function(self.function, prefix=prefix + f'f{function_index}.'))
return Background(peak=peak, background=background, functions=other_functions)
@property
def PhysicalProperty(self):
return self._physprop
@PhysicalProperty.setter
def PhysicalProperty(self, value):
vlist = value if islistlike(value) else [value]
if all([hasattr(pp, 'TypeID') for pp in vlist]):
nOldPP = len(self._physprop) if islistlike(self._physprop) else (0 if self._physprop is None else 1)
self._physprop = value
else:
errmsg = 'PhysicalProperty input must be a PhysicalProperties'
errmsg += ' instance or a list of such instances'
raise ValueError(errmsg)
# If a spectrum (temperature) is already defined, or multiple physical properties
# given, redefine the CrystalFieldMultiSpectrum function.
if not self.isPhysicalPropertyOnly or islistlike(self.PhysicalProperty):
tt = self.Temperature if islistlike(self.Temperature) else [self.Temperature]
ww = list(self.FWHM) if islistlike(self.FWHM) else [self.FWHM]
# Last n-set of temperatures correspond to PhysicalProperties
if nOldPP > 0:
tt = tt[:-nOldPP]
# Removes 'negative' temperature, which is a flag for no INS dataset
tt = [val for val in tt if val > 0]
pptt = [0 if val.Temperature is None else val.Temperature for val in vlist]
self._remakeFunction(list(tt)+pptt)
self.FWHM = ww
self._setPeaks()
ppids = [pp.TypeID for pp in vlist]
self.function.setAttributeValue('PhysicalProperties', [0]*len(tt)+ppids)
for attribs in [pp.getAttributes(i+len(tt)) for i, pp in enumerate(vlist)]:
for item in attribs.items():
if 'Lambda' in item[0] or 'Chi0' in item[0]:
self.function.setParameter(item[0], item[1])
if item[1] == 0.:
self.function.tie(item[0], '0.')
else:
self.function.removeTie(item[0])
else:
self.function.setAttributeValue(item[0], item[1])
@property
def isPhysicalPropertyOnly(self):
return (not islistlike(self.Temperature) and self.Temperature < 0
and self.PhysicalProperty is not None)
def ties(self, **kwargs):
"""Set ties on the field parameters.
@param kwargs: Ties as name=value pairs: name is a parameter name,
the value is a tie string or a number. For example:
tie(B20 = 0.1, IB23 = '2*B23')
"""
for tie in kwargs:
self.crystalFieldFunction.tie(tie, str(kwargs[tie]))
def constraints(self, *args):
"""
Set constraints for the field parameters.
@param args: A list of constraints. For example:
constraints('B00 > 0', '0.1 < B43 < 0.9')
"""
self.crystalFieldFunction.addConstraints(','.join(args))
def getEigenvalues(self):
self._calcEigensystem()
return self._eigenvalues
def getEigenvectors(self):
self._calcEigensystem()
return self._eigenvectors
def getHamiltonian(self):
self._calcEigensystem()
return self._hamiltonian
def getPeakList(self, i=0):
"""Get the peak list for spectrum i as a numpy array"""
self._calcPeaksList(i)
peaks = np.array([self._peakList.column(0), self._peakList.column(1)])
return peaks
def getSpectrum(self, i=0, workspace=None, ws_index=0):
"""
Get the i-th spectrum calculated with the current field and peak parameters.
Alternatively can be called getSpectrum(workspace, ws_index). Spectrum index i is assumed zero.
Examples:
cf.getSpectrum() # Return the first spectrum calculated on a generated set of x-values.
cf.getSpectrum(1, ws, 5) # Calculate the second spectrum using the x-values from the 6th spectrum
# in workspace ws.
cf.getSpectrum(ws) # Calculate the first spectrum using the x-values from the 1st spectrum
# in workspace ws.
cf.getSpectrum(ws, 3) # Calculate the first spectrum using the x-values from the 4th spectrum
# in workspace ws.
@param i: Index of a spectrum to get.
@param workspace: A workspace to base on. If not given the x-values of the output spectrum will be
generated.
@param ws_index: An index of a spectrum from workspace to use.
@return: A tuple of (x, y) arrays
"""
wksp = workspace
# Allow to call getSpectrum with a workspace as the first argument.
if not isinstance(i, int):
if wksp is not None:
if not isinstance(wksp, int):
raise RuntimeError('Spectrum index is expected to be int. Got %s' % i.__class__.__name__)
ws_index = wksp
wksp = i
i = 0
if (self.Temperature[i] if islistlike(self.Temperature) else self.Temperature) < 0:
raise RuntimeError('You must first define a temperature for the spectrum')
# Workspace is given, always calculate
if wksp is None:
xArray = None
elif isinstance(wksp, list) or isinstance(wksp, np.ndarray):
xArray = wksp
else:
return self._calcSpectrum(i, wksp, ws_index)
if xArray is None:
x_min, x_max = self.calc_xmin_xmax(i)
xArray = np.linspace(x_min, x_max, self.default_spectrum_size)
yArray = np.zeros_like(xArray)
wksp = makeWorkspace(xArray, yArray)
return self._calcSpectrum(i, wksp, 0)
def getHeatCapacity(self, workspace=None, ws_index=0):
"""
Get the heat cacpacity calculated with the current crystal field parameters
Examples:
cf.getHeatCapacity() # Returns the heat capacity from 1 < T < 300 K in 1 K steps
cf.getHeatCapacity(ws) # Returns the heat capacity with temperatures given by ws.
cf.getHeatCapacity(ws, ws_index) # Use the spectrum indicated by ws_index for x-values
@param workspace: Either a Mantid workspace whose x-values will be used as the temperatures
to calculate the heat capacity; or a list of numpy ndarray of temperatures.
Temperatures are in Kelvin.
@param ws_index: The index of a spectrum in workspace to use (default=0).
"""
return self._getPhysProp(PhysicalProperties('Cv'), workspace, ws_index)
def getSusceptibility(self, *args, **kwargs):
"""
Get the magnetic susceptibility calculated with the current crystal field parameters.
The susceptibility is calculated using Van Vleck's formula (2nd order perturbation theory)
Examples:
cf.getSusceptibility() # Returns the susceptibility || [001] for 1<T<300 K in 1 K steps
cf.getSusceptibility(T) # Returns the susceptibility with temperatures given by T.
cf.getSusceptibility(ws, 0) # Use x-axis of spectrum 0 of workspace ws as temperature
cf.getSusceptibility(T, [1, 1, 1]) # Returns the susceptibility along [111].
cf.getSusceptibility(T, 'powder') # Returns the powder averaged susceptibility
cf.getSusceptibility(T, 'cgs') # Returns the susceptibility || [001] in cgs normalisation
cf.getSusceptibility(..., Inverse=True) # Calculates the inverse susceptibility instead
cf.getSusceptibility(Temperature=ws, ws_index=0, Hdir=[1, 1, 0], Unit='SI', Inverse=True)
@param Temperature: Either a Mantid workspace whose x-values will be used as the temperatures
to calculate the heat capacity; or a list or numpy ndarray of temperatures.
Temperatures are in Kelvin.
@param ws_index: The index of a spectrum to use (default=0) if using a workspace for x.
@param Hdir: The magnetic field direction to calculate the susceptibility along. Either a
Cartesian vector with z along the quantisation axis of the CF parameters, or the
string 'powder' (case insensitive) to get the powder averaged susceptibility
default: [0, 0, 1]
@param Unit: Any one of the strings 'bohr', 'SI' or 'cgs' (case insensitive) to indicate whether
to output in atomic (bohr magneton/Tesla/ion), SI (m^3/mol) or cgs (cm^3/mol) units.
default: 'cgs'
@param Inverse: Whether to calculate the susceptibility (Inverse=False, default) or inverse
susceptibility (Inverse=True).
"""
# Sets defaults / parses keyword arguments
workspace = kwargs['Temperature'] if 'Temperature' in kwargs.keys() else None
ws_index = kwargs['ws_index'] if 'ws_index' in kwargs.keys() else 0
# Parses argument list
args = list(args)
if len(args) > 0:
workspace = args.pop(0)
if 'mantid' in str(type(workspace)) and len(args) > 0:
ws_index = args.pop(0)
# _calcSpectrum updates parameters and susceptibility has a 'Lambda' parameter which other
# CF functions don't have. This causes problems if you want to calculate another quantity after
x, y = self._getPhysProp(PhysicalProperties('chi', *args, **kwargs), workspace, ws_index)
return x, y
def getMagneticMoment(self, *args, **kwargs):
"""
Get the magnetic moment calculated with the current crystal field parameters.
The moment is calculated by adding a Zeeman term to the CF Hamiltonian and then diagonlising
the result. This function calculates either M(H) [default] or M(T), but can only calculate
a 1D function (e.g. not M(H,T) simultaneously).
Examples:
cf.getMagneticMoment() # Returns M(H) for H||[001] from 0 to 30 T in 0.1 T steps
cf.getMagneticMoment(H) # Returns M(H) for H||[001] at specified values of H (in Tesla)
cf.getMagneticMoment(ws, 0) # Use x-axis of spectrum 0 of ws as applied field magnitude.
cf.getMagneticMoment(H, [1, 1, 1]) # Returns the magnetic moment along [111].
cf.getMagneticMoment(H, 'powder') # Returns the powder averaged M(H)
cf.getMagneticMoment(H, 'cgs') # Returns the moment || [001] in cgs units (emu/mol)
cf.getMagneticMoment(Temperature=T) # Returns M(T) for H=1T || [001] at specified T (in K)
cf.getMagneticMoment(10, [1, 1, 0], Temperature=T) # Returns M(T) for H=10T || [110].
cf.getMagneticMoment(..., Inverse=True) # Calculates 1/M instead (keyword only)
cf.getMagneticMoment(Hmag=ws, ws_index=0, Hdir=[1, 1, 0], Unit='SI', Temperature=T, Inverse=True)
@param Hmag: The magnitude of the applied magnetic field in Tesla, specified either as a Mantid
workspace whose x-values will be used; or a list or numpy ndarray of field points.
If Temperature is specified as a list / array / workspace, Hmag must be scalar.
(default: 0-30T in 0.1T steps, or 1T if temperature vector specified)
@param Temperature: The temperature in Kelvin at which to calculate the moment.
Temperature is a keyword argument only. Can be a list, ndarray or workspace.
If Hmag is a list / array / workspace, Temperature must be scalar.
(default=1K)
@param ws_index: The index of a spectrum to use (default=0) if using a workspace for x.
@param Hdir: The magnetic field direction to calculate the susceptibility along. Either a
Cartesian vector with z along the quantisation axis of the CF parameters, or the
string 'powder' (case insensitive) to get the powder averaged susceptibility
default: [0, 0, 1]
@param Unit: Any one of the strings 'bohr', 'SI' or 'cgs' (case insensitive) to indicate whether
to output in atomic (bohr magneton/ion), SI (Am^2/mol) or cgs (emu/mol) units.
default: 'bohr'
@param Inverse: Whether to calculate the susceptibility (Inverse=False, default) or inverse
susceptibility (Inverse=True). Inverse is a keyword argument only.
"""
# Sets defaults / parses keyword arguments
workspace = None
ws_index = kwargs['ws_index'] if 'ws_index' in kwargs.keys() else 0
hmag = kwargs['Hmag'] if 'Hmag' in kwargs.keys() else 1.
temperature = kwargs['Temperature'] if 'Temperature' in kwargs.keys() else 1.
# Checks whether to calculate M(H) or M(T)
hmag_isscalar = (not islistlike(hmag) or len(hmag) == 1)
hmag_isvector = (islistlike(hmag) and len(hmag) > 1)
t_isscalar = (not islistlike(temperature) or len(temperature) == 1)
t_isvector = (islistlike(temperature) and len(temperature) > 1)
if hmag_isscalar and (t_isvector or 'mantid' in str(type(temperature))):
typeid = 4
workspace = temperature
kwargs['Hmag'] = hmag[0] if islistlike(hmag) else hmag
else:
typeid = 3
if t_isscalar and (hmag_isvector or 'mantid' in str(type(hmag))):
workspace = hmag
kwargs['Temperature'] = temperature[0] if islistlike(temperature) else temperature
# Parses argument list
args = list(args)
if len(args) > 0:
if typeid == 4:
kwargs['Hmag'] = args.pop(0)
else:
workspace = args.pop(0)
if 'mantid' in str(type(workspace)) and len(args) > 0:
ws_index = args.pop(0)
pptype = 'M(T)' if (typeid == 4) else 'M(H)'
self._typeid = self._str2id(typeid) if isinstance(typeid, str) else int(typeid)
return self._getPhysProp(PhysicalProperties(pptype, *args, **kwargs), workspace, ws_index)
def _calc_gJuB(self):
gj = 2.0 if (self._nre < 1) else self.lande_g[self._nre - 1]
gJuB = gj * physical_constants['Bohr magneton in eV/T'][0] * 1000.
return gJuB
def getDipoleMatrixComponent(self, nComponent, gJuB = None):
self._calcEigensystem() #will not recalculate if already called (unless _dirty_eigensystem)
if gJuB is None:
gJuB = self._calc_gJuB()
if nComponent == 'X' or nComponent == 'x':
_, _, h_n = energies(self._nre, BextX=1.0)
elif nComponent == 'Y' or nComponent == 'y':
_, _, h_n = energies(self._nre, BextY=1.0)
elif nComponent == 'Z' or nComponent == 'z':
_, _, h_n = energies(self._nre, BextZ=1.0)
else:
raise Exception('Invalid Argument, nComponent must be: X, Y or Z (case insensitive)')
i_n = np.dot(np.conj(np.transpose(self._eigenvectors)), np.dot(h_n, self._eigenvectors))
return np.multiply(i_n, np.conj(i_n))/(gJuB ** 2)
def getDipoleMatrix(self):
"""Returns the dipole transition matrix as a numpy array"""
gJuB = self._calc_gJuB()
trans = self.getDipoleMatrixComponent('X', gJuB) + self.getDipoleMatrixComponent('Y', gJuB) \
+ self.getDipoleMatrixComponent('Z', gJuB)
return trans
def plot(self, i=0, workspace=None, ws_index=0, name=None):
"""Plot a spectrum. Parameters are the same as in getSpectrum(...)"""
createWS = AlgorithmManager.createUnmanaged('CreateWorkspace')
createWS.initialize()
xArray, yArray = self.getSpectrum(i, workspace, ws_index)
ws_name = name if name is not None else 'CrystalField_%s' % self.Ion
if isinstance(i, int):
if workspace is None:
if i > 0:
ws_name += '_%s' % i
createWS.setProperty('DataX', xArray)
createWS.setProperty('DataY', yArray)
createWS.setProperty('OutputWorkspace', ws_name)
createWS.execute()
plot_window = self._plot_window[i] if i in self._plot_window else None
self._plot_window[i] = plotSpectrum(ws_name, 0, window=plot_window, clearWindow=True)
else:
ws_name += '_%s' % workspace
if i > 0:
ws_name += '_%s' % i
createWS.setProperty('DataX', xArray)
createWS.setProperty('DataY', yArray)
createWS.setProperty('OutputWorkspace', ws_name)
createWS.execute()
plotSpectrum(ws_name, 0)
else:
ws_name += '_%s' % i
createWS.setProperty('DataX', xArray)
createWS.setProperty('DataY', yArray)
createWS.setProperty('OutputWorkspace', ws_name)
createWS.execute()
plotSpectrum(ws_name, 0)
def update(self, func, index=0):
"""
Update values of the fitting parameters.
@param func: A IFunction object containing new parameter values.
@param index: The index of the function to update in the Background object.
"""
self.function = func
if self._background is not None:
if isinstance(self._background, list):
for background in self._background:
background.update(func, index)
else:
self._background.update(func, index)
self._setPeaks()
def calc_xmin_xmax(self, i):
"""Calculate the x-range containing interesting features of a spectrum (for plotting)
@param i: If an integer is given then calculate the x-range for the i-th spectrum.
If None given get the range covering all the spectra.
@return: Tuple (xmin, xmax)
"""
peaks = self.getPeakList(i)
x_min = np.min(peaks[0])
x_max = np.max(peaks[0])