D47calib
Generate, combine, display and apply Δ47 calibrations
This library provides support for:
- computing Δ47 calibrations by applying OGLS regression to sets of (T, Δ47) observations
- combining Δ47 datasets to produce a combined calibration
- various methods useful for creating Δ47 calibration plots
- Using Δ47 calibrations to convert between T and Δ47, keeping track of covariance between inputs
and/or uncertainties/covariance originating from calibration uncertainties. This may be done within
Python code or by using a simple command-line interface (e.g.,
D47calib input.csv > output.csv).
1. Installation
1.1 Library
This should do the trick:
pip install D47calib
Alternatively:
- download or clone the source from https://github.com/mdaeron/D47calib
- chose one of one of the following options:
- copy/move the
/src/D47calibdirectory to somewhere in your Python path - copy/move the
/src/D47calibdirectory to your current working directory - copy/move the
/src/D47calibdirectory to any other location (e.g.,/foo/bar) and include the following code snippet in your scripts:
- copy/move the
import sys
sys.path.append('/foo/bar')
I you don't install from pip, you will probably need to install the requirements listed in pyproject.toml.
1.2 Only install command-line interface using pipx
If you only want to install the CLI, one easy option is to do so using pipx:
pipx install D47calib
Then reopen a shell window and try D47calib --help.
2. Calibrations included
D47calib provides the following pre-built calibrations. All code and raw data used to compute these calibrations can be found in the build_calibs directory.
breitenbach_2018:
Cave pearls analyzed by Breitenbach et al. (2018).
Raw data were obtained from the original study’s supplementary information. The original publication processed data according to two sessions, each 4-5 months long, separated by 2 months. After reprocessing the original raw data using D47crunch, visual inspection of the standardization residuals defined revealed the presence of substantial drifts in both sessions. We thus assigned modified session boundaries defining four continuous measurement periods separated by 21 to 52 days, with new session lengths ranging from 24 to 80 days. The original data was not modified in any other way. Formation temperatures are from table 1 of the original study. We assigned arbitrary 95 % uncertainties of ±1 °C, which seem reasonable for cave environments.
peral_2018:
Planktic foraminifera analyzed by Peral et al. (2018), reprocessed by Daëron & Gray (2023).
Peral et al. [2018] reported Δ47 values of foraminifera from core-tops, both planktic and benthic, whose calcification temperature estimates were recently reassessed by Daëron & Gray (2023). Here we only consider Peral et al.’s planktic data, excluding two benthic samples (cf Daëron & Gray for reasons why we only consider planktic samples for now). In our reprocessing, as in the original study, “samples” are defined by default as a unique combination of core site, species, and size fraction. Δ47 values are then standardized in the usual way, before using D47crunch.combine_samples() to combine all size fractions with the same core and species, except for G. inflata samples (cf Daëron & Gray and accompanying GitHub repository). By properly accounting for analytical error covariance between the Δ47 values to combine, this two-step approach avoids underestimating the final standardization errors.
jautzy_2020:
Synthetic calcites analyzed by Jautzy et al. (2020).
Jautzy et al. reported data from a continuous period spanning 10 months, and used a moving-window approach to standardize their measurements. We assigned sessions defined, whenever possible, as periods of one or more complete weeks enclosing one of more unknown sample analyses. The resulting Δ47 residuals, on the order of 40 ppm (1SD), do not display evidence of instrumental drift. Formation temperatures are from table S2 of the original study. We assigned arbitrary 95 % uncertainties of ±1 °C, which seem reasonable for laboratory experiments.
anderson_2021_mit:
Various natural and synthetic carbonates analyzed at MIT by Anderson et al. (2021).
Raw IRMS data and temperature constraints were obtained from the original study’s supple- mentary information (tables S01 and S02). When reprocesseded the IRMS data we made no changes to the session defintions, but we excluded sessions 5 and 25 because they did not include any unknown sample analyses.
anderson_2021_lsce:
Slow-growing mammillary calcite from Devils Hole and Laghetto Basso analyzed at LSCE by Anderson et al. (2021).
Raw IRMS data is from the original study’s supplementary information (SI-S02). Temperature contraints are from table 1 in Daëron et al. (2019).
fiebig_2021:
Inorganic calcites analyzed by Fiebig et al. (2021).
Temperature contraints are duplicated from the earlier publications where the corresponding samples were first described [Daëron et al., 2019; Jautzy et al., 2020; Anderson et al., 2021]. Raw IRMS data and were obtained from the original study’s supplementary information, and processed as described by Fiebig et al. [2021], jointly using (a) heated and 25 °C-equilibrated CO2 to constrain the scrambling effect and compositional nonlinearity associated with each session, and (b) ETH-1 and ETH-2 reference materials to anchor unknown samples to the I-CDES scale.
huyghe_2022:
Marine calcitic bivalves analyzed by Huyghe et al. (2022).
Huyghe et al. reported Δ47 values of modern calcitic bivalves collected from localities with good environmental constraints. As was done in the original publication, different bivalve individuals were initially treated as distinct analytical samples. In some sites with strong seasonality, individuals were sub-sampled into winter-calcified a summer-calcified fractions. Δ47 values were then standardized in the usual way, before using D47crunch.combine_samples() method to combine all samples from the same locality. Calcification temperature estimates are from the original study.
devils_laghetto_2023:
Combined data set of slow-growing mammillary calcite from Devils Hole and Laghetto Basso, analyzed both at LSCE by Anderson et al. (2021) and at GU by Fiebig et al. (2021).
OGLS23:
Combined data set including all of the above. For a detailed discussion of goodness-of-fit and regression uncertainties, see Daëron & Vermeesch (2024). Also aliased as ogls_2023 for backward-compatibility.
3. Command-line interface
D47calib also provides a command-line interface (CLI) for converting between Δ47 and temperature values, computing uncertainties for each computed value (and how these uncertainties are correlated with each other) from different sources (from calibration errors alone, from measurement errors alone, and from both). The computed uncertainties are provided as standard errors, correlation matrix and/or covariance matrix. Input and output files may be comma-separated, tab-separated, or printed out as visually aligned data columns.
3.1 Simple examples
Start with the simplest input file possible (here named input.csv):
D47
0.567
Then process it:
D47calib input.csv
This prints out:
D47 T T_SE_from_calib T_correl_from_calib T_SE_from_input T_correl_from_input T_SE_from_both T_correl_from_both
0.567 34.20 0.38 1.000 0.00 1.000 0.38 1.000
Tis the temperature corresponding to aD47value of 0.567 ‰ according to the default calibration (OGLS23).T_SE_from_calibis the standard error onTfrom the calibration uncertaintyT_correl_from_calibis the correlation matrix for theT_SE_from_calibvalues. Because here there is only one value, this is a 1-by-1 matrix with a single value of one, which is not very exciting.T_SE_from_inputis the standard error onTfrom the measurement uncertainties onD47. Because these are not specified here,T_SE_from_inputis equal to zero.T_correl_from_inputis, predictably, the correlation matrix for theT_SE_from_inputvalues. Because here there is only one value, this is a 1-by-1 matrix with a single value of one, you know the drill.T_SE_from_bothis the standard error onTobtained by combining the two previously considered sources of uncertainties.T_correl_from_bothis what you expect it to be if you've been reading so far. Can you guess why it is a 1-by-1 matrix with a single value of one?
3.1.1 Adding D47 measurement uncertainties
This can be done by adding a column to input.csv:
D47,D47_SE
0.567,0.008
Because this is not very human-friendly, we'll replace the comma separators by whitespace. We'll also add a column listing sample names:
Sample D47 D47_SE
FOO-1 0.567 0.008
Then process it. We're adding an option (-i ' ', or --delimiter-in ' ') specifying that we're no longer using commas but whitespaces as delimiters:
D47calib -i ' ' input.csv
This yields:
Sample D47 D47_SE T T_SE_from_calib T_correl_from_calib T_SE_from_input T_correl_from_input T_SE_from_both T_correl_from_both
FOO-1 0.567 0.008 34.20 0.38 1.000 2.91 1.000 2.94 1.000
You can see that T_SE_from_input is now much larger than T_SE_from_calib, and that the combined T_SE_from_both is equal to the quadratic sum of T_SE_from_calib and T_SE_from_input.
3.1.2 Converting more than one measurement
Let's add lines to our input file:
Sample D47 D47_SE
FOO-1 0.567 0.008
BAR-2 0.575 0.009
BAZ-3 0.582 0.007
Which yields:
Sample D47 D47_SE T T_SE_from_calib T_correl_from_calib T_SE_from_input T_correl_from_input T_SE_from_both T_correl_from_both
FOO-1 0.567 0.008 34.20 0.38 1.000 0.996 0.987 2.91 1.000 0.000 0.000 2.94 1.000 0.015 0.019
BAR-2 0.575 0.009 31.33 0.37 0.996 1.000 0.997 3.18 0.000 1.000 0.000 3.21 0.015 1.000 0.017
BAZ-3 0.582 0.007 28.89 0.36 0.987 0.997 1.000 2.42 0.000 0.000 1.000 2.44 0.019 0.017 1.000
A notable change are the 3-by-3 correlation matrices, which tell us how the T errors or these three measurements covary. T_correl_from_calib shows that the T_SE_from_calib errors are strongly correlated, because the three D47 values are close to each other. T_correl_from_input indicates statistically independent T_SE_from_input errors. This may be true or not, but it is the expected result because our input file does not include any information on how the D47_SE errors may covary (see below how this additional info may be specified). Thus in this case D47calib assumes that the D47 values are statistically independent (gentle reminder: this is often not the case, see below).
Note that because T_SE_from_input errors are much larger than T_SE_from_calib errors, the combined T_SE_from_both errors are only weakly correlated, as seen in T_correl_from_both.
3.1.3 Accounting for correlations in D47 errors
Because Δ47 measurements performed in the same analytical session(s) are not statistically independent, we may add to input.csv a correlation matrix describing how D47_SE errors covary.
One simple way to compute this correlation matrix is to use the save_D47_correl() method from the D47crunch library (PyPI, GitHub, Zenodo) described by Daëron (2021).
Sample D47 D47_SE D47_correl
FOO-1 0.567 0.008 1.00 0.25 0.25
BAR-2 0.575 0.009 0.25 1.00 0.25
BAZ-3 0.582 0.007 0.25 0.25 1.00
This yields:
Sample D47 D47_SE D47_correl T T_SE_from_calib T_correl_from_calib T_SE_from_input T_correl_from_input T_SE_from_both T_correl_from_both
FOO-1 0.567 0.008 1.00 0.25 0.25 34.20 0.38 1.000 0.996 0.987 2.91 1.000 0.250 0.250 2.94 1.000 0.261 0.264
BAR-2 0.575 0.009 0.25 1.00 0.25 31.33 0.37 0.996 1.000 0.997 3.18 0.250 1.000 0.250 3.21 0.261 1.000 0.263
BAZ-3 0.582 0.007 0.25 0.25 1.00 28.89 0.36 0.987 0.997 1.000 2.42 0.250 0.250 1.000 2.44 0.264 0.263 1.000
What changed ? We now have propagated D47_correl into T_correl_from_input, and this is accounted for in the combined correlation matrix T_correl_from_both. Within the framework of our initial assumptions (multivariate Gaussian errors, first-order linear propagation of uncertainties...), this constitutes the “best” (or rather, the most “information-complete”) description of uncertainties constraining our final T estimates.
With increasing number of measurements, these correlation matrices become quite large, so that it becomes useless to print them out visually. To facilitate using the output of D47calib as an input to another piece of software, one may use the -j or --delimiter-out option to use machine-readable delimiters such as commas or tabs, and the '-o' or --output-file option to save the output as a file instead of printing it out:
D47calib -i ' ' -j ',' -o output.csv input.csv
This will create the following output.csv file:
Sample,D47,D47_SE,D47_correl,,,T,T_SE_from_calib,T_correl_from_calib,,,T_SE_from_input,T_correl_from_input,,,T_SE_from_both,T_correl_from_both,,
FOO-1,0.567,0.008,1.00,0.25,0.25,34.20,0.38,1.000,0.996,0.987,2.91,1.000,0.250,0.250,2.94,1.000,0.261,0.264
BAR-2,0.575,0.009,0.25,1.00,0.25,31.33,0.37,0.996,1.000,0.997,3.18,0.250,1.000,0.250,3.21,0.261,1.000,0.263
BAZ-3,0.582,0.007,0.25,0.25,1.00,28.89,0.36,0.987,0.997,1.000,2.42,0.250,0.250,1.000,2.44,0.264,0.263,1.000
Hint for Mac users: Quick Look (or “spacebar preview”, i.e. what happens when you select a file in the Finder and press the spacebar once) provides you with a nice view of a csv file when you just want to check the results visually, as long as you use a comma delimiter.
3.1.4 Converting from T to D47
Everything described above works in the other direction as well, without changing anything to the command-line instruction:
T T_SE
0 0.5
10 1.0
20 2.0
Yields:
T T_SE D47 D47_SE_from_calib D47_correl_from_calib D47_SE_from_input D47_correl_from_input D47_SE_from_both D47_correl_from_both
0 0.5 0.6798 0.0016 1.000 0.969 0.848 0.0020 1.000 0.000 0.000 0.0025 1.000 0.210 0.091
10 1.0 0.6424 0.0013 0.969 1.000 0.952 0.0035 0.000 1.000 0.000 0.0038 0.210 1.000 0.056
20 2.0 0.6090 0.0011 0.848 0.952 1.000 0.0063 0.000 0.000 1.000 0.0064 0.091 0.056 1.000
3.2 Integration with D47crunch
Starting with the following input file rawdata.csv:
UID Sample Session d45 d46 d47 d48 d49
1 BAZ-3 Session_01 -5.394955 7.938127 1.887702 15.671857 9.739724
2 ETH-1 Session_01 6.050787 10.800497 16.232192 21.429453 27.780042
3 ETH-3 Session_01 5.726647 11.144602 16.634507 22.275401 28.306614
4 ETH-1 Session_01 6.009875 10.711152 16.087994 21.275325 27.780042
5 FOO-1 Session_01 -0.848380 2.872996 1.535960 5.431873 4.665655
6 ETH-2 Session_01 -5.973121 -5.894178 -12.599303 -12.037594 -18.023381
7 BAR-2 Session_01 -9.941731 10.985508 0.208381 21.832546 10.707292
8 ETH-3 Session_01 5.755765 11.179326 16.689870 22.294132 28.306614
9 ETH-2 Session_01 -5.981599 -6.011356 -12.728742 -12.335559 -18.023381
10 ETH-2 Session_01 -5.991066 -5.968789 -12.685206 -12.219412 -18.023381
11 ETH-3 Session_01 5.734551 11.043266 16.556578 22.005474 28.306614
12 ETH-1 Session_01 5.994522 10.810755 16.181039 21.445703 27.780042
13 ETH-2 Session_02 -5.993596 -6.086738 -12.816033 -12.452397 -18.023381
14 ETH-3 Session_02 5.720443 11.187928 16.689130 22.308575 28.306614
15 ETH-1 Session_02 6.019913 10.773652 16.165645 21.309536 27.780042
16 ETH-1 Session_02 5.995178 10.702934 16.073811 21.169324 27.780042
17 BAR-2 Session_02 -9.951732 10.964153 0.171443 21.779172 10.707292
18 FOO-1 Session_02 -0.835316 2.868058 1.539147 5.483628 4.665655
19 ETH-2 Session_02 -5.952660 -5.887239 -12.582248 -12.081588 -18.023381
20 ETH-2 Session_02 -5.983050 -5.953794 -12.679809 -12.224144 -18.023381
21 ETH-3 Session_02 5.717666 11.175023 16.654427 22.251497 28.306614
22 ETH-3 Session_02 5.756394 11.109906 16.635806 22.189179 28.306614
23 ETH-1 Session_02 6.029949 10.682183 16.080186 21.182471 27.780042
24 BAZ-3 Session_02 -5.337619 7.818571 1.818572 15.400966 9.739724
The following script will read thart raw data, fully process it, convert the standardized output to temperatures, and save the final results to a file named output.csv:
D47crunch rawdata.csv
D47calib -o output.csv -j '>' output/D47_correl.csv
With the contents of output.csv being:
Sample D47 D47_SE D47_correl T T_SE_from_calib T_correl_from_calib T_SE_from_input T_correl_from_input T_SE_from_both T_correl_from_both
BAR-2 0.6777 0.0066 1.0000 0.3586 0.2798 0.51 0.41 1.000 0.731 -0.036 1.68 1.000 0.359 0.280 1.73 1.000 0.373 0.262
BAZ-3 0.5894 0.0060 0.3586 1.0000 0.2473 26.34 0.35 0.731 1.000 0.654 2.02 0.359 1.000 0.247 2.05 0.373 1.000 0.263
FOO-1 0.4873 0.0056 0.2798 0.2473 1.0000 68.21 0.69 -0.036 0.654 1.000 2.82 0.280 0.247 1.000 2.90 0.262 0.263 1.000
If a simpler output is required, just add the --ignore-correl or -g option to the second line above, which should yield:
Sample D47 D47_SE T T_SE_from_calib T_SE_from_input T_SE_from_both
BAR-2 0.6777 0.0066 0.51 0.41 1.68 1.73
BAZ-3 0.5894 0.0060 26.34 0.35 2.02 2.05
FOO-1 0.4873 0.0056 68.21 0.69 2.82 2.90
3.3 Further customizing the CLI
A complete list of options is provided by D47calib --help.
3.3.1 Using covariance instead of correlation matrix as input
Just provide D47_covar (or T_covar when converting in the other direction) in the input file instead of D47_SE and D47_correl.
3.3.2 Reporting covariance instead of correlation matrices in the output
Use the --return-covar option.
3.3.3 Reporting neither covariances nor correlations in the output
If you don't care about all this covariance nonsense, or just wish for an output that does't hurt your eyes, you can use the --ignore-correl option. Standard errors will still be reported.
3.3.4 Excluding or only including certain lines (samples) from the input
To filter the samples (lines) to process using --exclude-samples and --include-samples, first add a Sample column to the input data, assign a sample name to each line.
Then to exclude some samples, provide the --exclude-samples option with the name of a file where each line is one sample to exclude.
To exclude all samples except those listed in a file, provide the --include-samples option with the name of that file, where each line is one sample to include.
3.3.5 Changing the numerical precision of the output
This is controlled by the following options:
--T-precisionor-p(default: 2): AllTandT_SE_*values--D47-precisionor-q(default: 4): AllD47andD47_SE_*values--correl-precisionor-r(default: 3): All*_correl_*values--covar-precisionor-s(default: 3): All*_covar_*values
3.3.6 Using a different Δ47 calibration
You may use a different calibration than the default OGLS23 using the --calib or -c option. Any predefeined calibration from the D47calib library is a valid option.
You may also specify an arbitrary polynomial function of inverse T, by creating a file (e.g., calib.csv) with the following format:
degree coef covar
0 0.1741 2.4395e-05 -0.0262821 5.634934
1 -17.889 -0.0262821 32.17712 -7223.86
2 42614 5.634934 -7223.86 1654633.996
Then using -c calib.csv will use this new calibration.
If you don't know/care about the covariance of the calibration coefficients, just leave out the covar terms:
degree coef
0 0.1741
1 -17.889
2 42614
In this case, all *_SE_from_calib outputs will be equal to zero, but the *_from_input uncertainties will still be valid (and identical to *_from_both uncertainties, since we are ignoring calibration uncertainties).
1""" 2Generate, combine, display and apply Δ47 calibrations 3 4This library provides support for: 5 6- computing Δ47 calibrations by applying OGLS regression to sets of (T, Δ47) observations 7- combining Δ47 datasets to produce a combined calibration 8- various methods useful for creating Δ47 calibration plots 9- Using Δ47 calibrations to convert between T and Δ47, keeping track of covariance between inputs 10and/or uncertainties/covariance originating from calibration uncertainties. This may be done within 11Python code or by using a simple command-line interface (e.g., `D47calib input.csv > output.csv`). 12 13.. include:: ../../docpages/install.md 14.. include:: ../../docpages/calibs.md 15.. include:: ../../docpages/cli.md 16 17* * * 18""" 19 20__author__ = 'Mathieu Daëron' 21__contact__ = 'daeron@lsce.ipsl.fr' 22__copyright__ = 'Copyright (c) 2025 Mathieu Daëron' 23__license__ = 'MIT License - https://opensource.org/licenses/MIT' 24__date__ = '2025-12-15' 25__version__ = '1.4.0' 26 27 28import typer 29import typer.rich_utils 30import sys 31from typing_extensions import Annotated 32import ogls as _ogls 33import numpy as _np 34from scipy.linalg import block_diag as _block_diag 35from scipy.interpolate import interp1d as _interp1d 36from matplotlib import pyplot as _ppl 37 38typer.rich_utils.STYLE_HELPTEXT = '' 39 40class D47calib(_ogls.InverseTPolynomial): 41 """ 42 Δ47 calibration class based on OGLS regression 43 of Δ47 as a polynomial function of inverse T. 44 """ 45 46 def __init__(self, 47 samples, T, D47, 48 sT = None, 49 sD47 = None, 50 degrees = [0,2], 51 xpower = 2, 52 name = '', 53 label = '', 54 description = '', 55 **kwargs, 56 ): 57 """ 58 ### Parameters 59 60 + **samples**: a list of N sample names. 61 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 62 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 63 + **sT**: uncertainties on `T`. If specified as: 64 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 65 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 66 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 67 + **sD47**: uncertainties on `D47`. If specified as: 68 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 69 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 70 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 71 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 72 + **name**: a human-readable, short name assigned to the calibration. 73 + **label**: a short description of the calibration, e.g., to be used in legends. 74 + **description**: a longer description, including relevant references/DOIs. 75 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 76 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 77 78 ### Notable attributes 79 80 + **N**: 81 The total number of observations (samples) in the calibration data. 82 + **samples**: 83 The list sample names. 84 + **T**: 85 1-D `ndarray` of temperatures in degrees C. 86 + **D47**: 87 1-D `ndarray` of Δ47 values in permil. 88 + **sT**: 89 2-D `ndarray` equal to the full (co)variance matrix for `T`. 90 + **D47**: 91 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 92 + **xpower**: 93 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 94 It is possible to change this behavior to use a different power of 1/T. 95 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 96 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 97 + **bfp**: 98 The best-fit parameters of the regression. 99 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 100 + **bff()**: 101 The best-fit polynomial function of inverse T, defined as: 102 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 103 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 104 105 106 ### Examples 107 108 A very simple example: 109 110 ````py 111 .. include:: ../../code_examples/D47calib_init/example.py 112 ```` 113 114 Should yield: 115 116 ```` 117 .. include:: ../../code_examples/D47calib_init/output.txt 118 ```` 119 120 """ 121 122 self.samples = samples[:] 123 self.name = name 124 self.label = label 125 self.description = description 126 self.D47 = _np.asarray(D47, dtype = 'float') 127 self.N = self.D47.size 128 129 if sD47 is None: 130 self.sD47 = _np.zeros((self.N, self.N)) 131 else: 132 self.sD47 = _np.asarray(sD47) 133 if len(self.sD47.shape) == 1: 134 self.sD47 = _np.diag(self.sD47**2) 135 elif len(self.sD47.shape) == 0: 136 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 137 138 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 139 140 if self.bfp is None: 141 self.regress() 142 143 self._bff_deriv = lambda x: _np.array([k * self.bfp[f'a{k}'] * x**(k-1) for k in degrees if k > 0]).sum(axis = 0) 144 145 xi = _np.linspace(0,200**-1,1001) 146 self._inv_bff = _interp1d(self.bff(xi), xi) 147 148 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 149 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 150 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 151 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1 152 153 def __repr__(self): 154 return f'<D47calib: {self.name}>' 155 156 def invT_xaxis(self, 157 xlabel = None, 158 Ti = [0,20,50,100,250,1000], 159 ): 160 """ 161 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 162 but labeled in degrees Celsius. 163 164 ### Parameters 165 166 + **xlabel**: 167 Custom label for X axis (`r'$1\\,/\\,T^2$'` by default) 168 + **Ti**: 169 Specify tick locations for X axis, in degrees C. 170 171 ### Returns 172 173 + an `matplotlib.axes.Axes` instance 174 175 ### Examples 176 177 ````py 178 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 179 ```` 180 181 This should result in something like this: 182 183 <img align="center" src="example_invT_xaxis_1.png"> 184 185 It is also possible to define the X axis using a different power of 1/T 186 by first redefining the `xpower` attribute: 187 188 ````py 189 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 190 ```` 191 192 This should result in something like this: 193 194 <img align="center" src="example_invT_xaxis_2.png"> 195 """ 196 if xlabel is None: 197 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 198 _ppl.xlabel(xlabel) 199 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 200 ax = _ppl.gca() 201 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 202 ax.tick_params(which="major") 203 204 return ax 205 206 207 def plot_data(self, label = False, **kwargs): 208 """ 209 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 210 211 ### Parameters 212 213 + **label**: 214 + If `label` is a string, use this string as `label` for the underlyig 215 `matplotlib.pyplot.plot()` call. 216 + If `label = True`, use the caller's `label` attribute instead. 217 + If `label = False`, no label is specified (default behavior). 218 + **kwargs**: 219 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 220 221 ### Returns 222 223 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 224 225 ### Example 226 227 ````py 228 from matplotlib import pyplot as ppl 229 from D47calib import huyghe_2022 as calib 230 231 fig = ppl.figure(figsize = (5,3)) 232 ppl.subplots_adjust(bottom = .25, left = .15) 233 calib.invT_xaxis(Ti = [0,10,25]) 234 calib.plot_data(label = True) 235 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 236 ppl.legend() 237 ppl.savefig('example_plot_data.png', dpi = 100) 238 ````` 239 240 This should result in something like this: 241 242 <img align="center" src="example_plot_data.png"> 243 """ 244# if 'mec' not in kwargs: 245# kwargs['mec'] = self.color 246 if label is not False: 247 kwargs['label'] = self.label if label is True else label 248 return _ogls.InverseTPolynomial.plot_data(self, **kwargs) 249 250 251 def plot_error_bars(self, **kwargs): 252 """ 253 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 254 255 ### Parameters 256 257 + **kwargs**: 258 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 259 260 ### Returns 261 262 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 263 264 ### Example 265 266 ````py 267 from matplotlib import pyplot as ppl 268 from D47calib import huyghe_2022 as calib 269 270 fig = ppl.figure(figsize = (5,3)) 271 ppl.subplots_adjust(bottom = .25, left = .15) 272 calib.invT_xaxis(Ti = [0,10,25]) 273 calib.plot_error_bars(alpha = .4) 274 calib.plot_data(label = True) 275 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 276 ppl.legend() 277 ppl.savefig('example_plot_error_bars.png', dpi = 100) 278 ````` 279 280 This should result in something like this: 281 282 <img align="center" src="example_plot_error_bars.png"> 283 """ 284# if 'ecolor' not in kwargs: 285# kwargs['ecolor'] = self.color 286 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs) 287 288 289 def plot_error_ellipses(self, **kwargs): 290 """ 291 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 292 293 ### Parameters 294 295 + **kwargs**: 296 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 297 298 ### Returns 299 300 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 301 302 ### Example 303 304 ````py 305 from matplotlib import pyplot as ppl 306 from D47calib import huyghe_2022 as calib 307 308 fig = ppl.figure(figsize = (5,3)) 309 ppl.subplots_adjust(bottom = .25, left = .15) 310 calib.invT_xaxis(Ti = [0,10,25]) 311 calib.plot_error_ellipses(alpha = .4) 312 calib.plot_data(label = True) 313 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 314 ppl.legend() 315 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 316 ````` 317 318 This should result in something like this: 319 320 <img align="center" src="example_plot_error_ellipses.png"> 321 """ 322# if 'ec' not in kwargs: 323# kwargs['ec'] = self.color 324 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs) 325 326 327 def plot_bff(self, label = False, **kwargs): 328 """ 329 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 330 331 ### Parameters 332 333 + **label**: 334 + If `label` is a string, use this string as `label` for the underlyig 335 `matplotlib.pyplot.plot()` call. 336 + If `label = True`, use the caller's `label` attribute instead. 337 + If `label = False`, no label is specified (default behavior). 338 + **kwargs**: 339 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 340 341 ### Returns 342 343 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 344 345 ### Example 346 347 ````py 348 from matplotlib import pyplot as ppl 349 from D47calib import huyghe_2022 as calib 350 351 fig = ppl.figure(figsize = (5,3)) 352 ppl.subplots_adjust(bottom = .25, left = .15) 353 calib.invT_xaxis(Ti = [0,10,25]) 354 calib.plot_bff(label = True, dashes = (8,2,2,2)) 355 calib.plot_data() 356 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 357 ppl.legend() 358 ppl.savefig('example_plot_bff.png', dpi = 100) 359 ````` 360 361 This should result in something like this: 362 363 <img align="center" src="example_plot_bff.png"> 364 """ 365# if 'color' not in kwargs: 366# kwargs['color'] = self.color 367 if label is not False: 368 kwargs['label'] = self.label if label is True else label 369 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs) 370 371 372 def plot_bff_ci(self, **kwargs): 373 """ 374 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 375 376 ### Parameters 377 378 + **label**: 379 + **kwargs**: 380 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 381 382 ### Returns 383 384 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 385 386 ### Example 387 388 ````py 389 from matplotlib import pyplot as ppl 390 from D47calib import huyghe_2022 as calib 391 392 fig = ppl.figure(figsize = (5,3)) 393 ppl.subplots_adjust(bottom = .25, left = .15) 394 calib.invT_xaxis(Ti = [0,10,25]) 395 calib.plot_bff_ci(alpha = .15) 396 calib.plot_bff(label = True, dashes = (8,2,2,2)) 397 calib.plot_data() 398 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 399 ppl.legend() 400 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 401 ````` 402 403 This should result in something like this: 404 405 <img align="center" src="example_plot_bff_ci.png"> 406 """ 407# if 'color' not in kwargs: 408# kwargs['color'] = self.color 409 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs) 410 411 def T47(self, 412 D47 = None, 413 sD47 = None, 414 T=None, 415 sT = None, 416 error_from = 'both', 417 return_covar = False, 418 ): 419 ''' 420 When `D47` is input, computes corresponding T value(s). 421 `D47` input may be specified as a scalar, or as a 1-D array. 422 `T` output will then have the same type and size as `D47`. 423 424 When `T` is input, computes corresponding Δ47 value(s). 425 `T` input may be specified as a scalar, or as a 1-D array. 426 `D47` output will then have the same type and size as `T`. 427 428 Only one of either `D47` or `T` may be specified as input. 429 430 **Arguments:** 431 432 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 433 * `sD47`: Δ47 uncertainties, which may be: 434 - `None` (default) 435 - `float` or `int` (uniform standard error on `D47`) 436 - 1-D array (standard errors on `D47`) 437 - 2-D array (covariance matrix for `D47`) 438 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 439 * `sT`: T uncertainties, which may be: 440 - `None` (default) 441 - `float` or `int` (uniform standard error on `T`) 442 - 1-D array (standard errors on `T`) 443 - 2-D array (variance-covariance matrix for `T`) 444 * `error_from`: if set to `'both'` (default), returned errors take into account 445 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 446 if set to `'calib'`, only calibration uncertainties are accounted for; 447 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 448 * `return_covar`: (False by default) whether to return the full covariance matrix 449 for returned `T` or `D47` values, otherwise return standard errors for the returned 450 `T` or `D47` values instead. 451 452 **Returns (with `D47` input):** 453 454 * `T`: temperature value(s) computed from `D47` 455 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 456 457 **Returns (with `T` input):** 458 459 * `D47`: Δ47 value(s) computed from `D47` 460 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 461 462 ### Example 463 464 ````py 465 import numpy as np 466 from matplotlib import pyplot as ppl 467 from D47calib import OGLS23 as calib 468 469 X = np.linspace(1473**-2, 270**-2) 470 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 471 472 fig = ppl.figure(figsize = (5,3)) 473 ppl.subplots_adjust(bottom = .25, left = .15) 474 calib.invT_xaxis() 475 ppl.plot(X, 1000 * sD47, 'r-') 476 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 477 ppl.savefig('example_SE47.png', dpi = 100) 478 ````` 479 480 This should result in something like this: 481 482 <img src="example_SE47.png"> 483 ''' 484 485 if D47 is None and T is None: 486 raise ValueError('Either D47 or T must be specified, but both are undefined.') 487 488 if D47 is not None and T is not None: 489 raise ValueError('Either D47 or T must be specified, but not both.') 490 491 if T is not None: 492 493 D47 = self._D47_from_T(T) 494 Np = len(self.degrees) 495 N = D47.size 496 497 ### Compute covariance matrix of (*bfp, *T): 498 CM = _np.zeros((Np+N, Np+N)) 499 500 if error_from in ['calib', 'both']: 501 CM[:Np, :Np] = self.bfp_CM[:,:] 502 503 if (sT is not None) and error_from in ['sT', 'both']: 504 _sT = _np.asarray(sT) 505 if _sT.ndim == 0: 506 for k in range(N): 507 CM[Np+k, Np+k] = _sT**2 508 elif _sT.ndim == 1: 509 for k in range(N): 510 CM[Np+k, Np+k] = _sT[k]**2 511 elif _sT.ndim == 2: 512 CM[-N:, -N:] = _sT[:,:] 513 514 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 515 _T = _np.asarray(T) 516 if _T.ndim == 0: 517 _T = _np.expand_dims(_T, 0) 518 J = _np.zeros((N, Np+N)) 519 520 if (sT is not None) and error_from in ['sT', 'both']: 521 for k in range(N): 522 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 523 524 if error_from in ['calib', 'both']: 525 526 for k in range(Np): 527 528 p1 = {_: self.bfp[_] for _ in self.bfp} 529 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 530 531 p2 = {_: self.bfp[_] for _ in self.bfp} 532 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 533 534 J[:, k] = (self.model_fun(p1, (_T+273.15)**-1) - self.model_fun(p2, (_T+273.15)**-1)) / (0.002 * self.bfp_CM[k,k]**.5) 535 536 ### Error propagation: 537 CM_D47 = J @ CM @ J.T 538 539 if return_covar: 540 return D47, CM_D47 541 else: 542 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 543 544 if D47 is not None: 545 546 T = self._T_from_D47(D47) 547 Np = len(self.degrees) 548 N = T.size 549 550 ### Compute covariance matrix of (*bfp, *T): 551 CM = _np.zeros((Np+N, Np+N)) 552 553 if error_from in ['calib', 'both']: 554 CM[:Np, :Np] = self.bfp_CM[:,:] 555 556 if (sD47 is not None) and error_from in ['sD47', 'both']: 557 _sD47 = _np.asarray(sD47) 558 if _sD47.ndim == 0: 559 for k in range(N): 560 CM[Np+k, Np+k] = _sD47**2 561 elif _sD47.ndim == 1: 562 for k in range(N): 563 CM[Np+k, Np+k] = _sD47[k]**2 564 elif _sD47.ndim == 2: 565 CM[-N:, -N:] = _sD47[:,:] 566 567 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 568 _D47 = _np.asarray(D47) 569 if _D47.ndim == 0: 570 _D47 = _np.expand_dims(_D47, 0) 571 J = _np.zeros((N, Np+N)) 572 if (sD47 is not None) and error_from in ['sD47', 'both']: 573 for k in range(N): 574 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 575 if error_from in ['calib', 'both']: 576 577 xi = _np.linspace(0,200**-1,1001)[1:] 578 for k in range(Np): 579 580 if self.bfp_CM[k,k]: 581 _epsilon_ = self.bfp_CM[k,k]**.5 582 else: 583 _epsilon_ = 1e-6 584 585 p1 = {_: self.bfp[_] for _ in self.bfp} 586 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 587 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 588 589 p2 = {_: self.bfp[_] for _ in self.bfp} 590 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 591 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 592 593 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 594 595 ### Error propagation: 596 CM_T = J @ CM @ J.T 597 598 if return_covar: 599 return T, CM_T 600 else: 601 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5 602 603 604 def plot_T47_errors( 605 self, 606 calibname = None, 607 rD47 = 0.010, 608 Nr = [2,4,8,12,20], 609 Tmin = 0, 610 Tmax = 120, 611 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 612 yscale = 'lin', 613 ): 614 """ 615 Plot SE of T reconstructed using the calibration as a function of T for various 616 combinations of analytical precision and number of analytical replicates. 617 618 **Arguments** 619 620 + **calibname**: 621 Which calibration name to display. By default, use `label` attribute. 622 + **rD47**: 623 Analytical precision of a single analysis. 624 + **Nr**: 625 A list of lines to plot, each corresponding to a given number of replicates. 626 + **Tmin**: 627 Minimum T to plot. 628 + **Tmax**: 629 Maximum T to plot. 630 + **colors**: 631 A list of colors to distinguish the plotted lines. 632 + **yscale**: 633 + If `'lin'`, the Y axis uses a linear scale. 634 + If `'log'`, the Y axis uses a logarithmic scale. 635 636 **Example** 637 638 ````py 639 from matplotlib import pyplot as ppl 640 from D47calib import devils_laghetto_2023 as calib 641 642 fig = ppl.figure(figsize = (3.5,4)) 643 ppl.subplots_adjust(bottom = .2, left = .15) 644 calib.plot_T47_errors( 645 calibname = 'Devils Laghetto calibration', 646 Nr = [1,2,4,16], 647 Tmin =0, 648 Tmax = 40, 649 ) 650 ppl.savefig('example_SE_T.png', dpi = 100) 651 ```` 652 653 This should result in something like this: 654 655 <img src="example_SE_T.png"> 656 """ 657 658 if calibname is None: 659 calibname = self.label 660 661 Nr = _np.array(Nr) 662 if len(colors) < Nr.size: 663 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 664 from colorsys import hsv_to_rgb 665 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 666 colors = [hsv_to_rgb(*x) for x in hsv] 667 668 Ti = _np.linspace(Tmin, Tmax) 669 D47i, _ = self.T47(T = Ti) 670 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 671 672 ymax, ymin = 0, 1e6 673 for N,c in zip(Nr, colors): 674 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 675 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 676 ymin = min(ymin, min(sT)) 677 ymax = max(ymax, max(sT)) 678 679 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 680 681 _ppl.legend(fontsize=9) 682 _ppl.xlabel("T (°C)") 683 684 _ppl.ylabel("Standard error on reconstructed T (°C)") 685 686 # yticks([0,.5,1,1.5,2]) 687 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 688 _ppl.grid( alpha = .25) 689 if yscale == 'lin': 690 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 691 t1, t2 = self.T.min(), self.T.max() 692 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 693 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 694 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 695 elif yscale == 'log': 696 ymin /= 2 697 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 698 _ppl.yscale('log') 699 t1, t2 = self.T.min(), self.T.max() 700 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 701 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 702 703 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 704 """ 705 Write calibration data to a csv file. 706 707 ### Parameters 708 709 + **csvfile**: 710 The filename to write data to. 711 + **sep**: 712 The separator between CSV fields. 713 + **label**: 714 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 715 + If specified as a `str`, include a `Dataset` column with that string. 716 + If specified as `False`, do not include a `Dataset` column. 717 + **T_correl**: 718 + If `True`, include correlations between all `T` values. 719 + **D47_correl**: 720 + If `True`, include correlations between all `D47` values. 721 722 ### Example 723 724 ````py 725 D47calib.huyghe_2022.export_data( 726 csvfile = 'example_export_data.csv', 727 T_correl = True, 728 D47_correl = True, 729 ) 730 ```` 731 732 This should result in something like this ([link](example_export_data.csv)): 733 734 .. include:: ../../docs/example_export_data.md 735 736 """ 737 n = len(str(self.N)) 738 739 with open(csvfile, 'w') as f: 740 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 741 742 if label: 743 f.write(f'{sep}Dataset') 744 745 if T_correl: 746 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 747 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 748 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 749 750 if D47_correl: 751 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 752 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 753 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 754 755 for k, (s, T, sT, D47, sD47) in enumerate(zip( 756 self.samples, 757 self.T, 758 _np.diag(self.sT)**.5, 759 self.D47, 760 _np.diag(self.sD47)**.5, 761 )): 762 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 763 if label: 764 if label is True: 765 f.write(f'{sep}{self.label}') 766 else: 767 f.write(f'{sep}{label}') 768 if T_correl: 769 f.write(sep.join(['']+[ 770 f'{Tcorrel[k,_]:.0f}' 771 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 772 else f'{Tcorrel[k,_]:.6f}' 773 for _ in range(self.N)])) 774 if D47_correl: 775 f.write(sep.join(['']+[ 776 f'{D47correl[k,_]:.0f}' 777 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 778 else f'{D47correl[k,_]:.6f}' 779 for _ in range(self.N)])) 780 781 782 def export(self, name, filename): 783 """ 784 Save `D47calib` object as an importable file. 785 786 ### Parameters 787 788 + **name**: 789 The name of the variable to export. 790 + **filename**: 791 The filename to write to. 792 793 ### Example 794 795 ````py 796 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 797 ```` 798 799 This should result in a `bar.py` file with the following contents: 800 801 ````py 802 foo = D47calib( 803 samples = ['LGB-2', 'DVH-2'], 804 T = [7.9, 33.7], 805 D47 = [0.6485720997671647, 0.5695972909966959], 806 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 807 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 808 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 809 label = 'Slow-growing calcites from Anderson et al. (2021)', 810 color = (0, 0.5, 0), 811 degrees = [0, 2], 812 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 813 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 814 chisq = 6.421311854486162e-27, 815 Nf = 0, 816 ) 817 ```` 818 """ 819 with open(filename, 'w') as f: 820 f.write(f''' 821{name} = D47calib( 822 samples = {self.samples}, 823 T = {self.T.tolist()}, 824 D47 = {self.D47.tolist()}, 825 sT = {self.sT.tolist()}, 826 sD47 = {self.sD47.tolist()}, 827 degrees = {self.degrees}, 828 description = {repr(self.description)}, 829 name = {repr(self.name)}, 830 label = {repr(self.label)}, 831 bfp = {({k: float(self.bfp[k]) for k in self.bfp})}, 832 bfp_CM = {self.bfp_CM.tolist()}, 833 chisq = {self.chisq}, 834 cholesky_residuals = {self.cholesky_residuals.tolist()}, 835 aic = {self.aic}, 836 bic = {self.bic}, 837 ks_pvalue = {self.ks_pvalue}, 838 ) 839''') 840 841def combine_D47calibs(calibs, degrees = [0,2], same_T = [], exclude_samples = []): 842 ''' 843 Combine data from several `D47calib` instances. 844 845 ### Parameters 846 847 + **calibs**: 848 A list of `D47calib` instances 849 + **degrees**: 850 The polynomial degrees of the combined regression. 851 + **same_T**: 852 Use this `list` to specify when samples from different calibrations are known/postulated 853 to have formed at the same temperature (e.g. `DVH-2` and `DHC2-8` from the `fiebig_2021` 854 and `anderson_2021_lsce` data sets). Each element of `same_T` is a `list` with the names 855 of two or more samples formed at the same temperature. 856 + **exclude_samples**: Use this `list` to specify the names of samples to exclude from 857 the combined calibration. 858 859 For example, the `OGLS23` calibration is computed with: 860 861 `same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]` 862 863 Note that when samples from different calibrations have the same name, 864 it is not necessary to explicitly list them in `same_T`. 865 866 Also note that the regression will fail if samples listed together in `same_T` 867 actually have different `T` values specified in the original calibrations. 868 869 ### Example 870 871 The `devils_laghetto_2023` calibration is computed using the following code: 872 873 ````py 874 K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']] 875 876 fiebig_temp = D47calib( 877 samples = [fiebig_2021.samples[_] for _ in K], 878 T = fiebig_2021.T[K], 879 D47 = fiebig_2021.D47[K], 880 sT = fiebig_2021.sT[K,:][:,K], 881 sD47 = fiebig_2021.sD47[K,:][:,K], 882 ) 883 884 devils_laghetto_2023 = combine_D47calibs( 885 calibs = [anderson_2021_lsce, fiebig_temp], 886 degrees = [0,2], 887 same_T = [{'DVH-2', 'DHC2-8'}], 888 ) 889 ```` 890 ''' 891 892 samples = [s for c in calibs for s in c.samples] 893 T = [t for c in calibs for t in c.T] 894 D47 = [x for c in calibs for x in c.D47] 895 sD47 = _block_diag(*[c.sD47 for c in calibs]) 896 sT = _block_diag(*[c.sT for c in calibs]) 897 898 for i in range(len(samples)): 899 for j in range(len(samples)): 900 if i != j: 901 if (samples[i] == samples[j] or 902 any([samples[i] in _ and samples[j] in _ for _ in same_T])): 903 904 sT[i,j] = (sT[i,i] * sT[j,j])**.5 905 906 k = [_ for _, s in enumerate(samples) if s not in exclude_samples] 907 908 calib = D47calib( 909 samples = [samples[_] for _ in k], 910 T = [T[_] for _ in k], 911 D47 = [D47[_] for _ in k], 912 sT = sT[k,:][:,k], 913 sD47 = sD47[k,:][:,k], 914 degrees = degrees, 915 ) 916 917 return calib 918 919from ._calibs import * 920 921def _covar2correl(C): 922 SE = _np.diag(C)**.5 923 return SE, _np.diag(SE**-1) @ C @ _np.diag(SE**-1) 924 925try: 926 app = typer.Typer( 927 add_completion = False, 928 context_settings={'help_option_names': ['-h', '--help']}, 929 rich_markup_mode = 'rich', 930 ) 931 932 @app.command() 933 def _cli( 934 input: Annotated[str, typer.Argument(help = "Specify either the path of an input file or just '-' to read input from stdin")] = '-', 935 include_samples: Annotated[str, typer.Option('--include-samples', '-u', help = 'Only include samples listed in this file')] = 'all', 936 exclude_samples: Annotated[str, typer.Option('--exclude-samples', '-x', help = 'Exclude samples listed in this file')] = 'none', 937 outfile: Annotated[str, typer.Option('--output-file', '-o', help = 'Write output to this file instead of printing to stdout')] = 'none', 938 calib: Annotated[str, typer.Option('--calib', '-c', help = 'D47 calibration function to use')] = 'OGLS23', 939 delim_in: Annotated[str, typer.Option('--delimiter-in', '-i', help = "Delimiter used in the input.")] = ',', 940 delim_out: Annotated[str, typer.Option('--delimiter-out', '-j', help = "Delimiter used in the output. Use '>' or '<' for elastic white space with right- or left-justified cells.")] = "',' when writing to output file, '>' when printing to stdout", 941 T_precision: Annotated[int, typer.Option('--T-precision', '-p', help = 'Precision for T output')] = 2, 942 D47_precision: Annotated[int, typer.Option('--D47-precision', '-q', help = 'Precision for D47 output')] = 4, 943 correl_precision: Annotated[int, typer.Option('--correl-precision', '-r', help = 'Precision for correlation output')] = 3, 944 covar_precision: Annotated[int, typer.Option('--covar-precision', '-s', help = 'Precision for covariance output')] = 3, 945 return_covar: Annotated[bool, typer.Option('--return-covar', '-v', help = 'Output covariance matrix instead of correlation matrix')] = False, 946 ignore_correl: Annotated[bool, typer.Option('--ignore-correl', '-g', help = 'Only consider and report standard errors without correlations')] = False, 947 uncertainty_sources: Annotated[bool, typer.Option('--uncertainty-sources', '-U', help = 'Output different sources of uncertainty')] = False, 948 version: Annotated[bool, typer.Option('--version', '-V', help = 'Print D47calib version')] = False, 949 ): 950 """ 951[b]Purpose:[/b] 952 953Reads data from an input file, converts between T and D47 values, and print out the results. 954 955The input file is a CSV, or any similar file with data sorted into lines and columns. The line separator must be a <newline>. The column separator, noted <sep> hereafter, is "," by default, or may be any other single character such as ";" or <tab>. 956 957The first line of the input file must be one of the following: 958 959- [b]Option 1:[/b] T 960- [b]Option 2:[/b] T<sep>T_SE 961- [b]Option 3:[/b] T<sep>T_SE<sep>T_correl 962- [b]Option 4:[/b] T<sep>T_covar 963- [b]Option 5:[/b] D47 964- [b]Option 6:[/b] D47<sep>D47_SE 965- [b]Option 7:[/b] D47<sep>D47_SE<sep>D47_correl 966- [b]Option 8:[/b] D47<sep>D47_covar 967 968The rest of the input must be any number of lines with float values corresponding to the fields in the first line, separated by <sep>. 969 970[bold]Example input file:[/bold] 971 972[italic]D47 D47_SE D47_correl[/italic] 973[italic]0.6324 0.0101 1.00 0.25 0.25[/italic] 974[italic]0.6281 0.0087 0.25 1.00 0.25[/italic] 975[italic]0.6385 0.0095 0.25 0.25 1.00[/italic] 976 977The corresponding D47 (options 1-4) or T (options 4-8) values are computed, along with the combined standard errors accounting for both calibration and any (correlated or uncorrelated) uncertainties in the input values. 978 979For options 2-4 and 5-8, which specify standard errors or covariances for the input values, one may obtain (using option `-U`) the separate components of uncertainty from (a) calibration uncertainties and (b) input uncertainties. 980 981The example above will thus result in an output with the following fields: 982 983[italic]- D47[/italic] 984[italic]- D47_SE[/italic] 985[italic]- D47_correl[/italic] 986[italic]- T[/italic] 987[italic]- T_SE_from_calib[/italic] 988[italic]- T_correl_from_calib[/italic] 989[italic]- T_SE_from_input[/italic] 990[italic]- T_correl_from_input[/italic] 991[italic]- T_SE_from_both[/italic] 992[italic]- T_correl_from_both[/italic] 993 994Results may also be saved to a file using [bold]--output-file <filename>[/bold] or [bold]-o <filename>[/bold]. 995 996To filter the samples (lines) to process using [b]--exclude-samples[/b] and [b]--include-samples[/b], first add a [b]Sample[/b] column to the input data, assign a sample name to each line. 997Then to exclude some samples, provide the [b]--exclude-samples[/b] option with the name of a file where each line is one sample to exclude. 998To exclude all samples except those listed in a file, provide the [b]--include-samples[/b] option with the name of that file, where each line is one sample to include. 999""" 1000 1001 if version: 1002 print(__version__) 1003 return None 1004 1005 ### INCOMPATIBILITY BETWEEN --ignore-correl AND --return-covar 1006 if ignore_correl: 1007 return_covar = False 1008 1009 ### USE WHITESPACE AS INPUT DELIMITER 1010 if delim_in == ' ': 1011 delim_in = None 1012 1013 ### SMART SELECTION OF OUTPUT DELIMITER 1014 if delim_out == "',' when writing to output file, '>' when printing to stdout": 1015 if outfile == 'none': 1016 delim_out = '>' 1017 else: 1018 delim_out = ',' 1019 1020 ### CALIBRATION 1021 if calib in globals() and type(globals()[calib]) == D47calib: 1022 calib = globals()[calib] 1023 else: 1024 with open(calib) as f: 1025 calibdata = _np.array([[c.strip() for c in l.strip().split(delim_in)] for l in f.readlines()[1:]], dtype = float) 1026 1027 degrees = [int(d) for d in calibdata[:,0]] 1028 bfp = {f'a{k}': a for k,a in zip(degrees, calibdata[:,1])} 1029 bfp_CM = calibdata[:,2:] 1030 if bfp_CM.shape[0] != bfp_CM.shape[1]: 1031 bfp_CM = _np.zeros((len(degrees), len(degrees))) 1032 calib = D47calib( 1033 samples = [], T = [], sT = [], D47 = [], sD47 = [], 1034 degrees = degrees, bfp = bfp, bfp_CM = bfp_CM, 1035 ) 1036 1037 ### READ INPUT STRINGS 1038 if input == '-': 1039 data = [[c.strip() for c in l.strip().split(delim_in)] for l in sys.stdin] 1040 else: 1041 with open(input) as f: 1042 data = [[c.strip() for c in l.strip().split(delim_in)] for l in f.readlines()] 1043 1044 if include_samples == 'all': 1045 samples_to_include = [] 1046 else: 1047 with open(include_samples) as f: 1048 samples_to_include = [l.strip() for l in f.readlines()] 1049 1050 if exclude_samples == 'none': 1051 samples_to_exclude = [] 1052 else: 1053 with open(exclude_samples) as f: 1054 samples_to_exclude = [l.strip() for l in f.readlines()] 1055 1056 if len(samples_to_include) > 0 or len(samples_to_exclude) > 0: 1057 try: 1058 k = data[0].index('Sample') 1059 except ValueError: 1060 raise KeyError("When using options --include-samples or --exclude-samples, the input file must have a column labeled 'Sample'.") 1061 1062 if len(samples_to_include) > 0: 1063 data = [data[0]] + [l for l in data[1:] if l[k] in samples_to_include] 1064 data = [data[0]] + [l for l in data[1:] if l[k] not in samples_to_exclude] 1065 1066 ### FIND FIRST DATA COLUMN 1067 k = 0 1068 while data[0][k] not in ['T', 'D47']: 1069 k += 1 1070 if k == len(data[0]): 1071 raise KeyError("None of the input column headers are 'T' or 'D47'.") 1072 data_out = [l[:k] for l in data] 1073 data = [l[k:] for l in data] 1074 1075 ### READ INPUT FIELDS 1076 fields = data[0] 1077 1078 ### CHECK FOR UNSUPPORTED FIELD COMBINATIONS 1079 if fields not in [ 1080 ['T'], 1081 ['T', 'T_SE'], 1082 ['T', 'T_covar'], 1083 ['T', 'T_SE', 'T_correl'], 1084 ['D47'], 1085 ['D47', 'D47_SE'], 1086 ['D47', 'D47_covar'], 1087 ['D47', 'D47_SE', 'D47_correl'], 1088 ]: 1089 raise KeyError("There is a problem with the combination of field names provided as input.") 1090 1091 ### BOOK-KEEPING 1092 infield = fields[0] 1093 X_precision = {'T': T_precision, 'D47': D47_precision}[infield] 1094 outfield = {'T': 'D47', 'D47': 'T'}[infield] 1095 Y_precision = {'T': T_precision, 'D47': D47_precision}[outfield] 1096 N = len(data)-1 1097 1098 ### READ INPUT DATA, ALSO SAVING ITS ORIGINAL STRINGS 1099 X_str = [l[0] for l in data[1:]] 1100 X = _np.array(X_str, dtype = float) 1101 1102 if len(fields) == 1: 1103 X_SE = X*0 1104 X_correl = _np.eye(N) 1105 X_covar = _np.zeros((N, N)) 1106 X_SE_str = [f'{c:.{X_precision}f}' for c in X_SE] 1107 X_correl_str = [[f'{c:.{correl_precision}f}' for c in l] for l in X_correl] 1108 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1109 if len(fields) == 2: 1110 if fields[1].endswith('_SE'): 1111 X_SE_str = [l[1] for l in data[1:]] 1112 X_SE = _np.array(X_SE_str, dtype = float) 1113 X_covar = _np.diag(X_SE**2) 1114 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1115 elif fields[1].endswith('_covar'): 1116 X_covar_str = [l[1:N+1] for l in data[1:]] 1117 X_covar = _np.array(X_covar_str, dtype = float) 1118 X_SE = _np.diag(X_covar)**.5 1119 X_SE_str = [f'{c:.{X_precision}f}' for c in X_SE] 1120 X_correl = _np.diag(X_SE**-1) @ X_covar @ _np.diag(X_SE**-1) 1121 X_correl_str = [[f'{c:.{correl_precision}f}' for c in l] for l in X_correl] 1122 elif len(fields) == 3: 1123 X_SE_str = [l[1] for l in data[1:]] 1124 X_SE = _np.array(X_SE_str, dtype = float) 1125 X_correl_str = [l[2:N+2] for l in data[1:]] 1126 X_correl = _np.array(X_correl_str, dtype = float) 1127 X_covar = _np.diag(X_SE) @ X_correl @ _np.diag(X_SE) 1128 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1129 1130 ### COMPUTE OUTPUT VALUES AND COVAR 1131 kwargs = {infield: X, f's{infield}': X_covar} 1132 Y, Y_covar_from_calib = calib.T47(**kwargs, error_from = 'calib', return_covar = True) 1133 Y, Y_covar_from_input = calib.T47(**kwargs, error_from = f's{infield}', return_covar = True) 1134 Y, Y_covar_from_both = calib.T47(**kwargs, error_from = 'both', return_covar = True) 1135 1136 Y_SE_from_calib = _np.diag(Y_covar_from_calib)**.5 1137 Y_SE_from_input = _np.diag(Y_covar_from_input)**.5 1138 Y_SE_from_both = _np.diag(Y_covar_from_both)**.5 1139 1140 if (Y_SE_from_calib**2).min(): 1141 Y_correl_from_calib = _np.diag(Y_SE_from_calib**-1) @ Y_covar_from_calib @ _np.diag(Y_SE_from_calib**-1) 1142 else: 1143 Y_correl_from_calib = _np.eye(N) 1144 1145 if (Y_SE_from_input**2).min(): 1146 Y_correl_from_input = _np.diag(Y_SE_from_input**-1) @ Y_covar_from_input @ _np.diag(Y_SE_from_input**-1) 1147 else: 1148 Y_correl_from_input = _np.eye(N) 1149 1150 if (Y_SE_from_both**2).min(): 1151 Y_correl_from_both = _np.diag(Y_SE_from_both**-1) @ Y_covar_from_both @ _np.diag(Y_SE_from_both**-1) 1152 else: 1153 Y_correl_from_both = _np.eye(N) 1154 1155 ### BUILD Y STRINGS 1156 Y_str = [f'{y:.{Y_precision}f}' for y in Y] 1157 1158 Y_SE_from_calib_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_calib] 1159 Y_SE_from_input_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_input] 1160 Y_SE_from_both_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_both] 1161 1162 Y_covar_from_calib_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_calib] 1163 Y_covar_from_input_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_input] 1164 Y_covar_from_both_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_both] 1165 1166 Y_correl_from_calib_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_calib] 1167 Y_correl_from_input_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_input] 1168 Y_correl_from_both_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_both] 1169 1170 ### ADD SE COLUMN TO INPUT 1171 if f'{infield}_covar' in fields: 1172 fields.insert(1, f'{infield}_SE') 1173 1174 ### ADD X COLUMNS TO OUTPUT DATA 1175 data_out[0] += [infield] 1176 for k in range(N): 1177 data_out[k+1] += [X_str[k]] 1178 for f in fields[1:]: 1179 if f.endswith('_SE'): 1180 data_out[0] += [f] 1181 for k in range(N): 1182 data_out[k+1] += [X_SE_str[k]] 1183 if f.endswith('_covar') or f.endswith('_correl'): 1184 if not ignore_correl: 1185 data_out[0] += [f] + ['' for _ in range(N-1)] 1186 for k in range(N): 1187 data_out[k+1] += (X_covar_str if f.endswith('_covar') else X_correl_str)[k][:] 1188 1189 ### ADD Y COLUMNS TO OUTPUT DATA 1190 data_out[0] += [outfield] 1191 for k in range(N): 1192 data_out[k+1] += [Y_str[k]] 1193 1194 if uncertainty_sources: 1195 data_out[0] += [f'{outfield}_SE_from_calib'] 1196 for k in range(N): 1197 data_out[k+1] += [Y_SE_from_calib_str[k]] 1198 if not ignore_correl: 1199 if return_covar: 1200 data_out[0] += [f'{outfield}_covar_from_calib'] + ['' for _ in range(N-1)] 1201 for k in range(N): 1202 data_out[k+1] += Y_covar_from_calib_str[k] 1203 else: 1204 data_out[0] += [f'{outfield}_correl_from_calib'] + ['' for _ in range(N-1)] 1205 for k in range(N): 1206 data_out[k+1] += Y_correl_from_calib_str[k] 1207 1208 data_out[0] += [f'{outfield}_SE_from_input'] 1209 for k in range(N): 1210 data_out[k+1] += [Y_SE_from_input_str[k]] 1211 if not ignore_correl: 1212 if return_covar: 1213 data_out[0] += [f'{outfield}_covar_from_input'] + ['' for _ in range(N-1)] 1214 for k in range(N): 1215 data_out[k+1] += Y_covar_from_input_str[k] 1216 else: 1217 data_out[0] += [f'{outfield}_correl_from_input'] + ['' for _ in range(N-1)] 1218 for k in range(N): 1219 data_out[k+1] += Y_correl_from_input_str[k] 1220 1221 data_out[0] += [f'{outfield}_SE_from_both' if uncertainty_sources else f'{outfield}_SE'] 1222 for k in range(N): 1223 data_out[k+1] += [Y_SE_from_both_str[k]] 1224 if not ignore_correl: 1225 if return_covar: 1226 data_out[0] += [f'{outfield}_covar_from_both' if uncertainty_sources else f'{outfield}_covar'] + ['' for _ in range(N-1)] 1227 for k in range(N): 1228 data_out[k+1] += Y_covar_from_both_str[k] 1229 else: 1230 data_out[0] += [f'{outfield}_correl_from_both' if uncertainty_sources else f'{outfield}_correl'] + ['' for _ in range(N-1)] 1231 for k in range(N): 1232 data_out[k+1] += Y_correl_from_both_str[k] 1233 1234 1235 ### PRINT OUTPUT TO STDOUT OR SAVE IT TO FILE 1236 if delim_out in '<>': 1237 lengths = [max([len(data_out[j][k]) for j in range(len(data_out))]) for k in range(len(data_out[0]))] 1238 1239 txt = '' 1240 for l in data_out: 1241 for k in range(len(l)): 1242 if k > 0: 1243 txt += ' ' 1244 txt += f'{l[k]:{delim_out}{lengths[k]}s}' 1245 txt += '\n' 1246 1247 txt = txt[:-1] 1248 1249 else: 1250 txt = '\n'.join([delim_out.join(l) for l in data_out]) 1251 1252 if outfile == 'none': 1253 print(txt) 1254 else: 1255 with open(outfile, 'w') as f: 1256 f.write(txt) 1257 1258 def __cli(): 1259 app() 1260 1261except NameError: 1262 pass 1263
Reference
class
D47calib(ogls.InverseTPolynomial):
41class D47calib(_ogls.InverseTPolynomial): 42 """ 43 Δ47 calibration class based on OGLS regression 44 of Δ47 as a polynomial function of inverse T. 45 """ 46 47 def __init__(self, 48 samples, T, D47, 49 sT = None, 50 sD47 = None, 51 degrees = [0,2], 52 xpower = 2, 53 name = '', 54 label = '', 55 description = '', 56 **kwargs, 57 ): 58 """ 59 ### Parameters 60 61 + **samples**: a list of N sample names. 62 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 63 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 64 + **sT**: uncertainties on `T`. If specified as: 65 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 66 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 67 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 68 + **sD47**: uncertainties on `D47`. If specified as: 69 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 70 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 71 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 72 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 73 + **name**: a human-readable, short name assigned to the calibration. 74 + **label**: a short description of the calibration, e.g., to be used in legends. 75 + **description**: a longer description, including relevant references/DOIs. 76 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 77 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 78 79 ### Notable attributes 80 81 + **N**: 82 The total number of observations (samples) in the calibration data. 83 + **samples**: 84 The list sample names. 85 + **T**: 86 1-D `ndarray` of temperatures in degrees C. 87 + **D47**: 88 1-D `ndarray` of Δ47 values in permil. 89 + **sT**: 90 2-D `ndarray` equal to the full (co)variance matrix for `T`. 91 + **D47**: 92 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 93 + **xpower**: 94 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 95 It is possible to change this behavior to use a different power of 1/T. 96 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 97 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 98 + **bfp**: 99 The best-fit parameters of the regression. 100 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 101 + **bff()**: 102 The best-fit polynomial function of inverse T, defined as: 103 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 104 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 105 106 107 ### Examples 108 109 A very simple example: 110 111 ````py 112 .. include:: ../../code_examples/D47calib_init/example.py 113 ```` 114 115 Should yield: 116 117 ```` 118 .. include:: ../../code_examples/D47calib_init/output.txt 119 ```` 120 121 """ 122 123 self.samples = samples[:] 124 self.name = name 125 self.label = label 126 self.description = description 127 self.D47 = _np.asarray(D47, dtype = 'float') 128 self.N = self.D47.size 129 130 if sD47 is None: 131 self.sD47 = _np.zeros((self.N, self.N)) 132 else: 133 self.sD47 = _np.asarray(sD47) 134 if len(self.sD47.shape) == 1: 135 self.sD47 = _np.diag(self.sD47**2) 136 elif len(self.sD47.shape) == 0: 137 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 138 139 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 140 141 if self.bfp is None: 142 self.regress() 143 144 self._bff_deriv = lambda x: _np.array([k * self.bfp[f'a{k}'] * x**(k-1) for k in degrees if k > 0]).sum(axis = 0) 145 146 xi = _np.linspace(0,200**-1,1001) 147 self._inv_bff = _interp1d(self.bff(xi), xi) 148 149 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 150 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 151 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 152 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1 153 154 def __repr__(self): 155 return f'<D47calib: {self.name}>' 156 157 def invT_xaxis(self, 158 xlabel = None, 159 Ti = [0,20,50,100,250,1000], 160 ): 161 """ 162 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 163 but labeled in degrees Celsius. 164 165 ### Parameters 166 167 + **xlabel**: 168 Custom label for X axis (`r'$1\\,/\\,T^2$'` by default) 169 + **Ti**: 170 Specify tick locations for X axis, in degrees C. 171 172 ### Returns 173 174 + an `matplotlib.axes.Axes` instance 175 176 ### Examples 177 178 ````py 179 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 180 ```` 181 182 This should result in something like this: 183 184 <img align="center" src="example_invT_xaxis_1.png"> 185 186 It is also possible to define the X axis using a different power of 1/T 187 by first redefining the `xpower` attribute: 188 189 ````py 190 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 191 ```` 192 193 This should result in something like this: 194 195 <img align="center" src="example_invT_xaxis_2.png"> 196 """ 197 if xlabel is None: 198 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 199 _ppl.xlabel(xlabel) 200 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 201 ax = _ppl.gca() 202 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 203 ax.tick_params(which="major") 204 205 return ax 206 207 208 def plot_data(self, label = False, **kwargs): 209 """ 210 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 211 212 ### Parameters 213 214 + **label**: 215 + If `label` is a string, use this string as `label` for the underlyig 216 `matplotlib.pyplot.plot()` call. 217 + If `label = True`, use the caller's `label` attribute instead. 218 + If `label = False`, no label is specified (default behavior). 219 + **kwargs**: 220 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 221 222 ### Returns 223 224 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 225 226 ### Example 227 228 ````py 229 from matplotlib import pyplot as ppl 230 from D47calib import huyghe_2022 as calib 231 232 fig = ppl.figure(figsize = (5,3)) 233 ppl.subplots_adjust(bottom = .25, left = .15) 234 calib.invT_xaxis(Ti = [0,10,25]) 235 calib.plot_data(label = True) 236 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 237 ppl.legend() 238 ppl.savefig('example_plot_data.png', dpi = 100) 239 ````` 240 241 This should result in something like this: 242 243 <img align="center" src="example_plot_data.png"> 244 """ 245# if 'mec' not in kwargs: 246# kwargs['mec'] = self.color 247 if label is not False: 248 kwargs['label'] = self.label if label is True else label 249 return _ogls.InverseTPolynomial.plot_data(self, **kwargs) 250 251 252 def plot_error_bars(self, **kwargs): 253 """ 254 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 255 256 ### Parameters 257 258 + **kwargs**: 259 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 260 261 ### Returns 262 263 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 264 265 ### Example 266 267 ````py 268 from matplotlib import pyplot as ppl 269 from D47calib import huyghe_2022 as calib 270 271 fig = ppl.figure(figsize = (5,3)) 272 ppl.subplots_adjust(bottom = .25, left = .15) 273 calib.invT_xaxis(Ti = [0,10,25]) 274 calib.plot_error_bars(alpha = .4) 275 calib.plot_data(label = True) 276 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 277 ppl.legend() 278 ppl.savefig('example_plot_error_bars.png', dpi = 100) 279 ````` 280 281 This should result in something like this: 282 283 <img align="center" src="example_plot_error_bars.png"> 284 """ 285# if 'ecolor' not in kwargs: 286# kwargs['ecolor'] = self.color 287 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs) 288 289 290 def plot_error_ellipses(self, **kwargs): 291 """ 292 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 293 294 ### Parameters 295 296 + **kwargs**: 297 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 298 299 ### Returns 300 301 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 302 303 ### Example 304 305 ````py 306 from matplotlib import pyplot as ppl 307 from D47calib import huyghe_2022 as calib 308 309 fig = ppl.figure(figsize = (5,3)) 310 ppl.subplots_adjust(bottom = .25, left = .15) 311 calib.invT_xaxis(Ti = [0,10,25]) 312 calib.plot_error_ellipses(alpha = .4) 313 calib.plot_data(label = True) 314 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 315 ppl.legend() 316 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 317 ````` 318 319 This should result in something like this: 320 321 <img align="center" src="example_plot_error_ellipses.png"> 322 """ 323# if 'ec' not in kwargs: 324# kwargs['ec'] = self.color 325 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs) 326 327 328 def plot_bff(self, label = False, **kwargs): 329 """ 330 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 331 332 ### Parameters 333 334 + **label**: 335 + If `label` is a string, use this string as `label` for the underlyig 336 `matplotlib.pyplot.plot()` call. 337 + If `label = True`, use the caller's `label` attribute instead. 338 + If `label = False`, no label is specified (default behavior). 339 + **kwargs**: 340 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 341 342 ### Returns 343 344 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 345 346 ### Example 347 348 ````py 349 from matplotlib import pyplot as ppl 350 from D47calib import huyghe_2022 as calib 351 352 fig = ppl.figure(figsize = (5,3)) 353 ppl.subplots_adjust(bottom = .25, left = .15) 354 calib.invT_xaxis(Ti = [0,10,25]) 355 calib.plot_bff(label = True, dashes = (8,2,2,2)) 356 calib.plot_data() 357 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 358 ppl.legend() 359 ppl.savefig('example_plot_bff.png', dpi = 100) 360 ````` 361 362 This should result in something like this: 363 364 <img align="center" src="example_plot_bff.png"> 365 """ 366# if 'color' not in kwargs: 367# kwargs['color'] = self.color 368 if label is not False: 369 kwargs['label'] = self.label if label is True else label 370 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs) 371 372 373 def plot_bff_ci(self, **kwargs): 374 """ 375 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 376 377 ### Parameters 378 379 + **label**: 380 + **kwargs**: 381 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 382 383 ### Returns 384 385 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 386 387 ### Example 388 389 ````py 390 from matplotlib import pyplot as ppl 391 from D47calib import huyghe_2022 as calib 392 393 fig = ppl.figure(figsize = (5,3)) 394 ppl.subplots_adjust(bottom = .25, left = .15) 395 calib.invT_xaxis(Ti = [0,10,25]) 396 calib.plot_bff_ci(alpha = .15) 397 calib.plot_bff(label = True, dashes = (8,2,2,2)) 398 calib.plot_data() 399 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 400 ppl.legend() 401 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 402 ````` 403 404 This should result in something like this: 405 406 <img align="center" src="example_plot_bff_ci.png"> 407 """ 408# if 'color' not in kwargs: 409# kwargs['color'] = self.color 410 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs) 411 412 def T47(self, 413 D47 = None, 414 sD47 = None, 415 T=None, 416 sT = None, 417 error_from = 'both', 418 return_covar = False, 419 ): 420 ''' 421 When `D47` is input, computes corresponding T value(s). 422 `D47` input may be specified as a scalar, or as a 1-D array. 423 `T` output will then have the same type and size as `D47`. 424 425 When `T` is input, computes corresponding Δ47 value(s). 426 `T` input may be specified as a scalar, or as a 1-D array. 427 `D47` output will then have the same type and size as `T`. 428 429 Only one of either `D47` or `T` may be specified as input. 430 431 **Arguments:** 432 433 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 434 * `sD47`: Δ47 uncertainties, which may be: 435 - `None` (default) 436 - `float` or `int` (uniform standard error on `D47`) 437 - 1-D array (standard errors on `D47`) 438 - 2-D array (covariance matrix for `D47`) 439 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 440 * `sT`: T uncertainties, which may be: 441 - `None` (default) 442 - `float` or `int` (uniform standard error on `T`) 443 - 1-D array (standard errors on `T`) 444 - 2-D array (variance-covariance matrix for `T`) 445 * `error_from`: if set to `'both'` (default), returned errors take into account 446 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 447 if set to `'calib'`, only calibration uncertainties are accounted for; 448 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 449 * `return_covar`: (False by default) whether to return the full covariance matrix 450 for returned `T` or `D47` values, otherwise return standard errors for the returned 451 `T` or `D47` values instead. 452 453 **Returns (with `D47` input):** 454 455 * `T`: temperature value(s) computed from `D47` 456 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 457 458 **Returns (with `T` input):** 459 460 * `D47`: Δ47 value(s) computed from `D47` 461 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 462 463 ### Example 464 465 ````py 466 import numpy as np 467 from matplotlib import pyplot as ppl 468 from D47calib import OGLS23 as calib 469 470 X = np.linspace(1473**-2, 270**-2) 471 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 472 473 fig = ppl.figure(figsize = (5,3)) 474 ppl.subplots_adjust(bottom = .25, left = .15) 475 calib.invT_xaxis() 476 ppl.plot(X, 1000 * sD47, 'r-') 477 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 478 ppl.savefig('example_SE47.png', dpi = 100) 479 ````` 480 481 This should result in something like this: 482 483 <img src="example_SE47.png"> 484 ''' 485 486 if D47 is None and T is None: 487 raise ValueError('Either D47 or T must be specified, but both are undefined.') 488 489 if D47 is not None and T is not None: 490 raise ValueError('Either D47 or T must be specified, but not both.') 491 492 if T is not None: 493 494 D47 = self._D47_from_T(T) 495 Np = len(self.degrees) 496 N = D47.size 497 498 ### Compute covariance matrix of (*bfp, *T): 499 CM = _np.zeros((Np+N, Np+N)) 500 501 if error_from in ['calib', 'both']: 502 CM[:Np, :Np] = self.bfp_CM[:,:] 503 504 if (sT is not None) and error_from in ['sT', 'both']: 505 _sT = _np.asarray(sT) 506 if _sT.ndim == 0: 507 for k in range(N): 508 CM[Np+k, Np+k] = _sT**2 509 elif _sT.ndim == 1: 510 for k in range(N): 511 CM[Np+k, Np+k] = _sT[k]**2 512 elif _sT.ndim == 2: 513 CM[-N:, -N:] = _sT[:,:] 514 515 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 516 _T = _np.asarray(T) 517 if _T.ndim == 0: 518 _T = _np.expand_dims(_T, 0) 519 J = _np.zeros((N, Np+N)) 520 521 if (sT is not None) and error_from in ['sT', 'both']: 522 for k in range(N): 523 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 524 525 if error_from in ['calib', 'both']: 526 527 for k in range(Np): 528 529 p1 = {_: self.bfp[_] for _ in self.bfp} 530 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 531 532 p2 = {_: self.bfp[_] for _ in self.bfp} 533 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 534 535 J[:, k] = (self.model_fun(p1, (_T+273.15)**-1) - self.model_fun(p2, (_T+273.15)**-1)) / (0.002 * self.bfp_CM[k,k]**.5) 536 537 ### Error propagation: 538 CM_D47 = J @ CM @ J.T 539 540 if return_covar: 541 return D47, CM_D47 542 else: 543 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 544 545 if D47 is not None: 546 547 T = self._T_from_D47(D47) 548 Np = len(self.degrees) 549 N = T.size 550 551 ### Compute covariance matrix of (*bfp, *T): 552 CM = _np.zeros((Np+N, Np+N)) 553 554 if error_from in ['calib', 'both']: 555 CM[:Np, :Np] = self.bfp_CM[:,:] 556 557 if (sD47 is not None) and error_from in ['sD47', 'both']: 558 _sD47 = _np.asarray(sD47) 559 if _sD47.ndim == 0: 560 for k in range(N): 561 CM[Np+k, Np+k] = _sD47**2 562 elif _sD47.ndim == 1: 563 for k in range(N): 564 CM[Np+k, Np+k] = _sD47[k]**2 565 elif _sD47.ndim == 2: 566 CM[-N:, -N:] = _sD47[:,:] 567 568 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 569 _D47 = _np.asarray(D47) 570 if _D47.ndim == 0: 571 _D47 = _np.expand_dims(_D47, 0) 572 J = _np.zeros((N, Np+N)) 573 if (sD47 is not None) and error_from in ['sD47', 'both']: 574 for k in range(N): 575 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 576 if error_from in ['calib', 'both']: 577 578 xi = _np.linspace(0,200**-1,1001)[1:] 579 for k in range(Np): 580 581 if self.bfp_CM[k,k]: 582 _epsilon_ = self.bfp_CM[k,k]**.5 583 else: 584 _epsilon_ = 1e-6 585 586 p1 = {_: self.bfp[_] for _ in self.bfp} 587 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 588 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 589 590 p2 = {_: self.bfp[_] for _ in self.bfp} 591 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 592 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 593 594 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 595 596 ### Error propagation: 597 CM_T = J @ CM @ J.T 598 599 if return_covar: 600 return T, CM_T 601 else: 602 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5 603 604 605 def plot_T47_errors( 606 self, 607 calibname = None, 608 rD47 = 0.010, 609 Nr = [2,4,8,12,20], 610 Tmin = 0, 611 Tmax = 120, 612 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 613 yscale = 'lin', 614 ): 615 """ 616 Plot SE of T reconstructed using the calibration as a function of T for various 617 combinations of analytical precision and number of analytical replicates. 618 619 **Arguments** 620 621 + **calibname**: 622 Which calibration name to display. By default, use `label` attribute. 623 + **rD47**: 624 Analytical precision of a single analysis. 625 + **Nr**: 626 A list of lines to plot, each corresponding to a given number of replicates. 627 + **Tmin**: 628 Minimum T to plot. 629 + **Tmax**: 630 Maximum T to plot. 631 + **colors**: 632 A list of colors to distinguish the plotted lines. 633 + **yscale**: 634 + If `'lin'`, the Y axis uses a linear scale. 635 + If `'log'`, the Y axis uses a logarithmic scale. 636 637 **Example** 638 639 ````py 640 from matplotlib import pyplot as ppl 641 from D47calib import devils_laghetto_2023 as calib 642 643 fig = ppl.figure(figsize = (3.5,4)) 644 ppl.subplots_adjust(bottom = .2, left = .15) 645 calib.plot_T47_errors( 646 calibname = 'Devils Laghetto calibration', 647 Nr = [1,2,4,16], 648 Tmin =0, 649 Tmax = 40, 650 ) 651 ppl.savefig('example_SE_T.png', dpi = 100) 652 ```` 653 654 This should result in something like this: 655 656 <img src="example_SE_T.png"> 657 """ 658 659 if calibname is None: 660 calibname = self.label 661 662 Nr = _np.array(Nr) 663 if len(colors) < Nr.size: 664 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 665 from colorsys import hsv_to_rgb 666 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 667 colors = [hsv_to_rgb(*x) for x in hsv] 668 669 Ti = _np.linspace(Tmin, Tmax) 670 D47i, _ = self.T47(T = Ti) 671 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 672 673 ymax, ymin = 0, 1e6 674 for N,c in zip(Nr, colors): 675 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 676 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 677 ymin = min(ymin, min(sT)) 678 ymax = max(ymax, max(sT)) 679 680 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 681 682 _ppl.legend(fontsize=9) 683 _ppl.xlabel("T (°C)") 684 685 _ppl.ylabel("Standard error on reconstructed T (°C)") 686 687 # yticks([0,.5,1,1.5,2]) 688 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 689 _ppl.grid( alpha = .25) 690 if yscale == 'lin': 691 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 692 t1, t2 = self.T.min(), self.T.max() 693 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 694 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 695 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 696 elif yscale == 'log': 697 ymin /= 2 698 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 699 _ppl.yscale('log') 700 t1, t2 = self.T.min(), self.T.max() 701 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 702 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 703 704 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 705 """ 706 Write calibration data to a csv file. 707 708 ### Parameters 709 710 + **csvfile**: 711 The filename to write data to. 712 + **sep**: 713 The separator between CSV fields. 714 + **label**: 715 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 716 + If specified as a `str`, include a `Dataset` column with that string. 717 + If specified as `False`, do not include a `Dataset` column. 718 + **T_correl**: 719 + If `True`, include correlations between all `T` values. 720 + **D47_correl**: 721 + If `True`, include correlations between all `D47` values. 722 723 ### Example 724 725 ````py 726 D47calib.huyghe_2022.export_data( 727 csvfile = 'example_export_data.csv', 728 T_correl = True, 729 D47_correl = True, 730 ) 731 ```` 732 733 This should result in something like this ([link](example_export_data.csv)): 734 735 .. include:: ../../docs/example_export_data.md 736 737 """ 738 n = len(str(self.N)) 739 740 with open(csvfile, 'w') as f: 741 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 742 743 if label: 744 f.write(f'{sep}Dataset') 745 746 if T_correl: 747 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 748 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 749 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 750 751 if D47_correl: 752 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 753 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 754 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 755 756 for k, (s, T, sT, D47, sD47) in enumerate(zip( 757 self.samples, 758 self.T, 759 _np.diag(self.sT)**.5, 760 self.D47, 761 _np.diag(self.sD47)**.5, 762 )): 763 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 764 if label: 765 if label is True: 766 f.write(f'{sep}{self.label}') 767 else: 768 f.write(f'{sep}{label}') 769 if T_correl: 770 f.write(sep.join(['']+[ 771 f'{Tcorrel[k,_]:.0f}' 772 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 773 else f'{Tcorrel[k,_]:.6f}' 774 for _ in range(self.N)])) 775 if D47_correl: 776 f.write(sep.join(['']+[ 777 f'{D47correl[k,_]:.0f}' 778 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 779 else f'{D47correl[k,_]:.6f}' 780 for _ in range(self.N)])) 781 782 783 def export(self, name, filename): 784 """ 785 Save `D47calib` object as an importable file. 786 787 ### Parameters 788 789 + **name**: 790 The name of the variable to export. 791 + **filename**: 792 The filename to write to. 793 794 ### Example 795 796 ````py 797 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 798 ```` 799 800 This should result in a `bar.py` file with the following contents: 801 802 ````py 803 foo = D47calib( 804 samples = ['LGB-2', 'DVH-2'], 805 T = [7.9, 33.7], 806 D47 = [0.6485720997671647, 0.5695972909966959], 807 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 808 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 809 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 810 label = 'Slow-growing calcites from Anderson et al. (2021)', 811 color = (0, 0.5, 0), 812 degrees = [0, 2], 813 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 814 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 815 chisq = 6.421311854486162e-27, 816 Nf = 0, 817 ) 818 ```` 819 """ 820 with open(filename, 'w') as f: 821 f.write(f''' 822{name} = D47calib( 823 samples = {self.samples}, 824 T = {self.T.tolist()}, 825 D47 = {self.D47.tolist()}, 826 sT = {self.sT.tolist()}, 827 sD47 = {self.sD47.tolist()}, 828 degrees = {self.degrees}, 829 description = {repr(self.description)}, 830 name = {repr(self.name)}, 831 label = {repr(self.label)}, 832 bfp = {({k: float(self.bfp[k]) for k in self.bfp})}, 833 bfp_CM = {self.bfp_CM.tolist()}, 834 chisq = {self.chisq}, 835 cholesky_residuals = {self.cholesky_residuals.tolist()}, 836 aic = {self.aic}, 837 bic = {self.bic}, 838 ks_pvalue = {self.ks_pvalue}, 839 ) 840''')
Δ47 calibration class based on OGLS regression of Δ47 as a polynomial function of inverse T.
D47calib( samples, T, D47, sT=None, sD47=None, degrees=[0, 2], xpower=2, name='', label='', description='', **kwargs)
47 def __init__(self, 48 samples, T, D47, 49 sT = None, 50 sD47 = None, 51 degrees = [0,2], 52 xpower = 2, 53 name = '', 54 label = '', 55 description = '', 56 **kwargs, 57 ): 58 """ 59 ### Parameters 60 61 + **samples**: a list of N sample names. 62 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 63 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 64 + **sT**: uncertainties on `T`. If specified as: 65 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 66 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 67 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 68 + **sD47**: uncertainties on `D47`. If specified as: 69 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 70 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 71 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 72 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 73 + **name**: a human-readable, short name assigned to the calibration. 74 + **label**: a short description of the calibration, e.g., to be used in legends. 75 + **description**: a longer description, including relevant references/DOIs. 76 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 77 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 78 79 ### Notable attributes 80 81 + **N**: 82 The total number of observations (samples) in the calibration data. 83 + **samples**: 84 The list sample names. 85 + **T**: 86 1-D `ndarray` of temperatures in degrees C. 87 + **D47**: 88 1-D `ndarray` of Δ47 values in permil. 89 + **sT**: 90 2-D `ndarray` equal to the full (co)variance matrix for `T`. 91 + **D47**: 92 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 93 + **xpower**: 94 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 95 It is possible to change this behavior to use a different power of 1/T. 96 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 97 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 98 + **bfp**: 99 The best-fit parameters of the regression. 100 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 101 + **bff()**: 102 The best-fit polynomial function of inverse T, defined as: 103 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 104 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 105 106 107 ### Examples 108 109 A very simple example: 110 111 ````py 112 .. include:: ../../code_examples/D47calib_init/example.py 113 ```` 114 115 Should yield: 116 117 ```` 118 .. include:: ../../code_examples/D47calib_init/output.txt 119 ```` 120 121 """ 122 123 self.samples = samples[:] 124 self.name = name 125 self.label = label 126 self.description = description 127 self.D47 = _np.asarray(D47, dtype = 'float') 128 self.N = self.D47.size 129 130 if sD47 is None: 131 self.sD47 = _np.zeros((self.N, self.N)) 132 else: 133 self.sD47 = _np.asarray(sD47) 134 if len(self.sD47.shape) == 1: 135 self.sD47 = _np.diag(self.sD47**2) 136 elif len(self.sD47.shape) == 0: 137 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 138 139 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 140 141 if self.bfp is None: 142 self.regress() 143 144 self._bff_deriv = lambda x: _np.array([k * self.bfp[f'a{k}'] * x**(k-1) for k in degrees if k > 0]).sum(axis = 0) 145 146 xi = _np.linspace(0,200**-1,1001) 147 self._inv_bff = _interp1d(self.bff(xi), xi) 148 149 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 150 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 151 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 152 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1
Parameters
- samples: a list of N sample names.
- T: a 1-D array (or array-like) of temperatures values (in degrees C), of size N.
- D47: a 1-D array (or array-like) of Δ47 values (in permil), of size N.
- sT: uncertainties on
T. If specified as: - sD47: uncertainties on
D47. If specified as: - degrees: degrees of the polynomial regression, e.g.,
[0, 2]or[0, 1, 2, 3, 4]. - name: a human-readable, short name assigned to the calibration.
- label: a short description of the calibration, e.g., to be used in legends.
- description: a longer description, including relevant references/DOIs.
This is not necessary when
bfpandCM_bfpare specified at instantiation time. - kwargs: keyword arguments passed to the underlying
ogls.InverseTPolynomial()call.
Notable attributes
- N: The total number of observations (samples) in the calibration data.
- samples: The list sample names.
- T:
1-D
ndarrayof temperatures in degrees C. - D47:
1-D
ndarrayof Δ47 values in permil. - sT:
2-D
ndarrayequal to the full (co)variance matrix forT. - D47:
2-D
ndarrayequal to the full (co)variance matrix forD47. - xpower:
By default, all
D47calibgraphical methods plot Δ47 as a function of 1/T2. It is possible to change this behavior to use a different power of 1/T. This is done by redefining thexpowerattribute to a different, non-zerointvalue (e.g.foo.xpower = 1to plot as a function of 1/T instead of 1/T2). - bfp:
The best-fit parameters of the regression.
This is a
dictwith keys equal to the polynomial coefficients (seebffdefinition below) - bff():
The best-fit polynomial function of inverse T, defined as:
bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)Note thatbfftakesx = 1/(T+273.15)(instead ofT) as input.
Examples
A very simple example:
from D47calib import D47calib
mycalib = D47calib(
samples = ['FOO', 'BAR'],
T = [0. , 25. ],
D47 = [0.7 , 0.6 ],
sT = 1.,
sD47 = 0.01,
)
T, sT = mycalib.T47(D47 = 0.650)
print(f'T = {T:.1f}')
print(f'sT = {sT:.1f}')
Should yield:
T = 11.7
sT = 1.9
def
invT_xaxis(self, xlabel=None, Ti=[0, 20, 50, 100, 250, 1000]):
157 def invT_xaxis(self, 158 xlabel = None, 159 Ti = [0,20,50,100,250,1000], 160 ): 161 """ 162 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 163 but labeled in degrees Celsius. 164 165 ### Parameters 166 167 + **xlabel**: 168 Custom label for X axis (`r'$1\\,/\\,T^2$'` by default) 169 + **Ti**: 170 Specify tick locations for X axis, in degrees C. 171 172 ### Returns 173 174 + an `matplotlib.axes.Axes` instance 175 176 ### Examples 177 178 ````py 179 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 180 ```` 181 182 This should result in something like this: 183 184 <img align="center" src="example_invT_xaxis_1.png"> 185 186 It is also possible to define the X axis using a different power of 1/T 187 by first redefining the `xpower` attribute: 188 189 ````py 190 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 191 ```` 192 193 This should result in something like this: 194 195 <img align="center" src="example_invT_xaxis_2.png"> 196 """ 197 if xlabel is None: 198 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 199 _ppl.xlabel(xlabel) 200 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 201 ax = _ppl.gca() 202 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 203 ax.tick_params(which="major") 204 205 return ax
Create and return an Axes object with X values equal to 1/T2,
but labeled in degrees Celsius.
Parameters
- xlabel:
Custom label for X axis (
r'$1\,/\,T^2$'by default) - Ti: Specify tick locations for X axis, in degrees C.
Returns
- an
matplotlib.axes.Axesinstance
Examples
from matplotlib import pyplot as ppl
from D47calib import OGLS23 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
ax = calib.invT_xaxis()
ax.set_xlim((0, 270**-2))
ppl.savefig('example_invT_xaxis_1.png', dpi = 100)
This should result in something like this:
It is also possible to define the X axis using a different power of 1/T
by first redefining the xpower attribute:
from matplotlib import pyplot as ppl
from D47calib import OGLS23 as calib
calib.xpower = 4
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
ax = calib.invT_xaxis(Ti = [1000, 100, 50, 25, 0])
ax.set_xlim((0, 270**-4))
ppl.savefig('example_invT_xaxis_2.png', dpi = 100)
This should result in something like this:
def
plot_data(self, label=False, **kwargs):
208 def plot_data(self, label = False, **kwargs): 209 """ 210 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 211 212 ### Parameters 213 214 + **label**: 215 + If `label` is a string, use this string as `label` for the underlyig 216 `matplotlib.pyplot.plot()` call. 217 + If `label = True`, use the caller's `label` attribute instead. 218 + If `label = False`, no label is specified (default behavior). 219 + **kwargs**: 220 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 221 222 ### Returns 223 224 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 225 226 ### Example 227 228 ````py 229 from matplotlib import pyplot as ppl 230 from D47calib import huyghe_2022 as calib 231 232 fig = ppl.figure(figsize = (5,3)) 233 ppl.subplots_adjust(bottom = .25, left = .15) 234 calib.invT_xaxis(Ti = [0,10,25]) 235 calib.plot_data(label = True) 236 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 237 ppl.legend() 238 ppl.savefig('example_plot_data.png', dpi = 100) 239 ````` 240 241 This should result in something like this: 242 243 <img align="center" src="example_plot_data.png"> 244 """ 245# if 'mec' not in kwargs: 246# kwargs['mec'] = self.color 247 if label is not False: 248 kwargs['label'] = self.label if label is True else label 249 return _ogls.InverseTPolynomial.plot_data(self, **kwargs)
Plot Δ47 value of each sample as a function of 1/T2.
Parameters
- label:
- kwargs:
keyword arguments passed to the underlying
matplotlib.pyplot.plot()call.
Returns
- the return value(s) of the underlying
matplotlib.pyplot.plot()call.
Example
from matplotlib import pyplot as ppl
from D47calib import huyghe_2022 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis(Ti = [0,10,25])
calib.plot_data(label = True)
ppl.ylabel('$Δ_{47}$ (‰ I-CDES)')
ppl.legend()
ppl.savefig('example_plot_data.png', dpi = 100)
This should result in something like this:
def
plot_error_bars(self, **kwargs):
252 def plot_error_bars(self, **kwargs): 253 """ 254 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 255 256 ### Parameters 257 258 + **kwargs**: 259 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 260 261 ### Returns 262 263 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 264 265 ### Example 266 267 ````py 268 from matplotlib import pyplot as ppl 269 from D47calib import huyghe_2022 as calib 270 271 fig = ppl.figure(figsize = (5,3)) 272 ppl.subplots_adjust(bottom = .25, left = .15) 273 calib.invT_xaxis(Ti = [0,10,25]) 274 calib.plot_error_bars(alpha = .4) 275 calib.plot_data(label = True) 276 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 277 ppl.legend() 278 ppl.savefig('example_plot_error_bars.png', dpi = 100) 279 ````` 280 281 This should result in something like this: 282 283 <img align="center" src="example_plot_error_bars.png"> 284 """ 285# if 'ecolor' not in kwargs: 286# kwargs['ecolor'] = self.color 287 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs)
Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T2.
Parameters
- kwargs:
keyword arguments passed to the underlying
matplotlib.pyplot.errrobar()call.
Returns
- the return value(s) of the underlying
matplotlib.pyplot.errorbar()call.
Example
from matplotlib import pyplot as ppl
from D47calib import huyghe_2022 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis(Ti = [0,10,25])
calib.plot_error_bars(alpha = .4)
calib.plot_data(label = True)
ppl.ylabel('$Δ_{47}$ (‰ I-CDES)')
ppl.legend()
ppl.savefig('example_plot_error_bars.png', dpi = 100)
This should result in something like this:
def
plot_error_ellipses(self, **kwargs):
290 def plot_error_ellipses(self, **kwargs): 291 """ 292 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 293 294 ### Parameters 295 296 + **kwargs**: 297 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 298 299 ### Returns 300 301 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 302 303 ### Example 304 305 ````py 306 from matplotlib import pyplot as ppl 307 from D47calib import huyghe_2022 as calib 308 309 fig = ppl.figure(figsize = (5,3)) 310 ppl.subplots_adjust(bottom = .25, left = .15) 311 calib.invT_xaxis(Ti = [0,10,25]) 312 calib.plot_error_ellipses(alpha = .4) 313 calib.plot_data(label = True) 314 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 315 ppl.legend() 316 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 317 ````` 318 319 This should result in something like this: 320 321 <img align="center" src="example_plot_error_ellipses.png"> 322 """ 323# if 'ec' not in kwargs: 324# kwargs['ec'] = self.color 325 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs)
Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T2.
Parameters
- kwargs:
keyword arguments passed to the underlying
matplotlib.patches.Ellipse()call.
Returns
- the return value(s) of the underlying
matplotlib.patches.Ellipse()call.
Example
from matplotlib import pyplot as ppl
from D47calib import huyghe_2022 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis(Ti = [0,10,25])
calib.plot_error_ellipses(alpha = .4)
calib.plot_data(label = True)
ppl.ylabel('$Δ_{47}$ (‰ I-CDES)')
ppl.legend()
ppl.savefig('example_plot_error_ellipses.png', dpi = 100)
This should result in something like this:
def
plot_bff(self, label=False, **kwargs):
328 def plot_bff(self, label = False, **kwargs): 329 """ 330 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 331 332 ### Parameters 333 334 + **label**: 335 + If `label` is a string, use this string as `label` for the underlyig 336 `matplotlib.pyplot.plot()` call. 337 + If `label = True`, use the caller's `label` attribute instead. 338 + If `label = False`, no label is specified (default behavior). 339 + **kwargs**: 340 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 341 342 ### Returns 343 344 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 345 346 ### Example 347 348 ````py 349 from matplotlib import pyplot as ppl 350 from D47calib import huyghe_2022 as calib 351 352 fig = ppl.figure(figsize = (5,3)) 353 ppl.subplots_adjust(bottom = .25, left = .15) 354 calib.invT_xaxis(Ti = [0,10,25]) 355 calib.plot_bff(label = True, dashes = (8,2,2,2)) 356 calib.plot_data() 357 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 358 ppl.legend() 359 ppl.savefig('example_plot_bff.png', dpi = 100) 360 ````` 361 362 This should result in something like this: 363 364 <img align="center" src="example_plot_bff.png"> 365 """ 366# if 'color' not in kwargs: 367# kwargs['color'] = self.color 368 if label is not False: 369 kwargs['label'] = self.label if label is True else label 370 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs)
Plot best-fit regression of Δ47 as a function of 1/T2.
Parameters
- label:
- kwargs:
keyword arguments passed to the underlying
matplotlib.pyplot.plot()call.
Returns
- the return value(s) of the underlying
matplotlib.pyplot.plot()call.
Example
from matplotlib import pyplot as ppl
from D47calib import huyghe_2022 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis(Ti = [0,10,25])
calib.plot_bff(label = True, dashes = (8,2,2,2))
calib.plot_data()
ppl.ylabel('$Δ_{47}$ (‰ I-CDES)')
ppl.legend()
ppl.savefig('example_plot_bff.png', dpi = 100)
This should result in something like this:
def
plot_bff_ci(self, **kwargs):
373 def plot_bff_ci(self, **kwargs): 374 """ 375 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 376 377 ### Parameters 378 379 + **label**: 380 + **kwargs**: 381 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 382 383 ### Returns 384 385 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 386 387 ### Example 388 389 ````py 390 from matplotlib import pyplot as ppl 391 from D47calib import huyghe_2022 as calib 392 393 fig = ppl.figure(figsize = (5,3)) 394 ppl.subplots_adjust(bottom = .25, left = .15) 395 calib.invT_xaxis(Ti = [0,10,25]) 396 calib.plot_bff_ci(alpha = .15) 397 calib.plot_bff(label = True, dashes = (8,2,2,2)) 398 calib.plot_data() 399 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 400 ppl.legend() 401 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 402 ````` 403 404 This should result in something like this: 405 406 <img align="center" src="example_plot_bff_ci.png"> 407 """ 408# if 'color' not in kwargs: 409# kwargs['color'] = self.color 410 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs)
Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T2.
Parameters
- label:
- kwargs:
keyword arguments passed to the underlying
matplotlib.pyplot.fill_between()call.
Returns
- the return value(s) of the underlying
matplotlib.pyplot.fill_between()call.
Example
from matplotlib import pyplot as ppl
from D47calib import huyghe_2022 as calib
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis(Ti = [0,10,25])
calib.plot_bff_ci(alpha = .15)
calib.plot_bff(label = True, dashes = (8,2,2,2))
calib.plot_data()
ppl.ylabel('$Δ_{47}$ (‰ I-CDES)')
ppl.legend()
ppl.savefig('example_plot_bff_ci.png', dpi = 100)
This should result in something like this:
def
T47( self, D47=None, sD47=None, T=None, sT=None, error_from='both', return_covar=False):
412 def T47(self, 413 D47 = None, 414 sD47 = None, 415 T=None, 416 sT = None, 417 error_from = 'both', 418 return_covar = False, 419 ): 420 ''' 421 When `D47` is input, computes corresponding T value(s). 422 `D47` input may be specified as a scalar, or as a 1-D array. 423 `T` output will then have the same type and size as `D47`. 424 425 When `T` is input, computes corresponding Δ47 value(s). 426 `T` input may be specified as a scalar, or as a 1-D array. 427 `D47` output will then have the same type and size as `T`. 428 429 Only one of either `D47` or `T` may be specified as input. 430 431 **Arguments:** 432 433 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 434 * `sD47`: Δ47 uncertainties, which may be: 435 - `None` (default) 436 - `float` or `int` (uniform standard error on `D47`) 437 - 1-D array (standard errors on `D47`) 438 - 2-D array (covariance matrix for `D47`) 439 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 440 * `sT`: T uncertainties, which may be: 441 - `None` (default) 442 - `float` or `int` (uniform standard error on `T`) 443 - 1-D array (standard errors on `T`) 444 - 2-D array (variance-covariance matrix for `T`) 445 * `error_from`: if set to `'both'` (default), returned errors take into account 446 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 447 if set to `'calib'`, only calibration uncertainties are accounted for; 448 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 449 * `return_covar`: (False by default) whether to return the full covariance matrix 450 for returned `T` or `D47` values, otherwise return standard errors for the returned 451 `T` or `D47` values instead. 452 453 **Returns (with `D47` input):** 454 455 * `T`: temperature value(s) computed from `D47` 456 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 457 458 **Returns (with `T` input):** 459 460 * `D47`: Δ47 value(s) computed from `D47` 461 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 462 463 ### Example 464 465 ````py 466 import numpy as np 467 from matplotlib import pyplot as ppl 468 from D47calib import OGLS23 as calib 469 470 X = np.linspace(1473**-2, 270**-2) 471 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 472 473 fig = ppl.figure(figsize = (5,3)) 474 ppl.subplots_adjust(bottom = .25, left = .15) 475 calib.invT_xaxis() 476 ppl.plot(X, 1000 * sD47, 'r-') 477 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 478 ppl.savefig('example_SE47.png', dpi = 100) 479 ````` 480 481 This should result in something like this: 482 483 <img src="example_SE47.png"> 484 ''' 485 486 if D47 is None and T is None: 487 raise ValueError('Either D47 or T must be specified, but both are undefined.') 488 489 if D47 is not None and T is not None: 490 raise ValueError('Either D47 or T must be specified, but not both.') 491 492 if T is not None: 493 494 D47 = self._D47_from_T(T) 495 Np = len(self.degrees) 496 N = D47.size 497 498 ### Compute covariance matrix of (*bfp, *T): 499 CM = _np.zeros((Np+N, Np+N)) 500 501 if error_from in ['calib', 'both']: 502 CM[:Np, :Np] = self.bfp_CM[:,:] 503 504 if (sT is not None) and error_from in ['sT', 'both']: 505 _sT = _np.asarray(sT) 506 if _sT.ndim == 0: 507 for k in range(N): 508 CM[Np+k, Np+k] = _sT**2 509 elif _sT.ndim == 1: 510 for k in range(N): 511 CM[Np+k, Np+k] = _sT[k]**2 512 elif _sT.ndim == 2: 513 CM[-N:, -N:] = _sT[:,:] 514 515 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 516 _T = _np.asarray(T) 517 if _T.ndim == 0: 518 _T = _np.expand_dims(_T, 0) 519 J = _np.zeros((N, Np+N)) 520 521 if (sT is not None) and error_from in ['sT', 'both']: 522 for k in range(N): 523 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 524 525 if error_from in ['calib', 'both']: 526 527 for k in range(Np): 528 529 p1 = {_: self.bfp[_] for _ in self.bfp} 530 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 531 532 p2 = {_: self.bfp[_] for _ in self.bfp} 533 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 534 535 J[:, k] = (self.model_fun(p1, (_T+273.15)**-1) - self.model_fun(p2, (_T+273.15)**-1)) / (0.002 * self.bfp_CM[k,k]**.5) 536 537 ### Error propagation: 538 CM_D47 = J @ CM @ J.T 539 540 if return_covar: 541 return D47, CM_D47 542 else: 543 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 544 545 if D47 is not None: 546 547 T = self._T_from_D47(D47) 548 Np = len(self.degrees) 549 N = T.size 550 551 ### Compute covariance matrix of (*bfp, *T): 552 CM = _np.zeros((Np+N, Np+N)) 553 554 if error_from in ['calib', 'both']: 555 CM[:Np, :Np] = self.bfp_CM[:,:] 556 557 if (sD47 is not None) and error_from in ['sD47', 'both']: 558 _sD47 = _np.asarray(sD47) 559 if _sD47.ndim == 0: 560 for k in range(N): 561 CM[Np+k, Np+k] = _sD47**2 562 elif _sD47.ndim == 1: 563 for k in range(N): 564 CM[Np+k, Np+k] = _sD47[k]**2 565 elif _sD47.ndim == 2: 566 CM[-N:, -N:] = _sD47[:,:] 567 568 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 569 _D47 = _np.asarray(D47) 570 if _D47.ndim == 0: 571 _D47 = _np.expand_dims(_D47, 0) 572 J = _np.zeros((N, Np+N)) 573 if (sD47 is not None) and error_from in ['sD47', 'both']: 574 for k in range(N): 575 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 576 if error_from in ['calib', 'both']: 577 578 xi = _np.linspace(0,200**-1,1001)[1:] 579 for k in range(Np): 580 581 if self.bfp_CM[k,k]: 582 _epsilon_ = self.bfp_CM[k,k]**.5 583 else: 584 _epsilon_ = 1e-6 585 586 p1 = {_: self.bfp[_] for _ in self.bfp} 587 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 588 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 589 590 p2 = {_: self.bfp[_] for _ in self.bfp} 591 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 592 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 593 594 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 595 596 ### Error propagation: 597 CM_T = J @ CM @ J.T 598 599 if return_covar: 600 return T, CM_T 601 else: 602 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5
When D47 is input, computes corresponding T value(s).
D47 input may be specified as a scalar, or as a 1-D array.
T output will then have the same type and size as D47.
When T is input, computes corresponding Δ47 value(s).
T input may be specified as a scalar, or as a 1-D array.
D47 output will then have the same type and size as T.
Only one of either D47 or T may be specified as input.
Arguments:
D47: Δ47 value(s) to convert into temperature (floator 1-D array)sD47: Δ47 uncertainties, which may be:T: T value(s) to convert into Δ47 (floator 1-D array), in degrees CsT: T uncertainties, which may be:error_from: if set to'both'(default), returned errors take into account input uncertainties (sTorsD47) as well as calibration uncertainties; if set to'calib', only calibration uncertainties are accounted for; if set to'sT'or'sD47', calibration uncertainties are ignored.return_covar: (False by default) whether to return the full covariance matrix for returnedTorD47values, otherwise return standard errors for the returnedTorD47values instead.
Returns (with D47 input):
T: temperature value(s) computed fromD47sT: uncertainties onTvalue(s), whether as standard error(s) or covariance matrix
Returns (with T input):
D47: Δ47 value(s) computed fromD47sD47: uncertainties onD47value(s), whether as standard error(s) or covariance matrix
Example
import numpy as np
from matplotlib import pyplot as ppl
from D47calib import OGLS23 as calib
X = np.linspace(1473**-2, 270**-2)
D47, sD47 = calib.T47(T = X**-0.5 - 273.15)
fig = ppl.figure(figsize = (5,3))
ppl.subplots_adjust(bottom = .25, left = .15)
calib.invT_xaxis()
ppl.plot(X, 1000 * sD47, 'r-')
ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)')
ppl.savefig('example_SE47.png', dpi = 100)
This should result in something like this:
def
plot_T47_errors( self, calibname=None, rD47=0.01, Nr=[2, 4, 8, 12, 20], Tmin=0, Tmax=120, colors=[(1, 0, 0), (1, 0.5, 0), (0.25, 0.75, 0), (0, 0.5, 1), (0.5, 0.5, 0.5)], yscale='lin'):
605 def plot_T47_errors( 606 self, 607 calibname = None, 608 rD47 = 0.010, 609 Nr = [2,4,8,12,20], 610 Tmin = 0, 611 Tmax = 120, 612 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 613 yscale = 'lin', 614 ): 615 """ 616 Plot SE of T reconstructed using the calibration as a function of T for various 617 combinations of analytical precision and number of analytical replicates. 618 619 **Arguments** 620 621 + **calibname**: 622 Which calibration name to display. By default, use `label` attribute. 623 + **rD47**: 624 Analytical precision of a single analysis. 625 + **Nr**: 626 A list of lines to plot, each corresponding to a given number of replicates. 627 + **Tmin**: 628 Minimum T to plot. 629 + **Tmax**: 630 Maximum T to plot. 631 + **colors**: 632 A list of colors to distinguish the plotted lines. 633 + **yscale**: 634 + If `'lin'`, the Y axis uses a linear scale. 635 + If `'log'`, the Y axis uses a logarithmic scale. 636 637 **Example** 638 639 ````py 640 from matplotlib import pyplot as ppl 641 from D47calib import devils_laghetto_2023 as calib 642 643 fig = ppl.figure(figsize = (3.5,4)) 644 ppl.subplots_adjust(bottom = .2, left = .15) 645 calib.plot_T47_errors( 646 calibname = 'Devils Laghetto calibration', 647 Nr = [1,2,4,16], 648 Tmin =0, 649 Tmax = 40, 650 ) 651 ppl.savefig('example_SE_T.png', dpi = 100) 652 ```` 653 654 This should result in something like this: 655 656 <img src="example_SE_T.png"> 657 """ 658 659 if calibname is None: 660 calibname = self.label 661 662 Nr = _np.array(Nr) 663 if len(colors) < Nr.size: 664 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 665 from colorsys import hsv_to_rgb 666 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 667 colors = [hsv_to_rgb(*x) for x in hsv] 668 669 Ti = _np.linspace(Tmin, Tmax) 670 D47i, _ = self.T47(T = Ti) 671 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 672 673 ymax, ymin = 0, 1e6 674 for N,c in zip(Nr, colors): 675 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 676 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 677 ymin = min(ymin, min(sT)) 678 ymax = max(ymax, max(sT)) 679 680 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 681 682 _ppl.legend(fontsize=9) 683 _ppl.xlabel("T (°C)") 684 685 _ppl.ylabel("Standard error on reconstructed T (°C)") 686 687 # yticks([0,.5,1,1.5,2]) 688 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 689 _ppl.grid( alpha = .25) 690 if yscale == 'lin': 691 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 692 t1, t2 = self.T.min(), self.T.max() 693 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 694 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 695 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 696 elif yscale == 'log': 697 ymin /= 2 698 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 699 _ppl.yscale('log') 700 t1, t2 = self.T.min(), self.T.max() 701 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 702 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic')
Plot SE of T reconstructed using the calibration as a function of T for various combinations of analytical precision and number of analytical replicates.
Arguments
- calibname:
Which calibration name to display. By default, use
labelattribute. - rD47: Analytical precision of a single analysis.
- Nr: A list of lines to plot, each corresponding to a given number of replicates.
- Tmin: Minimum T to plot.
- Tmax: Maximum T to plot.
- colors: A list of colors to distinguish the plotted lines.
- yscale:
- If
'lin', the Y axis uses a linear scale. - If
'log', the Y axis uses a logarithmic scale.
- If
Example
from matplotlib import pyplot as ppl
from D47calib import devils_laghetto_2023 as calib
fig = ppl.figure(figsize = (3.5,4))
ppl.subplots_adjust(bottom = .2, left = .15)
calib.plot_T47_errors(
calibname = 'Devils Laghetto calibration',
Nr = [1,2,4,16],
Tmin =0,
Tmax = 40,
)
ppl.savefig('example_SE_T.png', dpi = 100)
This should result in something like this:
def
export_data( self, csvfile, sep=',', label=False, T_correl=False, D47_correl=False):
704 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 705 """ 706 Write calibration data to a csv file. 707 708 ### Parameters 709 710 + **csvfile**: 711 The filename to write data to. 712 + **sep**: 713 The separator between CSV fields. 714 + **label**: 715 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 716 + If specified as a `str`, include a `Dataset` column with that string. 717 + If specified as `False`, do not include a `Dataset` column. 718 + **T_correl**: 719 + If `True`, include correlations between all `T` values. 720 + **D47_correl**: 721 + If `True`, include correlations between all `D47` values. 722 723 ### Example 724 725 ````py 726 D47calib.huyghe_2022.export_data( 727 csvfile = 'example_export_data.csv', 728 T_correl = True, 729 D47_correl = True, 730 ) 731 ```` 732 733 This should result in something like this ([link](example_export_data.csv)): 734 735 .. include:: ../../docs/example_export_data.md 736 737 """ 738 n = len(str(self.N)) 739 740 with open(csvfile, 'w') as f: 741 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 742 743 if label: 744 f.write(f'{sep}Dataset') 745 746 if T_correl: 747 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 748 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 749 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 750 751 if D47_correl: 752 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 753 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 754 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 755 756 for k, (s, T, sT, D47, sD47) in enumerate(zip( 757 self.samples, 758 self.T, 759 _np.diag(self.sT)**.5, 760 self.D47, 761 _np.diag(self.sD47)**.5, 762 )): 763 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 764 if label: 765 if label is True: 766 f.write(f'{sep}{self.label}') 767 else: 768 f.write(f'{sep}{label}') 769 if T_correl: 770 f.write(sep.join(['']+[ 771 f'{Tcorrel[k,_]:.0f}' 772 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 773 else f'{Tcorrel[k,_]:.6f}' 774 for _ in range(self.N)])) 775 if D47_correl: 776 f.write(sep.join(['']+[ 777 f'{D47correl[k,_]:.0f}' 778 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 779 else f'{D47correl[k,_]:.6f}' 780 for _ in range(self.N)]))
Write calibration data to a csv file.
Parameters
- csvfile: The filename to write data to.
- sep: The separator between CSV fields.
- label:
- If specified as
True, include aDatasetcolumn with the calibration'slabelattribute. - If specified as a
str, include aDatasetcolumn with that string. - If specified as
False, do not include aDatasetcolumn.
- If specified as
- T_correl:
- If
True, include correlations between allTvalues.
- If
- D47_correl:
- If
True, include correlations between allD47values.
- If
Example
D47calib.huyghe_2022.export_data(
csvfile = 'example_export_data.csv',
T_correl = True,
D47_correl = True,
)
This should result in something like this (link):
| ID | Sample | T | SE_T | D47 | SE_D47 | D47correl_1 | D47correl_2 | D47correl_3 | D47correl_4 | D47correl_5 | D47correl_6 | D47correl_7 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Ad | -1.80 | 0.50 | 0.6893 | 0.0060 | 1 | 0.048049 | 0.028770 | 0.544016 | 0.093188 | 0.020880 | 0.471516 |
| 2 | BDV-S | 18.70 | 0.75 | 0.6121 | 0.0049 | 0.048049 | 1 | 0.650474 | 0.053281 | 0.011132 | 0.002494 | 0.050651 |
| 3 | BDV-W | 11.01 | 1.00 | 0.6349 | 0.0052 | 0.028770 | 0.650474 | 1 | 0.031903 | 0.006666 | 0.001494 | 0.030328 |
| 4 | PY | 13.44 | 0.06 | 0.6397 | 0.0049 | 0.544016 | 0.053281 | 0.031903 | 1 | 0.104392 | 0.023391 | 0.513257 |
| 5 | TES-S | 22.50 | 2.10 | 0.5972 | 0.0053 | 0.093188 | 0.011132 | 0.006666 | 0.104392 | 1 | 0.275629 | 0.100150 |
| 6 | TES-W | 12.23 | 1.00 | 0.6329 | 0.0102 | 0.020880 | 0.002494 | 0.001494 | 0.023391 | 0.275629 | 1 | 0.022440 |
| 7 | TW | 26.80 | 0.85 | 0.6001 | 0.0048 | 0.471516 | 0.050651 | 0.030328 | 0.513257 | 0.100150 | 0.022440 | 1 |
def
export(self, name, filename):
783 def export(self, name, filename): 784 """ 785 Save `D47calib` object as an importable file. 786 787 ### Parameters 788 789 + **name**: 790 The name of the variable to export. 791 + **filename**: 792 The filename to write to. 793 794 ### Example 795 796 ````py 797 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 798 ```` 799 800 This should result in a `bar.py` file with the following contents: 801 802 ````py 803 foo = D47calib( 804 samples = ['LGB-2', 'DVH-2'], 805 T = [7.9, 33.7], 806 D47 = [0.6485720997671647, 0.5695972909966959], 807 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 808 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 809 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 810 label = 'Slow-growing calcites from Anderson et al. (2021)', 811 color = (0, 0.5, 0), 812 degrees = [0, 2], 813 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 814 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 815 chisq = 6.421311854486162e-27, 816 Nf = 0, 817 ) 818 ```` 819 """ 820 with open(filename, 'w') as f: 821 f.write(f''' 822{name} = D47calib( 823 samples = {self.samples}, 824 T = {self.T.tolist()}, 825 D47 = {self.D47.tolist()}, 826 sT = {self.sT.tolist()}, 827 sD47 = {self.sD47.tolist()}, 828 degrees = {self.degrees}, 829 description = {repr(self.description)}, 830 name = {repr(self.name)}, 831 label = {repr(self.label)}, 832 bfp = {({k: float(self.bfp[k]) for k in self.bfp})}, 833 bfp_CM = {self.bfp_CM.tolist()}, 834 chisq = {self.chisq}, 835 cholesky_residuals = {self.cholesky_residuals.tolist()}, 836 aic = {self.aic}, 837 bic = {self.bic}, 838 ks_pvalue = {self.ks_pvalue}, 839 ) 840''')
Save D47calib object as an importable file.
Parameters
- name: The name of the variable to export.
- filename: The filename to write to.
Example
D47calib.anderson_2021_lsce.export('foo', 'bar.py')
This should result in a bar.py file with the following contents:
foo = D47calib(
samples = ['LGB-2', 'DVH-2'],
T = [7.9, 33.7],
D47 = [0.6485720997671647, 0.5695972909966959],
sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]],
sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]],
description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES',
label = 'Slow-growing calcites from Anderson et al. (2021)',
color = (0, 0.5, 0),
degrees = [0, 2],
bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721},
bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]],
chisq = 6.421311854486162e-27,
Nf = 0,
)
def
combine_D47calibs(calibs, degrees=[0, 2], same_T=[], exclude_samples=[]):
842def combine_D47calibs(calibs, degrees = [0,2], same_T = [], exclude_samples = []): 843 ''' 844 Combine data from several `D47calib` instances. 845 846 ### Parameters 847 848 + **calibs**: 849 A list of `D47calib` instances 850 + **degrees**: 851 The polynomial degrees of the combined regression. 852 + **same_T**: 853 Use this `list` to specify when samples from different calibrations are known/postulated 854 to have formed at the same temperature (e.g. `DVH-2` and `DHC2-8` from the `fiebig_2021` 855 and `anderson_2021_lsce` data sets). Each element of `same_T` is a `list` with the names 856 of two or more samples formed at the same temperature. 857 + **exclude_samples**: Use this `list` to specify the names of samples to exclude from 858 the combined calibration. 859 860 For example, the `OGLS23` calibration is computed with: 861 862 `same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]` 863 864 Note that when samples from different calibrations have the same name, 865 it is not necessary to explicitly list them in `same_T`. 866 867 Also note that the regression will fail if samples listed together in `same_T` 868 actually have different `T` values specified in the original calibrations. 869 870 ### Example 871 872 The `devils_laghetto_2023` calibration is computed using the following code: 873 874 ````py 875 K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']] 876 877 fiebig_temp = D47calib( 878 samples = [fiebig_2021.samples[_] for _ in K], 879 T = fiebig_2021.T[K], 880 D47 = fiebig_2021.D47[K], 881 sT = fiebig_2021.sT[K,:][:,K], 882 sD47 = fiebig_2021.sD47[K,:][:,K], 883 ) 884 885 devils_laghetto_2023 = combine_D47calibs( 886 calibs = [anderson_2021_lsce, fiebig_temp], 887 degrees = [0,2], 888 same_T = [{'DVH-2', 'DHC2-8'}], 889 ) 890 ```` 891 ''' 892 893 samples = [s for c in calibs for s in c.samples] 894 T = [t for c in calibs for t in c.T] 895 D47 = [x for c in calibs for x in c.D47] 896 sD47 = _block_diag(*[c.sD47 for c in calibs]) 897 sT = _block_diag(*[c.sT for c in calibs]) 898 899 for i in range(len(samples)): 900 for j in range(len(samples)): 901 if i != j: 902 if (samples[i] == samples[j] or 903 any([samples[i] in _ and samples[j] in _ for _ in same_T])): 904 905 sT[i,j] = (sT[i,i] * sT[j,j])**.5 906 907 k = [_ for _, s in enumerate(samples) if s not in exclude_samples] 908 909 calib = D47calib( 910 samples = [samples[_] for _ in k], 911 T = [T[_] for _ in k], 912 D47 = [D47[_] for _ in k], 913 sT = sT[k,:][:,k], 914 sD47 = sD47[k,:][:,k], 915 degrees = degrees, 916 ) 917 918 return calib
Combine data from several D47calib instances.
Parameters
- calibs:
A list of
D47calibinstances - degrees: The polynomial degrees of the combined regression.
- same_T:
Use this
listto specify when samples from different calibrations are known/postulated to have formed at the same temperature (e.g.DVH-2andDHC2-8from thefiebig_2021andanderson_2021_lscedata sets). Each element ofsame_Tis alistwith the names of two or more samples formed at the same temperature. - exclude_samples: Use this
listto specify the names of samples to exclude from the combined calibration.
For example, the OGLS23 calibration is computed with:
same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]
Note that when samples from different calibrations have the same name,
it is not necessary to explicitly list them in same_T.
Also note that the regression will fail if samples listed together in same_T
actually have different T values specified in the original calibrations.
Example
The devils_laghetto_2023 calibration is computed using the following code:
K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']]
fiebig_temp = D47calib(
samples = [fiebig_2021.samples[_] for _ in K],
T = fiebig_2021.T[K],
D47 = fiebig_2021.D47[K],
sT = fiebig_2021.sT[K,:][:,K],
sD47 = fiebig_2021.sD47[K,:][:,K],
)
devils_laghetto_2023 = combine_D47calibs(
calibs = [anderson_2021_lsce, fiebig_temp],
degrees = [0,2],
same_T = [{'DVH-2', 'DHC2-8'}],
)
breitenbach_2018 =
<D47calib: >
peral_2018 =
<D47calib: >
jautzy_2020 =
<D47calib: >
anderson_2021_mit =
<D47calib: >
anderson_2021_lsce =
<D47calib: anderson_2021_lsce>
fiebig_2021 =
<D47calib: >
huyghe_2022 =
<D47calib: >
devils_laghetto_2023 =
<D47calib: >
OGLS23 =
<D47calib: OGLS23>
ogls_2023 =
<D47calib: OGLS23>