Mercurial > repos > muon-spectroscopy-computational-project > larch_lcf

comparison larch_lcf.py @ 0:f59731986b61 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/MaterialsGalaxy/larch-tools/tree/main/larch_lcf commit 5be486890442dedfb327289d597e1c8110240735
author muon-spectroscopy-computational-project
date Tue, 14 Nov 2023 15:35:22 +0000
parents
children 6c28339b73f7
comparison
equal deleted inserted replaced
-1:000000000000 0:f59731986b61
1 import json
2 import sys
3
4 from common import read_group
5
6 from larch.math.lincombo_fitting import get_label, lincombo_fit
7 from larch.symboltable import Group
8
9 import matplotlib
10 import matplotlib.pyplot as plt
11
12
13 def plot(
14 group_to_fit: Group,
15 fit_group: Group,
16 energy_min: float,
17 energy_max: float,
18 ):
19 formatted_label = ""
20 for label, weight in fit_group.weights.items():
21 formatted_label += f"{label}: {weight:.3%}\n"
22
23 plt.figure()
24 plt.plot(
25 group_to_fit.energy,
26 group_to_fit.norm,
27 label=group_to_fit.filename,
28 linewidth=4,
29 color="blue",
30 )
31 plt.plot(
32 fit_group.xdata,
33 fit_group.ydata,
34 label=formatted_label[:-1],
35 linewidth=2,
36 color="orange",
37 linestyle="--",
38 )
39 plt.grid(color="black", linestyle=":", linewidth=1) # show and format grid
40 plt.xlim(energy_min, energy_max)
41 plt.xlabel("Energy (eV)")
42 plt.ylabel("normalised x$\mu$(E)") # noqa: W605
43 plt.legend()
44 plt.savefig("plot.png", format="png")
45 plt.close("all")
46
47
48 def set_label(component_group, label):
49 if label is not None:
50 component_group.filename = label
51 else:
52 component_group.filename = get_label(component_group)
53
54
55 if __name__ == "__main__":
56 # larch imports set this to an interactive backend, so need to change it
57 matplotlib.use("Agg")
58 prj_file = sys.argv[1]
59 input_values = json.load(open(sys.argv[2], "r", encoding="utf-8"))
60
61 group_to_fit = read_group(prj_file)
62 set_label(group_to_fit, input_values["label"])
63
64 component_groups = []
65 for component in input_values["components"]:
66 component_group = read_group(component["component_file"])
67 set_label(component_group, component["label"])
68 component_groups.append(component_group)
69
70 fit_group = lincombo_fit(group_to_fit, component_groups)
71 print(f"Goodness of fit (rfactor): {fit_group.rfactor:.6%}")
72
73 energy_min = input_values["energy_min"]
74 energy_max = input_values["energy_max"]
75 if input_values["energy_format"] == "relative":
76 e0 = group_to_fit.e0
77 if energy_min is not None:
78 energy_min += e0
79 if energy_max is not None:
80 energy_max += e0
81
82 plot(group_to_fit, fit_group, energy_min, energy_max)