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# __docformat__ = "restructuredtext" 25__date__ = '2025-09-10' 26__version__ = '1.3.3' 27 28 29import typer 30import typer.rich_utils 31import sys 32from typing_extensions import Annotated 33import ogls as _ogls 34import numpy as _np 35from scipy.linalg import block_diag as _block_diag 36from scipy.interpolate import interp1d as _interp1d 37from matplotlib import pyplot as _ppl 38 39typer.rich_utils.STYLE_HELPTEXT = '' 40 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''') 841 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 919 920from ._calibs import * 921 922def _covar2correl(C): 923 SE = _np.diag(C)**.5 924 return SE, _np.diag(SE**-1) @ C @ _np.diag(SE**-1) 925 926try: 927 app = typer.Typer( 928 add_completion = False, 929 context_settings={'help_option_names': ['-h', '--help']}, 930 rich_markup_mode = 'rich', 931 ) 932 933 @app.command() 934 def _cli( 935 input: Annotated[str, typer.Argument(help = "Specify either the path of an input file or just '-' to read input from stdin")] = '-', 936 include_samples: Annotated[str, typer.Option('--include-samples', '-u', help = 'Only include samples listed in this file')] = 'all', 937 exclude_samples: Annotated[str, typer.Option('--exclude-samples', '-x', help = 'Exclude samples listed in this file')] = 'none', 938 outfile: Annotated[str, typer.Option('--output-file', '-o', help = 'Write output to this file instead of printing to stdout')] = 'none', 939 calib: Annotated[str, typer.Option('--calib', '-c', help = 'D47 calibration function to use')] = 'OGLS23', 940 delim_in: Annotated[str, typer.Option('--delimiter-in', '-i', help = "Delimiter used in the input.")] = ',', 941 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", 942 T_precision: Annotated[int, typer.Option('--T-precision', '-p', help = 'Precision for T output')] = 2, 943 D47_precision: Annotated[int, typer.Option('--D47-precision', '-q', help = 'Precision for D47 output')] = 4, 944 correl_precision: Annotated[int, typer.Option('--correl-precision', '-r', help = 'Precision for correlation output')] = 3, 945 covar_precision: Annotated[int, typer.Option('--covar-precision', '-s', help = 'Precision for covariance output')] = 3, 946 return_covar: Annotated[bool, typer.Option('--return-covar', '-v', help = 'Output covariance matrix instead of correlation matrix')] = False, 947 ignore_correl: Annotated[bool, typer.Option('--ignore-correl', '-g', help = 'Only consider and report standard errors without correlations')] = 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 standard errors and error correlations from calibration uncertainties. 978 979For options 2-4 and 5-8, which specify standard errors or covariances for the input values, the standard errors and error correlations resulting from these input uncertainties are also computed, as well as the combined standard errors accounting for both calibration and 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 data_out[0] += [f'{outfield}_SE_from_calib'] 1195 for k in range(N): 1196 data_out[k+1] += [Y_SE_from_calib_str[k]] 1197 if not ignore_correl: 1198 if return_covar: 1199 data_out[0] += [f'{outfield}_covar_from_calib'] + ['' for _ in range(N-1)] 1200 for k in range(N): 1201 data_out[k+1] += Y_covar_from_calib_str[k] 1202 else: 1203 data_out[0] += [f'{outfield}_correl_from_calib'] + ['' for _ in range(N-1)] 1204 for k in range(N): 1205 data_out[k+1] += Y_correl_from_calib_str[k] 1206 1207 data_out[0] += [f'{outfield}_SE_from_input'] 1208 for k in range(N): 1209 data_out[k+1] += [Y_SE_from_input_str[k]] 1210 if not ignore_correl: 1211 if return_covar: 1212 data_out[0] += [f'{outfield}_covar_from_input'] + ['' for _ in range(N-1)] 1213 for k in range(N): 1214 data_out[k+1] += Y_covar_from_input_str[k] 1215 else: 1216 data_out[0] += [f'{outfield}_correl_from_input'] + ['' for _ in range(N-1)] 1217 for k in range(N): 1218 data_out[k+1] += Y_correl_from_input_str[k] 1219 1220 data_out[0] += [f'{outfield}_SE_from_both'] 1221 for k in range(N): 1222 data_out[k+1] += [Y_SE_from_both_str[k]] 1223 if not ignore_correl: 1224 if return_covar: 1225 data_out[0] += [f'{outfield}_covar_from_both'] + ['' for _ in range(N-1)] 1226 for k in range(N): 1227 data_out[k+1] += Y_covar_from_both_str[k] 1228 else: 1229 data_out[0] += [f'{outfield}_correl_from_both'] + ['' for _ in range(N-1)] 1230 for k in range(N): 1231 data_out[k+1] += Y_correl_from_both_str[k] 1232 1233 1234 ### PRINT OUTPUT TO STDOUT OR SAVE IT TO FILE 1235 if delim_out in '<>': 1236 lengths = [max([len(data_out[j][k]) for j in range(len(data_out))]) for k in range(len(data_out[0]))] 1237 1238 txt = '' 1239 for l in data_out: 1240 for k in range(len(l)): 1241 if k > 0: 1242 txt += ' ' 1243 txt += f'{l[k]:{delim_out}{lengths[k]}s}' 1244 txt += '\n' 1245 1246 txt = txt[:-1] 1247 1248 else: 1249 txt = '\n'.join([delim_out.join(l) for l in data_out]) 1250 1251 if outfile == 'none': 1252 print(txt) 1253 else: 1254 with open(outfile, 'w') as f: 1255 f.write(txt) 1256 1257 def __cli(): 1258 app() 1259 1260except NameError: 1261 pass 1262
Reference
class
D47calib(ogls.InverseTPolynomial):
42class D47calib(_ogls.InverseTPolynomial): 43 """ 44 Δ47 calibration class based on OGLS regression 45 of Δ47 as a polynomial function of inverse T. 46 """ 47 48 def __init__(self, 49 samples, T, D47, 50 sT = None, 51 sD47 = None, 52 degrees = [0,2], 53 xpower = 2, 54 name = '', 55 label = '', 56 description = '', 57 **kwargs, 58 ): 59 """ 60 ### Parameters 61 62 + **samples**: a list of N sample names. 63 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 64 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 65 + **sT**: uncertainties on `T`. If specified as: 66 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 67 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 68 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 69 + **sD47**: uncertainties on `D47`. If specified as: 70 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 71 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 72 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 73 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 74 + **name**: a human-readable, short name assigned to the calibration. 75 + **label**: a short description of the calibration, e.g., to be used in legends. 76 + **description**: a longer description, including relevant references/DOIs. 77 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 78 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 79 80 ### Notable attributes 81 82 + **N**: 83 The total number of observations (samples) in the calibration data. 84 + **samples**: 85 The list sample names. 86 + **T**: 87 1-D `ndarray` of temperatures in degrees C. 88 + **D47**: 89 1-D `ndarray` of Δ47 values in permil. 90 + **sT**: 91 2-D `ndarray` equal to the full (co)variance matrix for `T`. 92 + **D47**: 93 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 94 + **xpower**: 95 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 96 It is possible to change this behavior to use a different power of 1/T. 97 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 98 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 99 + **bfp**: 100 The best-fit parameters of the regression. 101 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 102 + **bff()**: 103 The best-fit polynomial function of inverse T, defined as: 104 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 105 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 106 107 108 ### Examples 109 110 A very simple example: 111 112 ````py 113 .. include:: ../../code_examples/D47calib_init/example.py 114 ```` 115 116 Should yield: 117 118 ```` 119 .. include:: ../../code_examples/D47calib_init/output.txt 120 ```` 121 122 """ 123 124 self.samples = samples[:] 125 self.name = name 126 self.label = label 127 self.description = description 128 self.D47 = _np.asarray(D47, dtype = 'float') 129 self.N = self.D47.size 130 131 if sD47 is None: 132 self.sD47 = _np.zeros((self.N, self.N)) 133 else: 134 self.sD47 = _np.asarray(sD47) 135 if len(self.sD47.shape) == 1: 136 self.sD47 = _np.diag(self.sD47**2) 137 elif len(self.sD47.shape) == 0: 138 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 139 140 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 141 142 if self.bfp is None: 143 self.regress() 144 145 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) 146 147 xi = _np.linspace(0,200**-1,1001) 148 self._inv_bff = _interp1d(self.bff(xi), xi) 149 150 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 151 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 152 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 153 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1 154 155 def __repr__(self): 156 return f'<D47calib: {self.name}>' 157 158 def invT_xaxis(self, 159 xlabel = None, 160 Ti = [0,20,50,100,250,1000], 161 ): 162 """ 163 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 164 but labeled in degrees Celsius. 165 166 ### Parameters 167 168 + **xlabel**: 169 Custom label for X axis (`r'$1\\,/\\,T^2$'` by default) 170 + **Ti**: 171 Specify tick locations for X axis, in degrees C. 172 173 ### Returns 174 175 + an `matplotlib.axes.Axes` instance 176 177 ### Examples 178 179 ````py 180 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 181 ```` 182 183 This should result in something like this: 184 185 <img align="center" src="example_invT_xaxis_1.png"> 186 187 It is also possible to define the X axis using a different power of 1/T 188 by first redefining the `xpower` attribute: 189 190 ````py 191 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 192 ```` 193 194 This should result in something like this: 195 196 <img align="center" src="example_invT_xaxis_2.png"> 197 """ 198 if xlabel is None: 199 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 200 _ppl.xlabel(xlabel) 201 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 202 ax = _ppl.gca() 203 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 204 ax.tick_params(which="major") 205 206 return ax 207 208 209 def plot_data(self, label = False, **kwargs): 210 """ 211 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 212 213 ### Parameters 214 215 + **label**: 216 + If `label` is a string, use this string as `label` for the underlyig 217 `matplotlib.pyplot.plot()` call. 218 + If `label = True`, use the caller's `label` attribute instead. 219 + If `label = False`, no label is specified (default behavior). 220 + **kwargs**: 221 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 222 223 ### Returns 224 225 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 226 227 ### Example 228 229 ````py 230 from matplotlib import pyplot as ppl 231 from D47calib import huyghe_2022 as calib 232 233 fig = ppl.figure(figsize = (5,3)) 234 ppl.subplots_adjust(bottom = .25, left = .15) 235 calib.invT_xaxis(Ti = [0,10,25]) 236 calib.plot_data(label = True) 237 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 238 ppl.legend() 239 ppl.savefig('example_plot_data.png', dpi = 100) 240 ````` 241 242 This should result in something like this: 243 244 <img align="center" src="example_plot_data.png"> 245 """ 246# if 'mec' not in kwargs: 247# kwargs['mec'] = self.color 248 if label is not False: 249 kwargs['label'] = self.label if label is True else label 250 return _ogls.InverseTPolynomial.plot_data(self, **kwargs) 251 252 253 def plot_error_bars(self, **kwargs): 254 """ 255 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 256 257 ### Parameters 258 259 + **kwargs**: 260 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 261 262 ### Returns 263 264 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 265 266 ### Example 267 268 ````py 269 from matplotlib import pyplot as ppl 270 from D47calib import huyghe_2022 as calib 271 272 fig = ppl.figure(figsize = (5,3)) 273 ppl.subplots_adjust(bottom = .25, left = .15) 274 calib.invT_xaxis(Ti = [0,10,25]) 275 calib.plot_error_bars(alpha = .4) 276 calib.plot_data(label = True) 277 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 278 ppl.legend() 279 ppl.savefig('example_plot_error_bars.png', dpi = 100) 280 ````` 281 282 This should result in something like this: 283 284 <img align="center" src="example_plot_error_bars.png"> 285 """ 286# if 'ecolor' not in kwargs: 287# kwargs['ecolor'] = self.color 288 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs) 289 290 291 def plot_error_ellipses(self, **kwargs): 292 """ 293 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 294 295 ### Parameters 296 297 + **kwargs**: 298 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 299 300 ### Returns 301 302 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 303 304 ### Example 305 306 ````py 307 from matplotlib import pyplot as ppl 308 from D47calib import huyghe_2022 as calib 309 310 fig = ppl.figure(figsize = (5,3)) 311 ppl.subplots_adjust(bottom = .25, left = .15) 312 calib.invT_xaxis(Ti = [0,10,25]) 313 calib.plot_error_ellipses(alpha = .4) 314 calib.plot_data(label = True) 315 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 316 ppl.legend() 317 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 318 ````` 319 320 This should result in something like this: 321 322 <img align="center" src="example_plot_error_ellipses.png"> 323 """ 324# if 'ec' not in kwargs: 325# kwargs['ec'] = self.color 326 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs) 327 328 329 def plot_bff(self, label = False, **kwargs): 330 """ 331 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 332 333 ### Parameters 334 335 + **label**: 336 + If `label` is a string, use this string as `label` for the underlyig 337 `matplotlib.pyplot.plot()` call. 338 + If `label = True`, use the caller's `label` attribute instead. 339 + If `label = False`, no label is specified (default behavior). 340 + **kwargs**: 341 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 342 343 ### Returns 344 345 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 346 347 ### Example 348 349 ````py 350 from matplotlib import pyplot as ppl 351 from D47calib import huyghe_2022 as calib 352 353 fig = ppl.figure(figsize = (5,3)) 354 ppl.subplots_adjust(bottom = .25, left = .15) 355 calib.invT_xaxis(Ti = [0,10,25]) 356 calib.plot_bff(label = True, dashes = (8,2,2,2)) 357 calib.plot_data() 358 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 359 ppl.legend() 360 ppl.savefig('example_plot_bff.png', dpi = 100) 361 ````` 362 363 This should result in something like this: 364 365 <img align="center" src="example_plot_bff.png"> 366 """ 367# if 'color' not in kwargs: 368# kwargs['color'] = self.color 369 if label is not False: 370 kwargs['label'] = self.label if label is True else label 371 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs) 372 373 374 def plot_bff_ci(self, **kwargs): 375 """ 376 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 377 378 ### Parameters 379 380 + **label**: 381 + **kwargs**: 382 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 383 384 ### Returns 385 386 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 387 388 ### Example 389 390 ````py 391 from matplotlib import pyplot as ppl 392 from D47calib import huyghe_2022 as calib 393 394 fig = ppl.figure(figsize = (5,3)) 395 ppl.subplots_adjust(bottom = .25, left = .15) 396 calib.invT_xaxis(Ti = [0,10,25]) 397 calib.plot_bff_ci(alpha = .15) 398 calib.plot_bff(label = True, dashes = (8,2,2,2)) 399 calib.plot_data() 400 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 401 ppl.legend() 402 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 403 ````` 404 405 This should result in something like this: 406 407 <img align="center" src="example_plot_bff_ci.png"> 408 """ 409# if 'color' not in kwargs: 410# kwargs['color'] = self.color 411 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs) 412 413 def T47(self, 414 D47 = None, 415 sD47 = None, 416 T=None, 417 sT = None, 418 error_from = 'both', 419 return_covar = False, 420 ): 421 ''' 422 When `D47` is input, computes corresponding T value(s). 423 `D47` input may be specified as a scalar, or as a 1-D array. 424 `T` output will then have the same type and size as `D47`. 425 426 When `T` is input, computes corresponding Δ47 value(s). 427 `T` input may be specified as a scalar, or as a 1-D array. 428 `D47` output will then have the same type and size as `T`. 429 430 Only one of either `D47` or `T` may be specified as input. 431 432 **Arguments:** 433 434 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 435 * `sD47`: Δ47 uncertainties, which may be: 436 - `None` (default) 437 - `float` or `int` (uniform standard error on `D47`) 438 - 1-D array (standard errors on `D47`) 439 - 2-D array (covariance matrix for `D47`) 440 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 441 * `sT`: T uncertainties, which may be: 442 - `None` (default) 443 - `float` or `int` (uniform standard error on `T`) 444 - 1-D array (standard errors on `T`) 445 - 2-D array (variance-covariance matrix for `T`) 446 * `error_from`: if set to `'both'` (default), returned errors take into account 447 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 448 if set to `'calib'`, only calibration uncertainties are accounted for; 449 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 450 * `return_covar`: (False by default) whether to return the full covariance matrix 451 for returned `T` or `D47` values, otherwise return standard errors for the returned 452 `T` or `D47` values instead. 453 454 **Returns (with `D47` input):** 455 456 * `T`: temperature value(s) computed from `D47` 457 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 458 459 **Returns (with `T` input):** 460 461 * `D47`: Δ47 value(s) computed from `D47` 462 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 463 464 ### Example 465 466 ````py 467 import numpy as np 468 from matplotlib import pyplot as ppl 469 from D47calib import OGLS23 as calib 470 471 X = np.linspace(1473**-2, 270**-2) 472 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 473 474 fig = ppl.figure(figsize = (5,3)) 475 ppl.subplots_adjust(bottom = .25, left = .15) 476 calib.invT_xaxis() 477 ppl.plot(X, 1000 * sD47, 'r-') 478 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 479 ppl.savefig('example_SE47.png', dpi = 100) 480 ````` 481 482 This should result in something like this: 483 484 <img src="example_SE47.png"> 485 ''' 486 487 if D47 is None and T is None: 488 raise ValueError('Either D47 or T must be specified, but both are undefined.') 489 490 if D47 is not None and T is not None: 491 raise ValueError('Either D47 or T must be specified, but not both.') 492 493 if T is not None: 494 495 D47 = self._D47_from_T(T) 496 Np = len(self.degrees) 497 N = D47.size 498 499 ### Compute covariance matrix of (*bfp, *T): 500 CM = _np.zeros((Np+N, Np+N)) 501 502 if error_from in ['calib', 'both']: 503 CM[:Np, :Np] = self.bfp_CM[:,:] 504 505 if (sT is not None) and error_from in ['sT', 'both']: 506 _sT = _np.asarray(sT) 507 if _sT.ndim == 0: 508 for k in range(N): 509 CM[Np+k, Np+k] = _sT**2 510 elif _sT.ndim == 1: 511 for k in range(N): 512 CM[Np+k, Np+k] = _sT[k]**2 513 elif _sT.ndim == 2: 514 CM[-N:, -N:] = _sT[:,:] 515 516 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 517 _T = _np.asarray(T) 518 if _T.ndim == 0: 519 _T = _np.expand_dims(_T, 0) 520 J = _np.zeros((N, Np+N)) 521 522 if (sT is not None) and error_from in ['sT', 'both']: 523 for k in range(N): 524 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 525 526 if error_from in ['calib', 'both']: 527 528 for k in range(Np): 529 530 p1 = {_: self.bfp[_] for _ in self.bfp} 531 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 532 533 p2 = {_: self.bfp[_] for _ in self.bfp} 534 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 535 536 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) 537 538 ### Error propagation: 539 CM_D47 = J @ CM @ J.T 540 541 if return_covar: 542 return D47, CM_D47 543 else: 544 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 545 546 if D47 is not None: 547 548 T = self._T_from_D47(D47) 549 Np = len(self.degrees) 550 N = T.size 551 552 ### Compute covariance matrix of (*bfp, *T): 553 CM = _np.zeros((Np+N, Np+N)) 554 555 if error_from in ['calib', 'both']: 556 CM[:Np, :Np] = self.bfp_CM[:,:] 557 558 if (sD47 is not None) and error_from in ['sD47', 'both']: 559 _sD47 = _np.asarray(sD47) 560 if _sD47.ndim == 0: 561 for k in range(N): 562 CM[Np+k, Np+k] = _sD47**2 563 elif _sD47.ndim == 1: 564 for k in range(N): 565 CM[Np+k, Np+k] = _sD47[k]**2 566 elif _sD47.ndim == 2: 567 CM[-N:, -N:] = _sD47[:,:] 568 569 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 570 _D47 = _np.asarray(D47) 571 if _D47.ndim == 0: 572 _D47 = _np.expand_dims(_D47, 0) 573 J = _np.zeros((N, Np+N)) 574 if (sD47 is not None) and error_from in ['sD47', 'both']: 575 for k in range(N): 576 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 577 if error_from in ['calib', 'both']: 578 579 xi = _np.linspace(0,200**-1,1001)[1:] 580 for k in range(Np): 581 582 if self.bfp_CM[k,k]: 583 _epsilon_ = self.bfp_CM[k,k]**.5 584 else: 585 _epsilon_ = 1e-6 586 587 p1 = {_: self.bfp[_] for _ in self.bfp} 588 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 589 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 590 591 p2 = {_: self.bfp[_] for _ in self.bfp} 592 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 593 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 594 595 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 596 597 ### Error propagation: 598 CM_T = J @ CM @ J.T 599 600 if return_covar: 601 return T, CM_T 602 else: 603 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5 604 605 606 def plot_T47_errors( 607 self, 608 calibname = None, 609 rD47 = 0.010, 610 Nr = [2,4,8,12,20], 611 Tmin = 0, 612 Tmax = 120, 613 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 614 yscale = 'lin', 615 ): 616 """ 617 Plot SE of T reconstructed using the calibration as a function of T for various 618 combinations of analytical precision and number of analytical replicates. 619 620 **Arguments** 621 622 + **calibname**: 623 Which calibration name to display. By default, use `label` attribute. 624 + **rD47**: 625 Analytical precision of a single analysis. 626 + **Nr**: 627 A list of lines to plot, each corresponding to a given number of replicates. 628 + **Tmin**: 629 Minimum T to plot. 630 + **Tmax**: 631 Maximum T to plot. 632 + **colors**: 633 A list of colors to distinguish the plotted lines. 634 + **yscale**: 635 + If `'lin'`, the Y axis uses a linear scale. 636 + If `'log'`, the Y axis uses a logarithmic scale. 637 638 **Example** 639 640 ````py 641 from matplotlib import pyplot as ppl 642 from D47calib import devils_laghetto_2023 as calib 643 644 fig = ppl.figure(figsize = (3.5,4)) 645 ppl.subplots_adjust(bottom = .2, left = .15) 646 calib.plot_T47_errors( 647 calibname = 'Devils Laghetto calibration', 648 Nr = [1,2,4,16], 649 Tmin =0, 650 Tmax = 40, 651 ) 652 ppl.savefig('example_SE_T.png', dpi = 100) 653 ```` 654 655 This should result in something like this: 656 657 <img src="example_SE_T.png"> 658 """ 659 660 if calibname is None: 661 calibname = self.label 662 663 Nr = _np.array(Nr) 664 if len(colors) < Nr.size: 665 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 666 from colorsys import hsv_to_rgb 667 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 668 colors = [hsv_to_rgb(*x) for x in hsv] 669 670 Ti = _np.linspace(Tmin, Tmax) 671 D47i, _ = self.T47(T = Ti) 672 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 673 674 ymax, ymin = 0, 1e6 675 for N,c in zip(Nr, colors): 676 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 677 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 678 ymin = min(ymin, min(sT)) 679 ymax = max(ymax, max(sT)) 680 681 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 682 683 _ppl.legend(fontsize=9) 684 _ppl.xlabel("T (°C)") 685 686 _ppl.ylabel("Standard error on reconstructed T (°C)") 687 688 # yticks([0,.5,1,1.5,2]) 689 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 690 _ppl.grid( alpha = .25) 691 if yscale == 'lin': 692 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 693 t1, t2 = self.T.min(), self.T.max() 694 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 695 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 696 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 697 elif yscale == 'log': 698 ymin /= 2 699 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 700 _ppl.yscale('log') 701 t1, t2 = self.T.min(), self.T.max() 702 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 703 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 704 705 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 706 """ 707 Write calibration data to a csv file. 708 709 ### Parameters 710 711 + **csvfile**: 712 The filename to write data to. 713 + **sep**: 714 The separator between CSV fields. 715 + **label**: 716 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 717 + If specified as a `str`, include a `Dataset` column with that string. 718 + If specified as `False`, do not include a `Dataset` column. 719 + **T_correl**: 720 + If `True`, include correlations between all `T` values. 721 + **D47_correl**: 722 + If `True`, include correlations between all `D47` values. 723 724 ### Example 725 726 ````py 727 D47calib.huyghe_2022.export_data( 728 csvfile = 'example_export_data.csv', 729 T_correl = True, 730 D47_correl = True, 731 ) 732 ```` 733 734 This should result in something like this ([link](example_export_data.csv)): 735 736 .. include:: ../../docs/example_export_data.md 737 738 """ 739 n = len(str(self.N)) 740 741 with open(csvfile, 'w') as f: 742 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 743 744 if label: 745 f.write(f'{sep}Dataset') 746 747 if T_correl: 748 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 749 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 750 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 751 752 if D47_correl: 753 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 754 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 755 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 756 757 for k, (s, T, sT, D47, sD47) in enumerate(zip( 758 self.samples, 759 self.T, 760 _np.diag(self.sT)**.5, 761 self.D47, 762 _np.diag(self.sD47)**.5, 763 )): 764 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 765 if label: 766 if label is True: 767 f.write(f'{sep}{self.label}') 768 else: 769 f.write(f'{sep}{label}') 770 if T_correl: 771 f.write(sep.join(['']+[ 772 f'{Tcorrel[k,_]:.0f}' 773 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 774 else f'{Tcorrel[k,_]:.6f}' 775 for _ in range(self.N)])) 776 if D47_correl: 777 f.write(sep.join(['']+[ 778 f'{D47correl[k,_]:.0f}' 779 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 780 else f'{D47correl[k,_]:.6f}' 781 for _ in range(self.N)])) 782 783 784 def export(self, name, filename): 785 """ 786 Save `D47calib` object as an importable file. 787 788 ### Parameters 789 790 + **name**: 791 The name of the variable to export. 792 + **filename**: 793 The filename to write to. 794 795 ### Example 796 797 ````py 798 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 799 ```` 800 801 This should result in a `bar.py` file with the following contents: 802 803 ````py 804 foo = D47calib( 805 samples = ['LGB-2', 'DVH-2'], 806 T = [7.9, 33.7], 807 D47 = [0.6485720997671647, 0.5695972909966959], 808 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 809 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 810 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 811 label = 'Slow-growing calcites from Anderson et al. (2021)', 812 color = (0, 0.5, 0), 813 degrees = [0, 2], 814 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 815 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 816 chisq = 6.421311854486162e-27, 817 Nf = 0, 818 ) 819 ```` 820 """ 821 with open(filename, 'w') as f: 822 f.write(f''' 823{name} = D47calib( 824 samples = {self.samples}, 825 T = {self.T.tolist()}, 826 D47 = {self.D47.tolist()}, 827 sT = {self.sT.tolist()}, 828 sD47 = {self.sD47.tolist()}, 829 degrees = {self.degrees}, 830 description = {repr(self.description)}, 831 name = {repr(self.name)}, 832 label = {repr(self.label)}, 833 bfp = {({k: float(self.bfp[k]) for k in self.bfp})}, 834 bfp_CM = {self.bfp_CM.tolist()}, 835 chisq = {self.chisq}, 836 cholesky_residuals = {self.cholesky_residuals.tolist()}, 837 aic = {self.aic}, 838 bic = {self.bic}, 839 ks_pvalue = {self.ks_pvalue}, 840 ) 841''')
Δ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)
48 def __init__(self, 49 samples, T, D47, 50 sT = None, 51 sD47 = None, 52 degrees = [0,2], 53 xpower = 2, 54 name = '', 55 label = '', 56 description = '', 57 **kwargs, 58 ): 59 """ 60 ### Parameters 61 62 + **samples**: a list of N sample names. 63 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 64 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 65 + **sT**: uncertainties on `T`. If specified as: 66 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 67 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 68 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 69 + **sD47**: uncertainties on `D47`. If specified as: 70 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 71 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 72 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 73 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 74 + **name**: a human-readable, short name assigned to the calibration. 75 + **label**: a short description of the calibration, e.g., to be used in legends. 76 + **description**: a longer description, including relevant references/DOIs. 77 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 78 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 79 80 ### Notable attributes 81 82 + **N**: 83 The total number of observations (samples) in the calibration data. 84 + **samples**: 85 The list sample names. 86 + **T**: 87 1-D `ndarray` of temperatures in degrees C. 88 + **D47**: 89 1-D `ndarray` of Δ47 values in permil. 90 + **sT**: 91 2-D `ndarray` equal to the full (co)variance matrix for `T`. 92 + **D47**: 93 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 94 + **xpower**: 95 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 96 It is possible to change this behavior to use a different power of 1/T. 97 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 98 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 99 + **bfp**: 100 The best-fit parameters of the regression. 101 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 102 + **bff()**: 103 The best-fit polynomial function of inverse T, defined as: 104 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 105 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 106 107 108 ### Examples 109 110 A very simple example: 111 112 ````py 113 .. include:: ../../code_examples/D47calib_init/example.py 114 ```` 115 116 Should yield: 117 118 ```` 119 .. include:: ../../code_examples/D47calib_init/output.txt 120 ```` 121 122 """ 123 124 self.samples = samples[:] 125 self.name = name 126 self.label = label 127 self.description = description 128 self.D47 = _np.asarray(D47, dtype = 'float') 129 self.N = self.D47.size 130 131 if sD47 is None: 132 self.sD47 = _np.zeros((self.N, self.N)) 133 else: 134 self.sD47 = _np.asarray(sD47) 135 if len(self.sD47.shape) == 1: 136 self.sD47 = _np.diag(self.sD47**2) 137 elif len(self.sD47.shape) == 0: 138 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 139 140 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 141 142 if self.bfp is None: 143 self.regress() 144 145 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) 146 147 xi = _np.linspace(0,200**-1,1001) 148 self._inv_bff = _interp1d(self.bff(xi), xi) 149 150 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 151 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 152 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 153 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]):
158 def invT_xaxis(self, 159 xlabel = None, 160 Ti = [0,20,50,100,250,1000], 161 ): 162 """ 163 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 164 but labeled in degrees Celsius. 165 166 ### Parameters 167 168 + **xlabel**: 169 Custom label for X axis (`r'$1\\,/\\,T^2$'` by default) 170 + **Ti**: 171 Specify tick locations for X axis, in degrees C. 172 173 ### Returns 174 175 + an `matplotlib.axes.Axes` instance 176 177 ### Examples 178 179 ````py 180 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 181 ```` 182 183 This should result in something like this: 184 185 <img align="center" src="example_invT_xaxis_1.png"> 186 187 It is also possible to define the X axis using a different power of 1/T 188 by first redefining the `xpower` attribute: 189 190 ````py 191 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 192 ```` 193 194 This should result in something like this: 195 196 <img align="center" src="example_invT_xaxis_2.png"> 197 """ 198 if xlabel is None: 199 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 200 _ppl.xlabel(xlabel) 201 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 202 ax = _ppl.gca() 203 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 204 ax.tick_params(which="major") 205 206 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):
209 def plot_data(self, label = False, **kwargs): 210 """ 211 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 212 213 ### Parameters 214 215 + **label**: 216 + If `label` is a string, use this string as `label` for the underlyig 217 `matplotlib.pyplot.plot()` call. 218 + If `label = True`, use the caller's `label` attribute instead. 219 + If `label = False`, no label is specified (default behavior). 220 + **kwargs**: 221 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 222 223 ### Returns 224 225 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 226 227 ### Example 228 229 ````py 230 from matplotlib import pyplot as ppl 231 from D47calib import huyghe_2022 as calib 232 233 fig = ppl.figure(figsize = (5,3)) 234 ppl.subplots_adjust(bottom = .25, left = .15) 235 calib.invT_xaxis(Ti = [0,10,25]) 236 calib.plot_data(label = True) 237 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 238 ppl.legend() 239 ppl.savefig('example_plot_data.png', dpi = 100) 240 ````` 241 242 This should result in something like this: 243 244 <img align="center" src="example_plot_data.png"> 245 """ 246# if 'mec' not in kwargs: 247# kwargs['mec'] = self.color 248 if label is not False: 249 kwargs['label'] = self.label if label is True else label 250 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):
253 def plot_error_bars(self, **kwargs): 254 """ 255 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 256 257 ### Parameters 258 259 + **kwargs**: 260 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 261 262 ### Returns 263 264 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 265 266 ### Example 267 268 ````py 269 from matplotlib import pyplot as ppl 270 from D47calib import huyghe_2022 as calib 271 272 fig = ppl.figure(figsize = (5,3)) 273 ppl.subplots_adjust(bottom = .25, left = .15) 274 calib.invT_xaxis(Ti = [0,10,25]) 275 calib.plot_error_bars(alpha = .4) 276 calib.plot_data(label = True) 277 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 278 ppl.legend() 279 ppl.savefig('example_plot_error_bars.png', dpi = 100) 280 ````` 281 282 This should result in something like this: 283 284 <img align="center" src="example_plot_error_bars.png"> 285 """ 286# if 'ecolor' not in kwargs: 287# kwargs['ecolor'] = self.color 288 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):
291 def plot_error_ellipses(self, **kwargs): 292 """ 293 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 294 295 ### Parameters 296 297 + **kwargs**: 298 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 299 300 ### Returns 301 302 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 303 304 ### Example 305 306 ````py 307 from matplotlib import pyplot as ppl 308 from D47calib import huyghe_2022 as calib 309 310 fig = ppl.figure(figsize = (5,3)) 311 ppl.subplots_adjust(bottom = .25, left = .15) 312 calib.invT_xaxis(Ti = [0,10,25]) 313 calib.plot_error_ellipses(alpha = .4) 314 calib.plot_data(label = True) 315 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 316 ppl.legend() 317 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 318 ````` 319 320 This should result in something like this: 321 322 <img align="center" src="example_plot_error_ellipses.png"> 323 """ 324# if 'ec' not in kwargs: 325# kwargs['ec'] = self.color 326 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):
329 def plot_bff(self, label = False, **kwargs): 330 """ 331 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 332 333 ### Parameters 334 335 + **label**: 336 + If `label` is a string, use this string as `label` for the underlyig 337 `matplotlib.pyplot.plot()` call. 338 + If `label = True`, use the caller's `label` attribute instead. 339 + If `label = False`, no label is specified (default behavior). 340 + **kwargs**: 341 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 342 343 ### Returns 344 345 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 346 347 ### Example 348 349 ````py 350 from matplotlib import pyplot as ppl 351 from D47calib import huyghe_2022 as calib 352 353 fig = ppl.figure(figsize = (5,3)) 354 ppl.subplots_adjust(bottom = .25, left = .15) 355 calib.invT_xaxis(Ti = [0,10,25]) 356 calib.plot_bff(label = True, dashes = (8,2,2,2)) 357 calib.plot_data() 358 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 359 ppl.legend() 360 ppl.savefig('example_plot_bff.png', dpi = 100) 361 ````` 362 363 This should result in something like this: 364 365 <img align="center" src="example_plot_bff.png"> 366 """ 367# if 'color' not in kwargs: 368# kwargs['color'] = self.color 369 if label is not False: 370 kwargs['label'] = self.label if label is True else label 371 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):
374 def plot_bff_ci(self, **kwargs): 375 """ 376 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 377 378 ### Parameters 379 380 + **label**: 381 + **kwargs**: 382 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 383 384 ### Returns 385 386 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 387 388 ### Example 389 390 ````py 391 from matplotlib import pyplot as ppl 392 from D47calib import huyghe_2022 as calib 393 394 fig = ppl.figure(figsize = (5,3)) 395 ppl.subplots_adjust(bottom = .25, left = .15) 396 calib.invT_xaxis(Ti = [0,10,25]) 397 calib.plot_bff_ci(alpha = .15) 398 calib.plot_bff(label = True, dashes = (8,2,2,2)) 399 calib.plot_data() 400 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 401 ppl.legend() 402 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 403 ````` 404 405 This should result in something like this: 406 407 <img align="center" src="example_plot_bff_ci.png"> 408 """ 409# if 'color' not in kwargs: 410# kwargs['color'] = self.color 411 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):
413 def T47(self, 414 D47 = None, 415 sD47 = None, 416 T=None, 417 sT = None, 418 error_from = 'both', 419 return_covar = False, 420 ): 421 ''' 422 When `D47` is input, computes corresponding T value(s). 423 `D47` input may be specified as a scalar, or as a 1-D array. 424 `T` output will then have the same type and size as `D47`. 425 426 When `T` is input, computes corresponding Δ47 value(s). 427 `T` input may be specified as a scalar, or as a 1-D array. 428 `D47` output will then have the same type and size as `T`. 429 430 Only one of either `D47` or `T` may be specified as input. 431 432 **Arguments:** 433 434 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 435 * `sD47`: Δ47 uncertainties, which may be: 436 - `None` (default) 437 - `float` or `int` (uniform standard error on `D47`) 438 - 1-D array (standard errors on `D47`) 439 - 2-D array (covariance matrix for `D47`) 440 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 441 * `sT`: T uncertainties, which may be: 442 - `None` (default) 443 - `float` or `int` (uniform standard error on `T`) 444 - 1-D array (standard errors on `T`) 445 - 2-D array (variance-covariance matrix for `T`) 446 * `error_from`: if set to `'both'` (default), returned errors take into account 447 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 448 if set to `'calib'`, only calibration uncertainties are accounted for; 449 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 450 * `return_covar`: (False by default) whether to return the full covariance matrix 451 for returned `T` or `D47` values, otherwise return standard errors for the returned 452 `T` or `D47` values instead. 453 454 **Returns (with `D47` input):** 455 456 * `T`: temperature value(s) computed from `D47` 457 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 458 459 **Returns (with `T` input):** 460 461 * `D47`: Δ47 value(s) computed from `D47` 462 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 463 464 ### Example 465 466 ````py 467 import numpy as np 468 from matplotlib import pyplot as ppl 469 from D47calib import OGLS23 as calib 470 471 X = np.linspace(1473**-2, 270**-2) 472 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 473 474 fig = ppl.figure(figsize = (5,3)) 475 ppl.subplots_adjust(bottom = .25, left = .15) 476 calib.invT_xaxis() 477 ppl.plot(X, 1000 * sD47, 'r-') 478 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 479 ppl.savefig('example_SE47.png', dpi = 100) 480 ````` 481 482 This should result in something like this: 483 484 <img src="example_SE47.png"> 485 ''' 486 487 if D47 is None and T is None: 488 raise ValueError('Either D47 or T must be specified, but both are undefined.') 489 490 if D47 is not None and T is not None: 491 raise ValueError('Either D47 or T must be specified, but not both.') 492 493 if T is not None: 494 495 D47 = self._D47_from_T(T) 496 Np = len(self.degrees) 497 N = D47.size 498 499 ### Compute covariance matrix of (*bfp, *T): 500 CM = _np.zeros((Np+N, Np+N)) 501 502 if error_from in ['calib', 'both']: 503 CM[:Np, :Np] = self.bfp_CM[:,:] 504 505 if (sT is not None) and error_from in ['sT', 'both']: 506 _sT = _np.asarray(sT) 507 if _sT.ndim == 0: 508 for k in range(N): 509 CM[Np+k, Np+k] = _sT**2 510 elif _sT.ndim == 1: 511 for k in range(N): 512 CM[Np+k, Np+k] = _sT[k]**2 513 elif _sT.ndim == 2: 514 CM[-N:, -N:] = _sT[:,:] 515 516 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 517 _T = _np.asarray(T) 518 if _T.ndim == 0: 519 _T = _np.expand_dims(_T, 0) 520 J = _np.zeros((N, Np+N)) 521 522 if (sT is not None) and error_from in ['sT', 'both']: 523 for k in range(N): 524 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 525 526 if error_from in ['calib', 'both']: 527 528 for k in range(Np): 529 530 p1 = {_: self.bfp[_] for _ in self.bfp} 531 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 532 533 p2 = {_: self.bfp[_] for _ in self.bfp} 534 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 535 536 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) 537 538 ### Error propagation: 539 CM_D47 = J @ CM @ J.T 540 541 if return_covar: 542 return D47, CM_D47 543 else: 544 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 545 546 if D47 is not None: 547 548 T = self._T_from_D47(D47) 549 Np = len(self.degrees) 550 N = T.size 551 552 ### Compute covariance matrix of (*bfp, *T): 553 CM = _np.zeros((Np+N, Np+N)) 554 555 if error_from in ['calib', 'both']: 556 CM[:Np, :Np] = self.bfp_CM[:,:] 557 558 if (sD47 is not None) and error_from in ['sD47', 'both']: 559 _sD47 = _np.asarray(sD47) 560 if _sD47.ndim == 0: 561 for k in range(N): 562 CM[Np+k, Np+k] = _sD47**2 563 elif _sD47.ndim == 1: 564 for k in range(N): 565 CM[Np+k, Np+k] = _sD47[k]**2 566 elif _sD47.ndim == 2: 567 CM[-N:, -N:] = _sD47[:,:] 568 569 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 570 _D47 = _np.asarray(D47) 571 if _D47.ndim == 0: 572 _D47 = _np.expand_dims(_D47, 0) 573 J = _np.zeros((N, Np+N)) 574 if (sD47 is not None) and error_from in ['sD47', 'both']: 575 for k in range(N): 576 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 577 if error_from in ['calib', 'both']: 578 579 xi = _np.linspace(0,200**-1,1001)[1:] 580 for k in range(Np): 581 582 if self.bfp_CM[k,k]: 583 _epsilon_ = self.bfp_CM[k,k]**.5 584 else: 585 _epsilon_ = 1e-6 586 587 p1 = {_: self.bfp[_] for _ in self.bfp} 588 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 589 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 590 591 p2 = {_: self.bfp[_] for _ in self.bfp} 592 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 593 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 594 595 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 596 597 ### Error propagation: 598 CM_T = J @ CM @ J.T 599 600 if return_covar: 601 return T, CM_T 602 else: 603 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'):
606 def plot_T47_errors( 607 self, 608 calibname = None, 609 rD47 = 0.010, 610 Nr = [2,4,8,12,20], 611 Tmin = 0, 612 Tmax = 120, 613 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 614 yscale = 'lin', 615 ): 616 """ 617 Plot SE of T reconstructed using the calibration as a function of T for various 618 combinations of analytical precision and number of analytical replicates. 619 620 **Arguments** 621 622 + **calibname**: 623 Which calibration name to display. By default, use `label` attribute. 624 + **rD47**: 625 Analytical precision of a single analysis. 626 + **Nr**: 627 A list of lines to plot, each corresponding to a given number of replicates. 628 + **Tmin**: 629 Minimum T to plot. 630 + **Tmax**: 631 Maximum T to plot. 632 + **colors**: 633 A list of colors to distinguish the plotted lines. 634 + **yscale**: 635 + If `'lin'`, the Y axis uses a linear scale. 636 + If `'log'`, the Y axis uses a logarithmic scale. 637 638 **Example** 639 640 ````py 641 from matplotlib import pyplot as ppl 642 from D47calib import devils_laghetto_2023 as calib 643 644 fig = ppl.figure(figsize = (3.5,4)) 645 ppl.subplots_adjust(bottom = .2, left = .15) 646 calib.plot_T47_errors( 647 calibname = 'Devils Laghetto calibration', 648 Nr = [1,2,4,16], 649 Tmin =0, 650 Tmax = 40, 651 ) 652 ppl.savefig('example_SE_T.png', dpi = 100) 653 ```` 654 655 This should result in something like this: 656 657 <img src="example_SE_T.png"> 658 """ 659 660 if calibname is None: 661 calibname = self.label 662 663 Nr = _np.array(Nr) 664 if len(colors) < Nr.size: 665 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 666 from colorsys import hsv_to_rgb 667 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 668 colors = [hsv_to_rgb(*x) for x in hsv] 669 670 Ti = _np.linspace(Tmin, Tmax) 671 D47i, _ = self.T47(T = Ti) 672 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 673 674 ymax, ymin = 0, 1e6 675 for N,c in zip(Nr, colors): 676 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 677 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 678 ymin = min(ymin, min(sT)) 679 ymax = max(ymax, max(sT)) 680 681 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 682 683 _ppl.legend(fontsize=9) 684 _ppl.xlabel("T (°C)") 685 686 _ppl.ylabel("Standard error on reconstructed T (°C)") 687 688 # yticks([0,.5,1,1.5,2]) 689 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 690 _ppl.grid( alpha = .25) 691 if yscale == 'lin': 692 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 693 t1, t2 = self.T.min(), self.T.max() 694 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 695 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 696 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 697 elif yscale == 'log': 698 ymin /= 2 699 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 700 _ppl.yscale('log') 701 t1, t2 = self.T.min(), self.T.max() 702 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 703 _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):
705 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 706 """ 707 Write calibration data to a csv file. 708 709 ### Parameters 710 711 + **csvfile**: 712 The filename to write data to. 713 + **sep**: 714 The separator between CSV fields. 715 + **label**: 716 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 717 + If specified as a `str`, include a `Dataset` column with that string. 718 + If specified as `False`, do not include a `Dataset` column. 719 + **T_correl**: 720 + If `True`, include correlations between all `T` values. 721 + **D47_correl**: 722 + If `True`, include correlations between all `D47` values. 723 724 ### Example 725 726 ````py 727 D47calib.huyghe_2022.export_data( 728 csvfile = 'example_export_data.csv', 729 T_correl = True, 730 D47_correl = True, 731 ) 732 ```` 733 734 This should result in something like this ([link](example_export_data.csv)): 735 736 .. include:: ../../docs/example_export_data.md 737 738 """ 739 n = len(str(self.N)) 740 741 with open(csvfile, 'w') as f: 742 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 743 744 if label: 745 f.write(f'{sep}Dataset') 746 747 if T_correl: 748 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 749 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 750 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 751 752 if D47_correl: 753 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 754 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 755 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 756 757 for k, (s, T, sT, D47, sD47) in enumerate(zip( 758 self.samples, 759 self.T, 760 _np.diag(self.sT)**.5, 761 self.D47, 762 _np.diag(self.sD47)**.5, 763 )): 764 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 765 if label: 766 if label is True: 767 f.write(f'{sep}{self.label}') 768 else: 769 f.write(f'{sep}{label}') 770 if T_correl: 771 f.write(sep.join(['']+[ 772 f'{Tcorrel[k,_]:.0f}' 773 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 774 else f'{Tcorrel[k,_]:.6f}' 775 for _ in range(self.N)])) 776 if D47_correl: 777 f.write(sep.join(['']+[ 778 f'{D47correl[k,_]:.0f}' 779 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 780 else f'{D47correl[k,_]:.6f}' 781 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):
784 def export(self, name, filename): 785 """ 786 Save `D47calib` object as an importable file. 787 788 ### Parameters 789 790 + **name**: 791 The name of the variable to export. 792 + **filename**: 793 The filename to write to. 794 795 ### Example 796 797 ````py 798 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 799 ```` 800 801 This should result in a `bar.py` file with the following contents: 802 803 ````py 804 foo = D47calib( 805 samples = ['LGB-2', 'DVH-2'], 806 T = [7.9, 33.7], 807 D47 = [0.6485720997671647, 0.5695972909966959], 808 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 809 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 810 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 811 label = 'Slow-growing calcites from Anderson et al. (2021)', 812 color = (0, 0.5, 0), 813 degrees = [0, 2], 814 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 815 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 816 chisq = 6.421311854486162e-27, 817 Nf = 0, 818 ) 819 ```` 820 """ 821 with open(filename, 'w') as f: 822 f.write(f''' 823{name} = D47calib( 824 samples = {self.samples}, 825 T = {self.T.tolist()}, 826 D47 = {self.D47.tolist()}, 827 sT = {self.sT.tolist()}, 828 sD47 = {self.sD47.tolist()}, 829 degrees = {self.degrees}, 830 description = {repr(self.description)}, 831 name = {repr(self.name)}, 832 label = {repr(self.label)}, 833 bfp = {({k: float(self.bfp[k]) for k in self.bfp})}, 834 bfp_CM = {self.bfp_CM.tolist()}, 835 chisq = {self.chisq}, 836 cholesky_residuals = {self.cholesky_residuals.tolist()}, 837 aic = {self.aic}, 838 bic = {self.bic}, 839 ks_pvalue = {self.ks_pvalue}, 840 ) 841''')
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=[]):
843def combine_D47calibs(calibs, degrees = [0,2], same_T = [], exclude_samples = []): 844 ''' 845 Combine data from several `D47calib` instances. 846 847 ### Parameters 848 849 + **calibs**: 850 A list of `D47calib` instances 851 + **degrees**: 852 The polynomial degrees of the combined regression. 853 + **same_T**: 854 Use this `list` to specify when samples from different calibrations are known/postulated 855 to have formed at the same temperature (e.g. `DVH-2` and `DHC2-8` from the `fiebig_2021` 856 and `anderson_2021_lsce` data sets). Each element of `same_T` is a `list` with the names 857 of two or more samples formed at the same temperature. 858 + **exclude_samples**: Use this `list` to specify the names of samples to exclude from 859 the combined calibration. 860 861 For example, the `OGLS23` calibration is computed with: 862 863 `same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]` 864 865 Note that when samples from different calibrations have the same name, 866 it is not necessary to explicitly list them in `same_T`. 867 868 Also note that the regression will fail if samples listed together in `same_T` 869 actually have different `T` values specified in the original calibrations. 870 871 ### Example 872 873 The `devils_laghetto_2023` calibration is computed using the following code: 874 875 ````py 876 K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']] 877 878 fiebig_temp = D47calib( 879 samples = [fiebig_2021.samples[_] for _ in K], 880 T = fiebig_2021.T[K], 881 D47 = fiebig_2021.D47[K], 882 sT = fiebig_2021.sT[K,:][:,K], 883 sD47 = fiebig_2021.sD47[K,:][:,K], 884 ) 885 886 devils_laghetto_2023 = combine_D47calibs( 887 calibs = [anderson_2021_lsce, fiebig_temp], 888 degrees = [0,2], 889 same_T = [{'DVH-2', 'DHC2-8'}], 890 ) 891 ```` 892 ''' 893 894 samples = [s for c in calibs for s in c.samples] 895 T = [t for c in calibs for t in c.T] 896 D47 = [x for c in calibs for x in c.D47] 897 sD47 = _block_diag(*[c.sD47 for c in calibs]) 898 sT = _block_diag(*[c.sT for c in calibs]) 899 900 for i in range(len(samples)): 901 for j in range(len(samples)): 902 if i != j: 903 if (samples[i] == samples[j] or 904 any([samples[i] in _ and samples[j] in _ for _ in same_T])): 905 906 sT[i,j] = (sT[i,i] * sT[j,j])**.5 907 908 k = [_ for _, s in enumerate(samples) if s not in exclude_samples] 909 910 calib = D47calib( 911 samples = [samples[_] for _ in k], 912 T = [T[_] for _ in k], 913 D47 = [D47[_] for _ in k], 914 sT = sT[k,:][:,k], 915 sD47 = sD47[k,:][:,k], 916 degrees = degrees, 917 ) 918 919 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>