<jupyter_start><jupyter_text>A brief demonstration of the double spike toolboxThis is a quick guide to the main features of the double spike toolbox for python. The package uses the numpy, scipy, and matplotlib libraries.<jupyter_code>import doublespike as ds
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>IsoDataA python object called `IsoData` is used to store data on the different isotope systems. It can be invoked with a string giving the element name. The object is initialised with a number of sensible default values, but these can then be altered by the user. Information in `IsoData` includes the isotope numbers, the atomic masses, the standard composition, and the default error model coefficients. The default information stored about Fe is shown below.<jupyter_code>isodata_fe = ds.IsoData("Fe")
print(isodata_fe)<jupyter_output>element: Fe
 isonum: [54 56 57 58]
 mass: [53.9396148 55.9349421 56.9353987 57.9332805]
 standard: [0.05845 0.91754 0.02119 0.00282]
 rawspike: [[9.83719674e-01 1.55003100e-02 7.00014000e-04 8.00016000e-05]
 [6.99300699e-04 9.98801199e-01 3.99600400e-04 9.99000999e-05]
 [1.00000000e-03 6.97000000e-02 9.28800000e-01 5.00000000e-04]
 [3.90000000e-03 1.37800000e-01 1.25000000e-02 8.45800000e-01]]
 spike: None
 isoinv: [54 56 57 58]
 errormodel: {'measured': {'type': 'fixed-total', 'intensity': 10.0, 'a': array([2.0709756e-10, 2.0709756e-10, 2.0709756e-10, 2.0709756e-10]), 'b': array([2.00272058e-09, 2.00272058e-09, 2.00272058e-09, 2.00272058e-09]), 'c': array([0., 0., 0., 0.])}, 'spike': {'type': 'fixed-total', 'intensity': 10.0, 'a': array([0., 0., 0., 0.]), 'b': array([0., 0., 0., 0.]), 'c': array([0., 0., 0., 0.])}, 'standard': {'type': 'fixed-total', 'intensity': 10.0, 'a': array([0., 0., 0., 0.]), 'b': array([0., 0., 0., 0.]), 'c': array([0., 0., 0., 0.])}}<jupyter_text>The stored data can be accessed and set through the object attributes. For example, shown below for Fe are the isotope numbers, standard composition, the atomic masses, and the third and fourth Oak Ridge National Labs (ORNL) spike compositions (corresponding to the 57Fe and 58Fe spikes). Note that all values are expressed as proportions of each isotope, rather than as isotopic ratios. Note also that python indexing starts from 0, so that the index 2 corresponds to the third spike.<jupyter_code>print(isodata_fe.isonum)
print(isodata_fe.standard)
print(isodata_fe.mass)
print(isodata_fe.rawspike[2, :])  # third ORNL single spike composition
print(isodata_fe.rawspike[3, :])  # fourth ORNL single spike composition<jupyter_output>[54 56 57 58]
[0.05845 0.91754 0.02119 0.00282]
[53.9396148 55.9349421 56.9353987 57.9332805]
[1.000e-03 6.970e-02 9.288e-01 5.000e-04]
[0.0039 0.1378 0.0125 0.8458]<jupyter_text>2D Error surfaces 1The function `errorcurve2d` plots a 2D error surface as a contour plot, showing how the error varies with both double spike to sample mixture proportions and the proportions in which the two single spikes are mixed to make the double spike. The example below is for a 57Fe-58Fe double spike, using pure spikes.<jupyter_code>ds.errorcurve2d(isodata_fe, "pure", [57, 58])<jupyter_output><empty_output><jupyter_text>2D Error surfaces 2Calculations can also be performed using real spikes rather than hypothetical pure spikes. Below is an example using the third and fourth ORNL spikes (corresponding to 57Fe and 58Fe spikes). This produces Figure 2 in the manuscript. Spike purity usually has only a little effect.<jupyter_code>ds.errorcurve2d(isodata_fe, "real", [2, 3])<jupyter_output><empty_output><jupyter_text>2D Error surfaces 3Some isotopes have more than four isotopes, and in these cases the isotopes to use in the inversion must be specified. Here is an example for Ca, with a 42Ca-48Ca double spike and 40Ca, 42Ca, 44Ca, 48Ca used in inversion. The second and sixth ORNL spikes are the 42Ca and 48Ca spikes. This produces Figure 5 in the manuscript.<jupyter_code>isodata_ca = ds.IsoData("Ca")
print("spike 2:", isodata_ca.rawspike[1, :])
print("spike 6:", isodata_ca.rawspike[5, :])
ds.errorcurve2d(isodata_ca, "real", [1, 5], [40, 42, 44, 48])<jupyter_output>spike 2: [4.87951205e-02 9.44805519e-01 6.99930007e-04 5.49945005e-03
 9.99900010e-05 9.99900010e-05]
spike 6: [2.09966405e-02 2.39961606e-04 9.99840026e-05 7.19884818e-04
 9.99840026e-05 9.77843545e-01]<jupyter_text>Error curves 1The function `errorcurve` plots the error in either the fractionation factor $\alpha$ or a chosen ratio as a function of the double spike to sample proportion. Note again that all compositions are specified by proportions of each isotope rather than by ratios. Here we look at the error curve for a double spike which is 50% 57Fe and 50% 58Fe.<jupyter_code>spike = [0, 0, 0.5, 0.5]
ds.errorcurve(isodata_fe, spike)<jupyter_output><empty_output><jupyter_text>Error curves 2The function `errorcurve2` plots the error in either the fractionation factor alpha or a chosen ratio as a function of the proportion of the two single spikes that make the double spike. The proportion in which double spike and sample are mixed must be specified, as must the single spikes to use. We give an example here for a pure 57Fe-58Fe double spike with 50% double spike to 50% sample.<jupyter_code>ds.errorcurve2(isodata_fe, "pure", 0.5, [57, 58])<jupyter_output><empty_output><jupyter_text>Optimal spikes 1The function `optimalspike` finds double spike compositions which minimise the error on alpha or a chosen ratio. This can be done either for pure spikes, or with the real spikes in the `rawspike` attribute. The following example finds the best 57Fe-58Fe spike using the real spikes available from ORNL. The 3rd and 4th spikes correspond to the 57Fe and 58Fe spikes. The optimal double spike turns out to be quite close to a 50-50 mix of the available spikes (optspikeprop). The actual double spike compositions are in optspike, the optimal double spike-sample mixing proportions in optprop, the error estimates in opterr, rescaled error estimates in optppmperamu, and the isotopes that were used in the inversion in optisoinv.<jupyter_code>ds.optimalspike(isodata_fe, "real", [2, 3])<jupyter_output><empty_output><jupyter_text>Optimal spikes 2If the isotopes to spike are not specified, the `optimalspike` function checks all possible combinations. An example for Fe pure spikes is shown below. This produces Table 1 in the manuscript.<jupyter_code>ds.optimalspike(isodata_fe, "pure")<jupyter_output><empty_output><jupyter_text>Optimal spikes 3An example for Fe ORNL spikes is shown below. This produces Table 2 in the manuscript.<jupyter_code>ds.optimalspike(isodata_fe, "real")<jupyter_output><empty_output><jupyter_text>Optimal spikes 4By default, `optimalspike` minimises the error on alpha, but for radiogenic work we often wish to minimise the error on a particular ratio. An example of this is Pb. Shown below is the result of minimising the error on 206Pb/204Pb. This produces part of Table 3 in the manuscript.<jupyter_code>isodata_pb = ds.IsoData("Pb")
ds.optimalspike(isodata_pb, "pure", errorratio=[206, 204])<jupyter_output><empty_output><jupyter_text>Optimal spikes 5For elements with more than four isotopes, such as Ca, `optimalspike` tries all combinations of four isotopes in the inversion. This produces Table 4 in the manuscript. We only display the first 31 rows of the optimal spike composition.<jupyter_code>os = ds.optimalspike(isodata_ca, "pure")
print(os["optspike"][0:31, :])<jupyter_output>[[0.         0.39639725 0.         0.         0.         0.60360275]
 [0.         0.46424839 0.         0.         0.53575161 0.        ]
 [0.         0.         0.43534086 0.         0.         0.56465914]
 [0.         0.         0.48983764 0.         0.51016236 0.        ]
 [0.         0.         0.         0.16050323 0.         0.83949677]
 [0.         0.         0.         0.15500241 0.         0.84499759]
 [0.         0.         0.         0.15407932 0.         0.84592068]
 [0.         0.         0.         0.         0.47249243 0.52750757]
 [0.         0.         0.         0.22975212 0.77024788 0.        ]
 [0.         0.         0.         0.22380252 0.77619748 0.        ]
 [0.         0.         0.         0.21621117 0.78378883 0.        ]
 [0.         0.48265354 0.51734646 0.         0.         0.        ]
 [0.         0.2126147  0.         0.         0.7873853  0.        ]
 [0.         0.         0.37656792 0.         0.62343208 0.        ]
 [0.         0.29328132 0.70671868[...]<jupyter_text>Optimal error curves 1The function `errorcurveoptimalspike` calculates the `optimalspike` and then plots the `errorcurve`.<jupyter_code>ds.errorcurveoptimalspike(isodata_fe, "real", [2, 3])<jupyter_output><empty_output><jupyter_text>Optimal error curves 2If the isotopes to spike are not specified, all possible combinations are shown. This produces Figure 3 in the manuscript.<jupyter_code>ds.errorcurveoptimalspike(isodata_fe, "real")<jupyter_output><empty_output><jupyter_text>Monte Carlo fake mass spec runs`monterun` performs Monte Carlo fake mass spec runs.The example below uses a 50-50 spike-sample mix, the pure 50-50 57Fe-58Fe spike used earlier, a natural fractionation of -0.2, an instrumental fractionation of 1.8, with 1000 Monte Carlo samples. The first 10 mixture measurements are shown.<jupyter_code>measured = ds.monterun(isodata_fe, 0.5, spike, -0.2, 1.8, 1000)
measured[0:10, :]<jupyter_output><empty_output><jupyter_text>Double spike inversions`dsinversion` performs the double spike inversion. Here we run the double spike inversion on the Monte-Carlo generated data with the chosen spike. We then produce a figure showing how the value of alpha varies over the run.<jupyter_code>out = ds.dsinversion(isodata_fe, measured, spike)
plt.plot(out["alpha"])
plt.ylabel(r"$\alpha$");<jupyter_output><empty_output><jupyter_text>Error estimatesThe `errorestimate` routine estimates the errors by linear error propagation. Here we compare the error obtained from the Monte-Carlo simulation with that predicted by linear error propagation. The fact that these are close is a good validation of the linear error propagation method.<jupyter_code>monteerror = np.std(out["alpha"])
predictederror = ds.errorestimate(isodata_fe, 0.5, spike, alpha=-0.2, beta=1.8)[0]
print("Monte Carlo error:", monteerror, "\nPredicted error:", predictederror)<jupyter_output>Monte Carlo error: 0.003707832053847423 
Predicted error: 0.0036364520647345615<jupyter_text>Error model 1The coefficients of the error model are contained in the `errormodel` attribute of the `IsoData` object. The coefficients can be specified for the measured, standard (unspiked run), and double spike compositions. See Appendix C of the manuscript for their definition. The choice of error model is a key part of determining what makes an optimal spike, and the user should carefully consider the assumptions implicit in any error model and decide whether they are appropriate for their system.<jupyter_code>isodata_fe.set_errormodel(R=1e11, deltat=8)
isodata_fe.errormodel["measured"]<jupyter_output><empty_output><jupyter_text>The error model set above fixes the total intensity of the mixture measurement to be 10 V, and the coefficients $a$, $b$, and $c$ correspond to those in (34) in the manuscript. The errors incorporate Johnson–Nyquist noise and counting statistics for a resistance $R = 10^{11}$  Ω and an integration time $\Delta t=8$ s.  Error model 2The function `set_errormodel` can be used to simply change the intensity for the error model. In the example below, a doubling of the intensity decreases the error by roughly a factor of $1/\sqrt{2}$.<jupyter_code>isodata_fe.set_errormodel(intensity=10)
error1 = ds.errorestimate(isodata_fe, 0.5, spike, alpha=-0.2, beta=1.8)[0]
print("Error 1:", error1)
isodata_fe.set_errormodel(intensity=20)
error2 = ds.errorestimate(isodata_fe, 0.5, spike, alpha=-0.2, beta=1.8)[0]
print("Error 2:", error2)
print("Error 2/Error 1:", error2 / error1)
isodata_fe.set_errormodel()  # return error model to defaults<jupyter_output>Error 1: 0.0036364520647345615
Error 2: 0.002457105485500852
Error 2/Error 1: 0.6756875772760133<jupyter_text>Error model 3The default error model of the double spike toolbox fixes the total voltage of the beams for the overall mixture at a given level. When sample-limited it may be more appropriate to consider an error model where the voltage for the sample is fixed (see [John 2012, J. Anal. At. Spectrom.](https://doi.org/10.1039/C2JA30215B)). In the toolbox this can be achieved by changing the error model type from `'fixed-total'` to `'fixed-sample'` The code below recreates Figure 8 of [Klaver and Coath 2019, Geostand. Geoanal. Res.](https://doi.org/10.1111/ggr.12248).<jupyter_code>isodata_ni = ds.IsoData("Ni")
isoinv = [58, 60, 61, 62]
spike6062 = [0.0132, 0.3295, 0.0014, 0.6547, 0.0012]
spike6162 = [0.0109, 0.0081, 0.4496, 0.5297, 0.0017]
# fix the sample intensity at 0.5 V for all runs
isodata_ni.set_errormodel(intensity=0.5, measured_type="fixed-sample")

ds.errorcurve(isodata_ni, spike6062, isoinv, label="$^{60}$Ni-$^{62}$Ni spike")
ds.errorcurve(isodata_ni, spike6162, isoinv, label="$^{61}$Ni-$^{62}$Ni spike")
plt.ylim([0, 0.06])
plt.legend()
isodata_ni.set_errormodel()  # return errormodel to defaults<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09825376601963576, 'fraction_numerical': 0.10469909315746084, 'mean_word_length': 4.3184535671582305, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 79, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '23156', 'license': 'MIT', 'n_tokens_mistral': 5381, 'n_tokens_neox': 4304, 'n_words': 1568, 'path': 'src/example.ipynb', 'repo_name': 'johnrudge/pyspike', 'size': '13343'}