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. This
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 (2023). 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 ETH-1 Session_01 5.994522 10.810755 16.181039 21.445703 27.780042
2 ETH-2 Session_01 -5.991066 -5.968789 -12.685206 -12.219412 -18.023381
3 ETH-3 Session_01 5.734551 11.043266 16.556578 22.005474 28.306614
4 ETH-2 Session_01 -5.981599 -6.011356 -12.728742 -12.335559 -18.023381
5 BAR-2 Session_01 -9.941731 10.985508 0.208381 21.832546 10.707292
6 ETH-3 Session_01 5.726647 11.144602 16.634507 22.275401 28.306614
7 ETH-1 Session_01 6.050787 10.800497 16.232192 21.429453 27.780042
8 FOO-1 Session_01 -0.848380 2.872996 1.535960 5.431873 4.665655
9 ETH-2 Session_01 -5.973121 -5.894178 -12.599303 -12.037594 -18.023381
10 ETH-3 Session_01 5.755765 11.179326 16.689870 22.294132 28.306614
11 BAZ-3 Session_01 -5.394955 7.938127 1.887702 15.671857 9.739724
12 ETH-1 Session_01 6.009875 10.711152 16.087994 21.275325 27.780042
13 ETH-2 Session_02 -5.993596 -6.086738 -12.816033 -12.452397 -18.023381
14 ETH-1 Session_02 6.019913 10.773652 16.165645 21.309536 27.780042
15 FOO-1 Session_02 -0.835316 2.868058 1.539147 5.483628 4.665655
16 ETH-3 Session_02 5.720443 11.187928 16.689130 22.308575 28.306614
17 ETH-1 Session_02 6.029949 10.682183 16.080186 21.182471 27.780042
18 BAZ-3 Session_02 -5.337619 7.818571 1.818572 15.400966 9.739724
19 BAR-2 Session_02 -9.951732 10.964153 0.171443 21.779172 10.707292
20 ETH-2 Session_02 -5.983050 -5.953794 -12.679809 -12.224144 -18.023381
21 ETH-3 Session_02 5.756394 11.109906 16.635806 22.189179 28.306614
22 ETH-3 Session_02 5.717666 11.175023 16.654427 22.251497 28.306614
23 ETH-2 Session_02 -5.952660 -5.887239 -12.582248 -12.081588 -18.023381
24 ETH-1 Session_02 5.995178 10.702934 16.073811 21.169324 27.780042
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) 2023 Mathieu Daëron' 23__license__ = 'MIT License - https://opensource.org/licenses/MIT' 24# __docformat__ = "restructuredtext" 25__date__ = '2023-11-15' 26__version__ = '1.3.1' 27 28 29import typer, sys 30from typing_extensions import Annotated 31import ogls as _ogls 32import numpy as _np 33from scipy.linalg import block_diag as _block_diag 34from scipy.interpolate import interp1d as _interp1d 35from matplotlib import pyplot as _ppl 36 37typer.rich_utils.STYLE_HELPTEXT = '' 38 39class D47calib(_ogls.InverseTPolynomial): 40 """ 41 Δ47 calibration class based on OGLS regression 42 of Δ47 as a polynomial function of inverse T. 43 """ 44 45 def __init__(self, 46 samples, T, D47, 47 sT = None, 48 sD47 = None, 49 degrees = [0,2], 50 xpower = 2, 51 name = '', 52 label = '', 53 description = '', 54 **kwargs, 55 ): 56 """ 57 ### Parameters 58 59 + **samples**: a list of N sample names. 60 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 61 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 62 + **sT**: uncertainties on `T`. If specified as: 63 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 64 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 65 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 66 + **sD47**: uncertainties on `D47`. If specified as: 67 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 68 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 69 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 70 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 71 + **name**: a human-readable, short name assigned to the calibration. 72 + **label**: a short description of the calibration, e.g., to be used in legends. 73 + **description**: a longer description, including relevant references/DOIs. 74 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 75 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 76 77 ### Notable attributes 78 79 + **N**: 80 The total number of observations (samples) in the calibration data. 81 + **samples**: 82 The list sample names. 83 + **T**: 84 1-D `ndarray` of temperatures in degrees C. 85 + **D47**: 86 1-D `ndarray` of Δ47 values in permil. 87 + **sT**: 88 2-D `ndarray` equal to the full (co)variance matrix for `T`. 89 + **D47**: 90 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 91 + **xpower**: 92 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 93 It is possible to change this behavior to use a different power of 1/T. 94 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 95 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 96 + **bfp**: 97 The best-fit parameters of the regression. 98 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 99 + **bff()**: 100 The best-fit polynomial function of inverse T, defined as: 101 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 102 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 103 104 105 ### Examples 106 107 A very simple example: 108 109 ````py 110 .. include:: ../../code_examples/D47calib_init/example.py 111 ```` 112 113 Should yield: 114 115 ```` 116 .. include:: ../../code_examples/D47calib_init/output.txt 117 ```` 118 119 """ 120 121 self.samples = samples[:] 122 self.name = name 123 self.label = label 124 self.description = description 125 self.D47 = _np.asarray(D47, dtype = 'float') 126 self.N = self.D47.size 127 128 if sD47 is None: 129 self.sD47 = _np.zeros((self.N, self.N)) 130 else: 131 self.sD47 = _np.asarray(sD47) 132 if len(self.sD47.shape) == 1: 133 self.sD47 = _np.diag(self.sD47**2) 134 elif len(self.sD47.shape) == 0: 135 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 136 137 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 138 139 if self.bfp is None: 140 self.regress() 141 142 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) 143 144 xi = _np.linspace(0,200**-1,1001) 145 self._inv_bff = _interp1d(self.bff(xi), xi) 146 147 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 148 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 149 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 150 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1 151 152 def __repr__(self): 153 return f'<D47calib: {self.name}>' 154 155 def invT_xaxis(self, 156 xlabel = None, 157 Ti = [0,20,50,100,250,1000], 158 ): 159 """ 160 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 161 but labeled in degrees Celsius. 162 163 ### Parameters 164 165 + **xlabel**: 166 Custom label for X axis (`r'$1\,/\,T^2$'` by default) 167 + **Ti**: 168 Specify tick locations for X axis, in degrees C. 169 170 ### Returns 171 172 + an `matplotlib.axes.Axes` instance 173 174 ### Examples 175 176 ````py 177 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 178 ```` 179 180 This should result in something like this: 181 182 <img align="center" src="example_invT_xaxis_1.png"> 183 184 It is also possible to define the X axis using a different power of 1/T 185 by first redefining the `xpower` attribute: 186 187 ````py 188 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 189 ```` 190 191 This should result in something like this: 192 193 <img align="center" src="example_invT_xaxis_2.png"> 194 """ 195 if xlabel is None: 196 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 197 _ppl.xlabel(xlabel) 198 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 199 ax = _ppl.gca() 200 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 201 ax.tick_params(which="major") 202 203 return ax 204 205 206 def plot_data(self, label = False, **kwargs): 207 """ 208 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 209 210 ### Parameters 211 212 + **label**: 213 + If `label` is a string, use this string as `label` for the underlyig 214 `matplotlib.pyplot.plot()` call. 215 + If `label = True`, use the caller's `label` attribute instead. 216 + If `label = False`, no label is specified (default behavior). 217 + **kwargs**: 218 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 219 220 ### Returns 221 222 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 223 224 ### Example 225 226 ````py 227 from matplotlib import pyplot as ppl 228 from D47calib import huyghe_2022 as calib 229 230 fig = ppl.figure(figsize = (5,3)) 231 ppl.subplots_adjust(bottom = .25, left = .15) 232 calib.invT_xaxis(Ti = [0,10,25]) 233 calib.plot_data(label = True) 234 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 235 ppl.legend() 236 ppl.savefig('example_plot_data.png', dpi = 100) 237 ````` 238 239 This should result in something like this: 240 241 <img align="center" src="example_plot_data.png"> 242 """ 243# if 'mec' not in kwargs: 244# kwargs['mec'] = self.color 245 if label is not False: 246 kwargs['label'] = self.label if label is True else label 247 return _ogls.InverseTPolynomial.plot_data(self, **kwargs) 248 249 250 def plot_error_bars(self, **kwargs): 251 """ 252 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 253 254 ### Parameters 255 256 + **kwargs**: 257 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 258 259 ### Returns 260 261 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 262 263 ### Example 264 265 ````py 266 from matplotlib import pyplot as ppl 267 from D47calib import huyghe_2022 as calib 268 269 fig = ppl.figure(figsize = (5,3)) 270 ppl.subplots_adjust(bottom = .25, left = .15) 271 calib.invT_xaxis(Ti = [0,10,25]) 272 calib.plot_error_bars(alpha = .4) 273 calib.plot_data(label = True) 274 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 275 ppl.legend() 276 ppl.savefig('example_plot_error_bars.png', dpi = 100) 277 ````` 278 279 This should result in something like this: 280 281 <img align="center" src="example_plot_error_bars.png"> 282 """ 283# if 'ecolor' not in kwargs: 284# kwargs['ecolor'] = self.color 285 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs) 286 287 288 def plot_error_ellipses(self, **kwargs): 289 """ 290 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 291 292 ### Parameters 293 294 + **kwargs**: 295 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 296 297 ### Returns 298 299 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 300 301 ### Example 302 303 ````py 304 from matplotlib import pyplot as ppl 305 from D47calib import huyghe_2022 as calib 306 307 fig = ppl.figure(figsize = (5,3)) 308 ppl.subplots_adjust(bottom = .25, left = .15) 309 calib.invT_xaxis(Ti = [0,10,25]) 310 calib.plot_error_ellipses(alpha = .4) 311 calib.plot_data(label = True) 312 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 313 ppl.legend() 314 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 315 ````` 316 317 This should result in something like this: 318 319 <img align="center" src="example_plot_error_ellipses.png"> 320 """ 321# if 'ec' not in kwargs: 322# kwargs['ec'] = self.color 323 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs) 324 325 326 def plot_bff(self, label = False, **kwargs): 327 """ 328 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 329 330 ### Parameters 331 332 + **label**: 333 + If `label` is a string, use this string as `label` for the underlyig 334 `matplotlib.pyplot.plot()` call. 335 + If `label = True`, use the caller's `label` attribute instead. 336 + If `label = False`, no label is specified (default behavior). 337 + **kwargs**: 338 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 339 340 ### Returns 341 342 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 343 344 ### Example 345 346 ````py 347 from matplotlib import pyplot as ppl 348 from D47calib import huyghe_2022 as calib 349 350 fig = ppl.figure(figsize = (5,3)) 351 ppl.subplots_adjust(bottom = .25, left = .15) 352 calib.invT_xaxis(Ti = [0,10,25]) 353 calib.plot_bff(label = True, dashes = (8,2,2,2)) 354 calib.plot_data() 355 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 356 ppl.legend() 357 ppl.savefig('example_plot_bff.png', dpi = 100) 358 ````` 359 360 This should result in something like this: 361 362 <img align="center" src="example_plot_bff.png"> 363 """ 364# if 'color' not in kwargs: 365# kwargs['color'] = self.color 366 if label is not False: 367 kwargs['label'] = self.label if label is True else label 368 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs) 369 370 371 def plot_bff_ci(self, **kwargs): 372 """ 373 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 374 375 ### Parameters 376 377 + **label**: 378 + **kwargs**: 379 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 380 381 ### Returns 382 383 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 384 385 ### Example 386 387 ````py 388 from matplotlib import pyplot as ppl 389 from D47calib import huyghe_2022 as calib 390 391 fig = ppl.figure(figsize = (5,3)) 392 ppl.subplots_adjust(bottom = .25, left = .15) 393 calib.invT_xaxis(Ti = [0,10,25]) 394 calib.plot_bff_ci(alpha = .15) 395 calib.plot_bff(label = True, dashes = (8,2,2,2)) 396 calib.plot_data() 397 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 398 ppl.legend() 399 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 400 ````` 401 402 This should result in something like this: 403 404 <img align="center" src="example_plot_bff_ci.png"> 405 """ 406# if 'color' not in kwargs: 407# kwargs['color'] = self.color 408 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs) 409 410 def T47(self, 411 D47 = None, 412 sD47 = None, 413 T=None, 414 sT = None, 415 error_from = 'both', 416 return_covar = False, 417 ): 418 ''' 419 When `D47` is input, computes corresponding T value(s). 420 `D47` input may be specified as a scalar, or as a 1-D array. 421 `T` output will then have the same type and size as `D47`. 422 423 When `T` is input, computes corresponding Δ47 value(s). 424 `T` input may be specified as a scalar, or as a 1-D array. 425 `D47` output will then have the same type and size as `T`. 426 427 Only one of either `D47` or `T` may be specified as input. 428 429 **Arguments:** 430 431 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 432 * `sD47`: Δ47 uncertainties, which may be: 433 - `None` (default) 434 - `float` or `int` (uniform standard error on `D47`) 435 - 1-D array (standard errors on `D47`) 436 - 2-D array (covariance matrix for `D47`) 437 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 438 * `sT`: T uncertainties, which may be: 439 - `None` (default) 440 - `float` or `int` (uniform standard error on `T`) 441 - 1-D array (standard errors on `T`) 442 - 2-D array (variance-covariance matrix for `T`) 443 * `error_from`: if set to `'both'` (default), returned errors take into account 444 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 445 if set to `'calib'`, only calibration uncertainties are accounted for; 446 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 447 * `return_covar`: (False by default) whether to return the full covariance matrix 448 for returned `T` or `D47` values, otherwise return standard errors for the returned 449 `T` or `D47` values instead. 450 451 **Returns (with `D47` input):** 452 453 * `T`: temperature value(s) computed from `D47` 454 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 455 456 **Returns (with `T` input):** 457 458 * `D47`: Δ47 value(s) computed from `D47` 459 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 460 461 ### Example 462 463 ````py 464 import numpy as np 465 from matplotlib import pyplot as ppl 466 from D47calib import OGLS23 as calib 467 468 X = np.linspace(1473**-2, 270**-2) 469 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 470 471 fig = ppl.figure(figsize = (5,3)) 472 ppl.subplots_adjust(bottom = .25, left = .15) 473 calib.invT_xaxis() 474 ppl.plot(X, 1000 * sD47, 'r-') 475 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 476 ppl.savefig('example_SE47.png', dpi = 100) 477 ````` 478 479 This should result in something like this: 480 481 <img src="example_SE47.png"> 482 ''' 483 484 if D47 is None and T is None: 485 raise ValueError('Either D47 or T must be specified, but both are undefined.') 486 487 if D47 is not None and T is not None: 488 raise ValueError('Either D47 or T must be specified, but not both.') 489 490 if T is not None: 491 492 D47 = self._D47_from_T(T) 493 Np = len(self.degrees) 494 N = D47.size 495 496 ### Compute covariance matrix of (*bfp, *T): 497 CM = _np.zeros((Np+N, Np+N)) 498 499 if error_from in ['calib', 'both']: 500 CM[:Np, :Np] = self.bfp_CM[:,:] 501 502 if (sT is not None) and error_from in ['sT', 'both']: 503 _sT = _np.asarray(sT) 504 if _sT.ndim == 0: 505 for k in range(N): 506 CM[Np+k, Np+k] = _sT**2 507 elif _sT.ndim == 1: 508 for k in range(N): 509 CM[Np+k, Np+k] = _sT[k]**2 510 elif _sT.ndim == 2: 511 CM[-N:, -N:] = _sT[:,:] 512 513 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 514 _T = _np.asarray(T) 515 if _T.ndim == 0: 516 _T = _np.expand_dims(_T, 0) 517 J = _np.zeros((N, Np+N)) 518 519 if (sT is not None) and error_from in ['sT', 'both']: 520 for k in range(N): 521 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 522 523 if error_from in ['calib', 'both']: 524 525 for k in range(Np): 526 527 p1 = {_: self.bfp[_] for _ in self.bfp} 528 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 529 530 p2 = {_: self.bfp[_] for _ in self.bfp} 531 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 532 533 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) 534 535 ### Error propagation: 536 CM_D47 = J @ CM @ J.T 537 538 if return_covar: 539 return D47, CM_D47 540 else: 541 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 542 543 if D47 is not None: 544 545 T = self._T_from_D47(D47) 546 Np = len(self.degrees) 547 N = T.size 548 549 ### Compute covariance matrix of (*bfp, *T): 550 CM = _np.zeros((Np+N, Np+N)) 551 552 if error_from in ['calib', 'both']: 553 CM[:Np, :Np] = self.bfp_CM[:,:] 554 555 if (sD47 is not None) and error_from in ['sD47', 'both']: 556 _sD47 = _np.asarray(sD47) 557 if _sD47.ndim == 0: 558 for k in range(N): 559 CM[Np+k, Np+k] = _sD47**2 560 elif _sD47.ndim == 1: 561 for k in range(N): 562 CM[Np+k, Np+k] = _sD47[k]**2 563 elif _sD47.ndim == 2: 564 CM[-N:, -N:] = _sD47[:,:] 565 566 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 567 _D47 = _np.asarray(D47) 568 if _D47.ndim == 0: 569 _D47 = _np.expand_dims(_D47, 0) 570 J = _np.zeros((N, Np+N)) 571 if (sD47 is not None) and error_from in ['sD47', 'both']: 572 for k in range(N): 573 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 574 if error_from in ['calib', 'both']: 575 576 xi = _np.linspace(0,200**-1,1001)[1:] 577 for k in range(Np): 578 579 if self.bfp_CM[k,k]: 580 _epsilon_ = self.bfp_CM[k,k]**.5 581 else: 582 _epsilon_ = 1e-6 583 584 p1 = {_: self.bfp[_] for _ in self.bfp} 585 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 586 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 587 588 p2 = {_: self.bfp[_] for _ in self.bfp} 589 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 590 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 591 592 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 593 594 ### Error propagation: 595 CM_T = J @ CM @ J.T 596 597 if return_covar: 598 return T, CM_T 599 else: 600 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5 601 602 603 def plot_T47_errors( 604 self, 605 calibname = None, 606 rD47 = 0.010, 607 Nr = [2,4,8,12,20], 608 Tmin = 0, 609 Tmax = 120, 610 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 611 yscale = 'lin', 612 ): 613 """ 614 Plot SE of T reconstructed using the calibration as a function of T for various 615 combinations of analytical precision and number of analytical replicates. 616 617 **Arguments** 618 619 + **calibname**: 620 Which calibration name to display. By default, use `label` attribute. 621 + **rD47**: 622 Analytical precision of a single analysis. 623 + **Nr**: 624 A list of lines to plot, each corresponding to a given number of replicates. 625 + **Tmin**: 626 Minimum T to plot. 627 + **Tmax**: 628 Maximum T to plot. 629 + **colors**: 630 A list of colors to distinguish the plotted lines. 631 + **yscale**: 632 + If `'lin'`, the Y axis uses a linear scale. 633 + If `'log'`, the Y axis uses a logarithmic scale. 634 635 **Example** 636 637 ````py 638 from matplotlib import pyplot as ppl 639 from D47calib import devils_laghetto_2023 as calib 640 641 fig = ppl.figure(figsize = (3.5,4)) 642 ppl.subplots_adjust(bottom = .2, left = .15) 643 calib.plot_T47_errors( 644 calibname = 'Devils Laghetto calibration', 645 Nr = [1,2,4,16], 646 Tmin =0, 647 Tmax = 40, 648 ) 649 ppl.savefig('example_SE_T.png', dpi = 100) 650 ```` 651 652 This should result in something like this: 653 654 <img src="example_SE_T.png"> 655 """ 656 657 if calibname is None: 658 calibname = self.label 659 660 Nr = _np.array(Nr) 661 if len(colors) < Nr.size: 662 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 663 from colorsys import hsv_to_rgb 664 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 665 colors = [hsv_to_rgb(*x) for x in hsv] 666 667 Ti = _np.linspace(Tmin, Tmax) 668 D47i, _ = self.T47(T = Ti) 669 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 670 671 ymax, ymin = 0, 1e6 672 for N,c in zip(Nr, colors): 673 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 674 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 675 ymin = min(ymin, min(sT)) 676 ymax = max(ymax, max(sT)) 677 678 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 679 680 _ppl.legend(fontsize=9) 681 _ppl.xlabel("T (°C)") 682 683 _ppl.ylabel("Standard error on reconstructed T (°C)") 684 685 # yticks([0,.5,1,1.5,2]) 686 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 687 _ppl.grid( alpha = .25) 688 if yscale == 'lin': 689 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 690 t1, t2 = self.T.min(), self.T.max() 691 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 692 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 693 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 694 elif yscale == 'log': 695 ymin /= 2 696 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 697 _ppl.yscale('log') 698 t1, t2 = self.T.min(), self.T.max() 699 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 700 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 701 702 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 703 """ 704 Write calibration data to a csv file. 705 706 ### Parameters 707 708 + **csvfile**: 709 The filename to write data to. 710 + **sep**: 711 The separator between CSV fields. 712 + **label**: 713 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 714 + If specified as a `str`, include a `Dataset` column with that string. 715 + If specified as `False`, do not include a `Dataset` column. 716 + **T_correl**: 717 + If `True`, include correlations between all `T` values. 718 + **D47_correl**: 719 + If `True`, include correlations between all `D47` values. 720 721 ### Example 722 723 ````py 724 D47calib.huyghe_2022.export_data( 725 csvfile = 'example_export_data.csv', 726 T_correl = True, 727 D47_correl = True, 728 ) 729 ```` 730 731 This should result in something like this ([link](example_export_data.csv)): 732 733 .. include:: ../../docs/example_export_data.md 734 735 """ 736 n = len(str(self.N)) 737 738 with open(csvfile, 'w') as f: 739 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 740 741 if label: 742 f.write(f'{sep}Dataset') 743 744 if T_correl: 745 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 746 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 747 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 748 749 if D47_correl: 750 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 751 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 752 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 753 754 for k, (s, T, sT, D47, sD47) in enumerate(zip( 755 self.samples, 756 self.T, 757 _np.diag(self.sT)**.5, 758 self.D47, 759 _np.diag(self.sD47)**.5, 760 )): 761 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 762 if label: 763 if label is True: 764 f.write(f'{sep}{self.label}') 765 else: 766 f.write(f'{sep}{label}') 767 if T_correl: 768 f.write(sep.join(['']+[ 769 f'{Tcorrel[k,_]:.0f}' 770 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 771 else f'{Tcorrel[k,_]:.6f}' 772 for _ in range(self.N)])) 773 if D47_correl: 774 f.write(sep.join(['']+[ 775 f'{D47correl[k,_]:.0f}' 776 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 777 else f'{D47correl[k,_]:.6f}' 778 for _ in range(self.N)])) 779 780 781 def export(self, name, filename): 782 """ 783 Save `D47calib` object as an importable file. 784 785 ### Parameters 786 787 + **name**: 788 The name of the variable to export. 789 + **filename**: 790 The filename to write to. 791 792 ### Example 793 794 ````py 795 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 796 ```` 797 798 This should result in a `bar.py` file with the following contents: 799 800 ````py 801 foo = D47calib( 802 samples = ['LGB-2', 'DVH-2'], 803 T = [7.9, 33.7], 804 D47 = [0.6485720997671647, 0.5695972909966959], 805 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 806 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 807 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 808 label = 'Slow-growing calcites from Anderson et al. (2021)', 809 color = (0, 0.5, 0), 810 degrees = [0, 2], 811 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 812 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 813 chisq = 6.421311854486162e-27, 814 Nf = 0, 815 ) 816 ```` 817 """ 818 with open(filename, 'w') as f: 819 f.write(f''' 820{name} = D47calib( 821 samples = {self.samples}, 822 T = {list(self.T)}, 823 D47 = {list(self.D47)}, 824 sT = {[list(l) for l in self.sT]}, 825 sD47 = {[list(l) for l in self.sD47]}, 826 degrees = {self.degrees}, 827 description = {repr(self.description)}, 828 name = {repr(self.name)}, 829 label = {repr(self.label)}, 830 bfp = {self.bfp}, 831 bfp_CM = {[list(l) for l in self.bfp_CM]}, 832 chisq = {self.chisq}, 833 cholesky_residuals = {list(self.cholesky_residuals)}, 834 aic = {self.aic}, 835 bic = {self.bic}, 836 ks_pvalue = {self.ks_pvalue}, 837 ) 838''') 839 840def combine_D47calibs(calibs, degrees = [0,2], same_T = [], exclude_samples = []): 841 ''' 842 Combine data from several `D47calib` instances. 843 844 ### Parameters 845 846 + **calibs**: 847 A list of `D47calib` instances 848 + **degrees**: 849 The polynomial degrees of the combined regression. 850 + **same_T**: 851 Use this `list` to specify when samples from different calibrations are known/postulated 852 to have formed at the same temperature (e.g. `DVH-2` and `DHC2-8` from the `fiebig_2021` 853 and `anderson_2021_lsce` data sets). Each element of `same_T` is a `list` with the names 854 of two or more samples formed at the same temperature. 855 + **exclude_samples**: Use this `list` to specify the names of samples to exclude from 856 the combined calibration. 857 858 For example, the `OGLS23` calibration is computed with: 859 860 `same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]` 861 862 Note that when samples from different calibrations have the same name, 863 it is not necessary to explicitly list them in `same_T`. 864 865 Also note that the regression will fail if samples listed together in `same_T` 866 actually have different `T` values specified in the original calibrations. 867 868 ### Example 869 870 The `devils_laghetto_2023` calibration is computed using the following code: 871 872 ````py 873 K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']] 874 875 fiebig_temp = D47calib( 876 samples = [fiebig_2021.samples[_] for _ in K], 877 T = fiebig_2021.T[K], 878 D47 = fiebig_2021.D47[K], 879 sT = fiebig_2021.sT[K,:][:,K], 880 sD47 = fiebig_2021.sD47[K,:][:,K], 881 ) 882 883 devils_laghetto_2023 = combine_D47calibs( 884 calibs = [anderson_2021_lsce, fiebig_temp], 885 degrees = [0,2], 886 same_T = [{'DVH-2', 'DHC2-8'}], 887 ) 888 ```` 889 ''' 890 891 samples = [s for c in calibs for s in c.samples] 892 T = [t for c in calibs for t in c.T] 893 D47 = [x for c in calibs for x in c.D47] 894 sD47 = _block_diag(*[c.sD47 for c in calibs]) 895 sT = _block_diag(*[c.sT for c in calibs]) 896 897 for i in range(len(samples)): 898 for j in range(len(samples)): 899 if i != j: 900 if (samples[i] == samples[j] or 901 any([samples[i] in _ and samples[j] in _ for _ in same_T])): 902 903 sT[i,j] = (sT[i,i] * sT[j,j])**.5 904 905 k = [_ for _, s in enumerate(samples) if s not in exclude_samples] 906 907 calib = D47calib( 908 samples = [samples[_] for _ in k], 909 T = [T[_] for _ in k], 910 D47 = [D47[_] for _ in k], 911 sT = sT[k,:][:,k], 912 sD47 = sD47[k,:][:,k], 913 degrees = degrees, 914 ) 915 916 return calib 917 918from ._calibs import * 919 920def _covar2correl(C): 921 SE = _np.diag(C)**.5 922 return SE, _np.diag(SE**-1) @ C @ _np.diag(SE**-1) 923 924try: 925 app = typer.Typer( 926 add_completion = False, 927 context_settings={'help_option_names': ['-h', '--help']}, 928 rich_markup_mode = 'rich', 929 ) 930 931 @app.command() 932 def _cli( 933 input: Annotated[str, typer.Argument(help = "Specify either the path of an input file or just '-' to read input from stdin")] = '-', 934 include_samples: Annotated[str, typer.Option('--include-samples', '-u', help = 'Only include samples listed in this file')] = 'all', 935 exclude_samples: Annotated[str, typer.Option('--exclude-samples', '-x', help = 'Exclude samples listed in this file')] = 'none', 936 outfile: Annotated[str, typer.Option('--output-file', '-o', help = 'Write output to this file instead of printing to stdout')] = 'none', 937 calib: Annotated[str, typer.Option('--calib', '-c', help = 'D47 calibration function to use')] = 'OGLS23', 938 delim_in: Annotated[str, typer.Option('--delimiter-in', '-i', help = "Delimiter used in the input.")] = ',', 939 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", 940 T_precision: Annotated[int, typer.Option('--T-precision', '-p', help = 'Precision for T output')] = 2, 941 D47_precision: Annotated[int, typer.Option('--D47-precision', '-q', help = 'Precision for D47 output')] = 4, 942 correl_precision: Annotated[int, typer.Option('--correl-precision', '-r', help = 'Precision for correlation output')] = 3, 943 covar_precision: Annotated[int, typer.Option('--covar-precision', '-s', help = 'Precision for covariance output')] = 3, 944 return_covar: Annotated[bool, typer.Option('--return-covar', '-v', help = 'Output covariance matrix instead of correlation matrix')] = False, 945 ignore_correl: Annotated[bool, typer.Option('--ignore-correl', '-g', help = 'Only consider and report standard errors without correlations')] = False, 946 version: Annotated[bool, typer.Option('--version', '-V', help = 'Print D47calib version')] = False, 947 ): 948 """ 949[b]Purpose:[/b] 950 951Reads data from an input file, converts between T and D47 values, and print out the results. 952 953The 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>. 954 955The first line of the input file must be one of the following: 956 957- [b]Option 1:[/b] T 958- [b]Option 2:[/b] T<sep>T_SE 959- [b]Option 3:[/b] T<sep>T_SE<sep>T_correl 960- [b]Option 4:[/b] T<sep>T_covar 961- [b]Option 5:[/b] D47 962- [b]Option 6:[/b] D47<sep>D47_SE 963- [b]Option 7:[/b] D47<sep>D47_SE<sep>D47_correl 964- [b]Option 8:[/b] D47<sep>D47_covar 965 966The rest of the input must be any number of lines with float values corresponding to the fields in the first line, separated by <sep>. 967 968[bold]Example input file:[/bold] 969 970[italic]D47 D47_SE D47_correl[/italic] 971[italic]0.6324 0.0101 1.00 0.25 0.25[/italic] 972[italic]0.6281 0.0087 0.25 1.00 0.25[/italic] 973[italic]0.6385 0.0095 0.25 0.25 1.00[/italic] 974 975The corresponding D47 (options 1-4) or T (options 4-8) values are computed, along with standard errors and error correlations from calibration uncertainties. 976 977For 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. 978 979The example above will thus result in an output with the following fields: 980 981[italic]- D47[/italic] 982[italic]- D47_SE[/italic] 983[italic]- D47_correl[/italic] 984[italic]- T[/italic] 985[italic]- T_SE_from_calib[/italic] 986[italic]- T_correl_from_calib[/italic] 987[italic]- T_SE_from_input[/italic] 988[italic]- T_correl_from_input[/italic] 989[italic]- T_SE_from_both[/italic] 990[italic]- T_correl_from_both[/italic] 991 992Results may also be saved to a file using [bold]--output-file <filename>[/bold] or [bold]-o <filename>[/bold]. 993 994To 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. 995Then 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. 996To 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. 997""" 998 999 if version: 1000 print(__version__) 1001 return None 1002 1003 ### INCOMPATIBILITY BETWEEN --ignore-correl AND --return-covar 1004 if ignore_correl: 1005 return_covar = False 1006 1007 ### USE WHITESPACE AS INPUT DELIMITER 1008 if delim_in == ' ': 1009 delim_in = None 1010 1011 ### SMART SELECTION OF OUTPUT DELIMITER 1012 if delim_out == "',' when writing to output file, '>' when printing to stdout": 1013 if outfile == 'none': 1014 delim_out = '>' 1015 else: 1016 delim_out = ',' 1017 1018 ### CALIBRATION 1019 if calib in globals() and type(globals()[calib]) == D47calib: 1020 calib = globals()[calib] 1021 else: 1022 with open(calib) as f: 1023 calibdata = _np.array([[c.strip() for c in l.strip().split(delim_in)] for l in f.readlines()[1:]], dtype = float) 1024 1025 degrees = [int(d) for d in calibdata[:,0]] 1026 bfp = {f'a{k}': a for k,a in zip(degrees, calibdata[:,1])} 1027 bfp_CM = calibdata[:,2:] 1028 if bfp_CM.shape[0] != bfp_CM.shape[1]: 1029 bfp_CM = _np.zeros((len(degrees), len(degrees))) 1030 calib = D47calib( 1031 samples = [], T = [], sT = [], D47 = [], sD47 = [], 1032 degrees = degrees, bfp = bfp, bfp_CM = bfp_CM, 1033 ) 1034 1035 ### READ INPUT STRINGS 1036 if input == '-': 1037 data = [[c.strip() for c in l.strip().split(delim_in)] for l in sys.stdin] 1038 else: 1039 with open(input) as f: 1040 data = [[c.strip() for c in l.strip().split(delim_in)] for l in f.readlines()] 1041 1042 if include_samples == 'all': 1043 samples_to_include = [] 1044 else: 1045 with open(include_samples) as f: 1046 samples_to_include = [l.strip() for l in f.readlines()] 1047 1048 if exclude_samples == 'none': 1049 samples_to_exclude = [] 1050 else: 1051 with open(exclude_samples) as f: 1052 samples_to_exclude = [l.strip() for l in f.readlines()] 1053 1054 if len(samples_to_include) > 0 or len(samples_to_exclude) > 0: 1055 try: 1056 k = data[0].index('Sample') 1057 except ValueError: 1058 raise KeyError("When using options --include-samples or --exclude-samples, the input file must have a column labeled 'Sample'.") 1059 1060 if len(samples_to_include) > 0: 1061 data = [data[0]] + [l for l in data[1:] if l[k] in samples_to_include] 1062 data = [data[0]] + [l for l in data[1:] if l[k] not in samples_to_exclude] 1063 1064 ### FIND FIRST DATA COLUMN 1065 k = 0 1066 while data[0][k] not in ['T', 'D47']: 1067 k += 1 1068 if k == len(data[0]): 1069 raise KeyError("None of the input column headers are 'T' or 'D47'.") 1070 data_out = [l[:k] for l in data] 1071 data = [l[k:] for l in data] 1072 1073 ### READ INPUT FIELDS 1074 fields = data[0] 1075 1076 ### CHECK FOR UNSUPPORTED FIELD COMBINATIONS 1077 if fields not in [ 1078 ['T'], 1079 ['T', 'T_SE'], 1080 ['T', 'T_covar'], 1081 ['T', 'T_SE', 'T_correl'], 1082 ['D47'], 1083 ['D47', 'D47_SE'], 1084 ['D47', 'D47_covar'], 1085 ['D47', 'D47_SE', 'D47_correl'], 1086 ]: 1087 raise KeyError("There is a problem with the combination of field names provided as input.") 1088 1089 ### BOOK-KEEPING 1090 infield = fields[0] 1091 X_precision = {'T': T_precision, 'D47': D47_precision}[infield] 1092 outfield = {'T': 'D47', 'D47': 'T'}[infield] 1093 Y_precision = {'T': T_precision, 'D47': D47_precision}[outfield] 1094 N = len(data)-1 1095 1096 ### READ INPUT DATA, ALSO SAVING ITS ORIGINAL STRINGS 1097 X_str = [l[0] for l in data[1:]] 1098 X = _np.array(X_str, dtype = float) 1099 1100 if len(fields) == 1: 1101 X_SE = X*0 1102 X_correl = _np.eye(N) 1103 X_covar = _np.zeros((N, N)) 1104 X_SE_str = [f'{c:.{X_precision}f}' for c in X_SE] 1105 X_correl_str = [[f'{c:.{correl_precision}f}' for c in l] for l in X_correl] 1106 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1107 if len(fields) == 2: 1108 if fields[1].endswith('_SE'): 1109 X_SE_str = [l[1] for l in data[1:]] 1110 X_SE = _np.array(X_SE_str, dtype = float) 1111 X_covar = _np.diag(X_SE**2) 1112 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1113 elif fields[1].endswith('_covar'): 1114 X_covar_str = [l[1:N+1] for l in data[1:]] 1115 X_covar = _np.array(X_covar_str, dtype = float) 1116 X_SE = _np.diag(X_covar)**.5 1117 X_SE_str = [f'{c:.{X_precision}f}' for c in X_SE] 1118 X_correl = _np.diag(X_SE**-1) @ X_covar @ _np.diag(X_SE**-1) 1119 X_correl_str = [[f'{c:.{correl_precision}f}' for c in l] for l in X_correl] 1120 elif len(fields) == 3: 1121 X_SE_str = [l[1] for l in data[1:]] 1122 X_SE = _np.array(X_SE_str, dtype = float) 1123 X_correl_str = [l[2:N+2] for l in data[1:]] 1124 X_correl = _np.array(X_correl_str, dtype = float) 1125 X_covar = _np.diag(X_SE) @ X_correl @ _np.diag(X_SE) 1126 X_covar_str = [[f'{c:.{covar_precision}e}' for c in l] for l in X_covar] 1127 1128 ### COMPUTE OUTPUT VALUES AND COVAR 1129 kwargs = {infield: X, f's{infield}': X_covar} 1130 Y, Y_covar_from_calib = calib.T47(**kwargs, error_from = 'calib', return_covar = True) 1131 Y, Y_covar_from_input = calib.T47(**kwargs, error_from = f's{infield}', return_covar = True) 1132 Y, Y_covar_from_both = calib.T47(**kwargs, error_from = 'both', return_covar = True) 1133 1134 Y_SE_from_calib = _np.diag(Y_covar_from_calib)**.5 1135 Y_SE_from_input = _np.diag(Y_covar_from_input)**.5 1136 Y_SE_from_both = _np.diag(Y_covar_from_both)**.5 1137 1138 if (Y_SE_from_calib**2).min(): 1139 Y_correl_from_calib = _np.diag(Y_SE_from_calib**-1) @ Y_covar_from_calib @ _np.diag(Y_SE_from_calib**-1) 1140 else: 1141 Y_correl_from_calib = _np.eye(N) 1142 1143 if (Y_SE_from_input**2).min(): 1144 Y_correl_from_input = _np.diag(Y_SE_from_input**-1) @ Y_covar_from_input @ _np.diag(Y_SE_from_input**-1) 1145 else: 1146 Y_correl_from_input = _np.eye(N) 1147 1148 if (Y_SE_from_both**2).min(): 1149 Y_correl_from_both = _np.diag(Y_SE_from_both**-1) @ Y_covar_from_both @ _np.diag(Y_SE_from_both**-1) 1150 else: 1151 Y_correl_from_both = _np.eye(N) 1152 1153 ### BUILD Y STRINGS 1154 Y_str = [f'{y:.{Y_precision}f}' for y in Y] 1155 1156 Y_SE_from_calib_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_calib] 1157 Y_SE_from_input_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_input] 1158 Y_SE_from_both_str = [f'{sy:.{Y_precision}f}' for sy in Y_SE_from_both] 1159 1160 Y_covar_from_calib_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_calib] 1161 Y_covar_from_input_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_input] 1162 Y_covar_from_both_str = [[f'{c:.{covar_precision}e}' for c in l] for l in Y_covar_from_both] 1163 1164 Y_correl_from_calib_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_calib] 1165 Y_correl_from_input_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_input] 1166 Y_correl_from_both_str = [[f'{c:.{correl_precision}f}' for c in l] for l in Y_correl_from_both] 1167 1168 ### ADD SE COLUMN TO INPUT 1169 if f'{infield}_covar' in fields: 1170 fields.insert(1, f'{infield}_SE') 1171 1172 ### ADD X COLUMNS TO OUTPUT DATA 1173 data_out[0] += [infield] 1174 for k in range(N): 1175 data_out[k+1] += [X_str[k]] 1176 for f in fields[1:]: 1177 if f.endswith('_SE'): 1178 data_out[0] += [f] 1179 for k in range(N): 1180 data_out[k+1] += [X_SE_str[k]] 1181 if f.endswith('_covar') or f.endswith('_correl'): 1182 if not ignore_correl: 1183 data_out[0] += [f] + ['' for _ in range(N-1)] 1184 for k in range(N): 1185 data_out[k+1] += (X_covar_str if f.endswith('_covar') else X_correl_str)[k][:] 1186 1187 ### ADD Y COLUMNS TO OUTPUT DATA 1188 data_out[0] += [outfield] 1189 for k in range(N): 1190 data_out[k+1] += [Y_str[k]] 1191 1192 data_out[0] += [f'{outfield}_SE_from_calib'] 1193 for k in range(N): 1194 data_out[k+1] += [Y_SE_from_calib_str[k]] 1195 if not ignore_correl: 1196 if return_covar: 1197 data_out[0] += [f'{outfield}_covar_from_calib'] + ['' for _ in range(N-1)] 1198 for k in range(N): 1199 data_out[k+1] += Y_covar_from_calib_str[k] 1200 else: 1201 data_out[0] += [f'{outfield}_correl_from_calib'] + ['' for _ in range(N-1)] 1202 for k in range(N): 1203 data_out[k+1] += Y_correl_from_calib_str[k] 1204 1205 data_out[0] += [f'{outfield}_SE_from_input'] 1206 for k in range(N): 1207 data_out[k+1] += [Y_SE_from_input_str[k]] 1208 if not ignore_correl: 1209 if return_covar: 1210 data_out[0] += [f'{outfield}_covar_from_input'] + ['' for _ in range(N-1)] 1211 for k in range(N): 1212 data_out[k+1] += Y_covar_from_input_str[k] 1213 else: 1214 data_out[0] += [f'{outfield}_correl_from_input'] + ['' for _ in range(N-1)] 1215 for k in range(N): 1216 data_out[k+1] += Y_correl_from_input_str[k] 1217 1218 data_out[0] += [f'{outfield}_SE_from_both'] 1219 for k in range(N): 1220 data_out[k+1] += [Y_SE_from_both_str[k]] 1221 if not ignore_correl: 1222 if return_covar: 1223 data_out[0] += [f'{outfield}_covar_from_both'] + ['' for _ in range(N-1)] 1224 for k in range(N): 1225 data_out[k+1] += Y_covar_from_both_str[k] 1226 else: 1227 data_out[0] += [f'{outfield}_correl_from_both'] + ['' for _ in range(N-1)] 1228 for k in range(N): 1229 data_out[k+1] += Y_correl_from_both_str[k] 1230 1231 1232 ### PRINT OUTPUT TO STDOUT OR SAVE IT TO FILE 1233 if delim_out in '<>': 1234 lengths = [max([len(data_out[j][k]) for j in range(len(data_out))]) for k in range(len(data_out[0]))] 1235 1236 txt = '' 1237 for l in data_out: 1238 for k in range(len(l)): 1239 if k > 0: 1240 txt += ' ' 1241 txt += f'{l[k]:{delim_out}{lengths[k]}s}' 1242 txt += '\n' 1243 1244 txt = txt[:-1] 1245 1246 else: 1247 txt = '\n'.join([delim_out.join(l) for l in data_out]) 1248 1249 if outfile == 'none': 1250 print(txt) 1251 else: 1252 with open(outfile, 'w') as f: 1253 f.write(txt) 1254 1255 def __cli(): 1256 app() 1257 1258except NameError: 1259 pass 1260
Reference
class
D47calib(ogls.InverseTPolynomial):
40class D47calib(_ogls.InverseTPolynomial): 41 """ 42 Δ47 calibration class based on OGLS regression 43 of Δ47 as a polynomial function of inverse T. 44 """ 45 46 def __init__(self, 47 samples, T, D47, 48 sT = None, 49 sD47 = None, 50 degrees = [0,2], 51 xpower = 2, 52 name = '', 53 label = '', 54 description = '', 55 **kwargs, 56 ): 57 """ 58 ### Parameters 59 60 + **samples**: a list of N sample names. 61 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 62 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 63 + **sT**: uncertainties on `T`. If specified as: 64 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 65 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 66 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 67 + **sD47**: uncertainties on `D47`. If specified as: 68 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 69 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 70 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 71 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 72 + **name**: a human-readable, short name assigned to the calibration. 73 + **label**: a short description of the calibration, e.g., to be used in legends. 74 + **description**: a longer description, including relevant references/DOIs. 75 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 76 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 77 78 ### Notable attributes 79 80 + **N**: 81 The total number of observations (samples) in the calibration data. 82 + **samples**: 83 The list sample names. 84 + **T**: 85 1-D `ndarray` of temperatures in degrees C. 86 + **D47**: 87 1-D `ndarray` of Δ47 values in permil. 88 + **sT**: 89 2-D `ndarray` equal to the full (co)variance matrix for `T`. 90 + **D47**: 91 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 92 + **xpower**: 93 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 94 It is possible to change this behavior to use a different power of 1/T. 95 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 96 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 97 + **bfp**: 98 The best-fit parameters of the regression. 99 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 100 + **bff()**: 101 The best-fit polynomial function of inverse T, defined as: 102 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 103 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 104 105 106 ### Examples 107 108 A very simple example: 109 110 ````py 111 .. include:: ../../code_examples/D47calib_init/example.py 112 ```` 113 114 Should yield: 115 116 ```` 117 .. include:: ../../code_examples/D47calib_init/output.txt 118 ```` 119 120 """ 121 122 self.samples = samples[:] 123 self.name = name 124 self.label = label 125 self.description = description 126 self.D47 = _np.asarray(D47, dtype = 'float') 127 self.N = self.D47.size 128 129 if sD47 is None: 130 self.sD47 = _np.zeros((self.N, self.N)) 131 else: 132 self.sD47 = _np.asarray(sD47) 133 if len(self.sD47.shape) == 1: 134 self.sD47 = _np.diag(self.sD47**2) 135 elif len(self.sD47.shape) == 0: 136 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 137 138 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 139 140 if self.bfp is None: 141 self.regress() 142 143 self._bff_deriv = lambda x: _np.array([k * self.bfp[f'a{k}'] * x**(k-1) for k in degrees if k > 0]).sum(axis = 0) 144 145 xi = _np.linspace(0,200**-1,1001) 146 self._inv_bff = _interp1d(self.bff(xi), xi) 147 148 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 149 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 150 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 151 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1 152 153 def __repr__(self): 154 return f'<D47calib: {self.name}>' 155 156 def invT_xaxis(self, 157 xlabel = None, 158 Ti = [0,20,50,100,250,1000], 159 ): 160 """ 161 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 162 but labeled in degrees Celsius. 163 164 ### Parameters 165 166 + **xlabel**: 167 Custom label for X axis (`r'$1\,/\,T^2$'` by default) 168 + **Ti**: 169 Specify tick locations for X axis, in degrees C. 170 171 ### Returns 172 173 + an `matplotlib.axes.Axes` instance 174 175 ### Examples 176 177 ````py 178 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 179 ```` 180 181 This should result in something like this: 182 183 <img align="center" src="example_invT_xaxis_1.png"> 184 185 It is also possible to define the X axis using a different power of 1/T 186 by first redefining the `xpower` attribute: 187 188 ````py 189 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 190 ```` 191 192 This should result in something like this: 193 194 <img align="center" src="example_invT_xaxis_2.png"> 195 """ 196 if xlabel is None: 197 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 198 _ppl.xlabel(xlabel) 199 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 200 ax = _ppl.gca() 201 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 202 ax.tick_params(which="major") 203 204 return ax 205 206 207 def plot_data(self, label = False, **kwargs): 208 """ 209 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 210 211 ### Parameters 212 213 + **label**: 214 + If `label` is a string, use this string as `label` for the underlyig 215 `matplotlib.pyplot.plot()` call. 216 + If `label = True`, use the caller's `label` attribute instead. 217 + If `label = False`, no label is specified (default behavior). 218 + **kwargs**: 219 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 220 221 ### Returns 222 223 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 224 225 ### Example 226 227 ````py 228 from matplotlib import pyplot as ppl 229 from D47calib import huyghe_2022 as calib 230 231 fig = ppl.figure(figsize = (5,3)) 232 ppl.subplots_adjust(bottom = .25, left = .15) 233 calib.invT_xaxis(Ti = [0,10,25]) 234 calib.plot_data(label = True) 235 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 236 ppl.legend() 237 ppl.savefig('example_plot_data.png', dpi = 100) 238 ````` 239 240 This should result in something like this: 241 242 <img align="center" src="example_plot_data.png"> 243 """ 244# if 'mec' not in kwargs: 245# kwargs['mec'] = self.color 246 if label is not False: 247 kwargs['label'] = self.label if label is True else label 248 return _ogls.InverseTPolynomial.plot_data(self, **kwargs) 249 250 251 def plot_error_bars(self, **kwargs): 252 """ 253 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 254 255 ### Parameters 256 257 + **kwargs**: 258 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 259 260 ### Returns 261 262 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 263 264 ### Example 265 266 ````py 267 from matplotlib import pyplot as ppl 268 from D47calib import huyghe_2022 as calib 269 270 fig = ppl.figure(figsize = (5,3)) 271 ppl.subplots_adjust(bottom = .25, left = .15) 272 calib.invT_xaxis(Ti = [0,10,25]) 273 calib.plot_error_bars(alpha = .4) 274 calib.plot_data(label = True) 275 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 276 ppl.legend() 277 ppl.savefig('example_plot_error_bars.png', dpi = 100) 278 ````` 279 280 This should result in something like this: 281 282 <img align="center" src="example_plot_error_bars.png"> 283 """ 284# if 'ecolor' not in kwargs: 285# kwargs['ecolor'] = self.color 286 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs) 287 288 289 def plot_error_ellipses(self, **kwargs): 290 """ 291 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 292 293 ### Parameters 294 295 + **kwargs**: 296 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 297 298 ### Returns 299 300 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 301 302 ### Example 303 304 ````py 305 from matplotlib import pyplot as ppl 306 from D47calib import huyghe_2022 as calib 307 308 fig = ppl.figure(figsize = (5,3)) 309 ppl.subplots_adjust(bottom = .25, left = .15) 310 calib.invT_xaxis(Ti = [0,10,25]) 311 calib.plot_error_ellipses(alpha = .4) 312 calib.plot_data(label = True) 313 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 314 ppl.legend() 315 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 316 ````` 317 318 This should result in something like this: 319 320 <img align="center" src="example_plot_error_ellipses.png"> 321 """ 322# if 'ec' not in kwargs: 323# kwargs['ec'] = self.color 324 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs) 325 326 327 def plot_bff(self, label = False, **kwargs): 328 """ 329 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 330 331 ### Parameters 332 333 + **label**: 334 + If `label` is a string, use this string as `label` for the underlyig 335 `matplotlib.pyplot.plot()` call. 336 + If `label = True`, use the caller's `label` attribute instead. 337 + If `label = False`, no label is specified (default behavior). 338 + **kwargs**: 339 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 340 341 ### Returns 342 343 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 344 345 ### Example 346 347 ````py 348 from matplotlib import pyplot as ppl 349 from D47calib import huyghe_2022 as calib 350 351 fig = ppl.figure(figsize = (5,3)) 352 ppl.subplots_adjust(bottom = .25, left = .15) 353 calib.invT_xaxis(Ti = [0,10,25]) 354 calib.plot_bff(label = True, dashes = (8,2,2,2)) 355 calib.plot_data() 356 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 357 ppl.legend() 358 ppl.savefig('example_plot_bff.png', dpi = 100) 359 ````` 360 361 This should result in something like this: 362 363 <img align="center" src="example_plot_bff.png"> 364 """ 365# if 'color' not in kwargs: 366# kwargs['color'] = self.color 367 if label is not False: 368 kwargs['label'] = self.label if label is True else label 369 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs) 370 371 372 def plot_bff_ci(self, **kwargs): 373 """ 374 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 375 376 ### Parameters 377 378 + **label**: 379 + **kwargs**: 380 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 381 382 ### Returns 383 384 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 385 386 ### Example 387 388 ````py 389 from matplotlib import pyplot as ppl 390 from D47calib import huyghe_2022 as calib 391 392 fig = ppl.figure(figsize = (5,3)) 393 ppl.subplots_adjust(bottom = .25, left = .15) 394 calib.invT_xaxis(Ti = [0,10,25]) 395 calib.plot_bff_ci(alpha = .15) 396 calib.plot_bff(label = True, dashes = (8,2,2,2)) 397 calib.plot_data() 398 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 399 ppl.legend() 400 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 401 ````` 402 403 This should result in something like this: 404 405 <img align="center" src="example_plot_bff_ci.png"> 406 """ 407# if 'color' not in kwargs: 408# kwargs['color'] = self.color 409 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs) 410 411 def T47(self, 412 D47 = None, 413 sD47 = None, 414 T=None, 415 sT = None, 416 error_from = 'both', 417 return_covar = False, 418 ): 419 ''' 420 When `D47` is input, computes corresponding T value(s). 421 `D47` input may be specified as a scalar, or as a 1-D array. 422 `T` output will then have the same type and size as `D47`. 423 424 When `T` is input, computes corresponding Δ47 value(s). 425 `T` input may be specified as a scalar, or as a 1-D array. 426 `D47` output will then have the same type and size as `T`. 427 428 Only one of either `D47` or `T` may be specified as input. 429 430 **Arguments:** 431 432 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 433 * `sD47`: Δ47 uncertainties, which may be: 434 - `None` (default) 435 - `float` or `int` (uniform standard error on `D47`) 436 - 1-D array (standard errors on `D47`) 437 - 2-D array (covariance matrix for `D47`) 438 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 439 * `sT`: T uncertainties, which may be: 440 - `None` (default) 441 - `float` or `int` (uniform standard error on `T`) 442 - 1-D array (standard errors on `T`) 443 - 2-D array (variance-covariance matrix for `T`) 444 * `error_from`: if set to `'both'` (default), returned errors take into account 445 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 446 if set to `'calib'`, only calibration uncertainties are accounted for; 447 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 448 * `return_covar`: (False by default) whether to return the full covariance matrix 449 for returned `T` or `D47` values, otherwise return standard errors for the returned 450 `T` or `D47` values instead. 451 452 **Returns (with `D47` input):** 453 454 * `T`: temperature value(s) computed from `D47` 455 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 456 457 **Returns (with `T` input):** 458 459 * `D47`: Δ47 value(s) computed from `D47` 460 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 461 462 ### Example 463 464 ````py 465 import numpy as np 466 from matplotlib import pyplot as ppl 467 from D47calib import OGLS23 as calib 468 469 X = np.linspace(1473**-2, 270**-2) 470 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 471 472 fig = ppl.figure(figsize = (5,3)) 473 ppl.subplots_adjust(bottom = .25, left = .15) 474 calib.invT_xaxis() 475 ppl.plot(X, 1000 * sD47, 'r-') 476 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 477 ppl.savefig('example_SE47.png', dpi = 100) 478 ````` 479 480 This should result in something like this: 481 482 <img src="example_SE47.png"> 483 ''' 484 485 if D47 is None and T is None: 486 raise ValueError('Either D47 or T must be specified, but both are undefined.') 487 488 if D47 is not None and T is not None: 489 raise ValueError('Either D47 or T must be specified, but not both.') 490 491 if T is not None: 492 493 D47 = self._D47_from_T(T) 494 Np = len(self.degrees) 495 N = D47.size 496 497 ### Compute covariance matrix of (*bfp, *T): 498 CM = _np.zeros((Np+N, Np+N)) 499 500 if error_from in ['calib', 'both']: 501 CM[:Np, :Np] = self.bfp_CM[:,:] 502 503 if (sT is not None) and error_from in ['sT', 'both']: 504 _sT = _np.asarray(sT) 505 if _sT.ndim == 0: 506 for k in range(N): 507 CM[Np+k, Np+k] = _sT**2 508 elif _sT.ndim == 1: 509 for k in range(N): 510 CM[Np+k, Np+k] = _sT[k]**2 511 elif _sT.ndim == 2: 512 CM[-N:, -N:] = _sT[:,:] 513 514 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 515 _T = _np.asarray(T) 516 if _T.ndim == 0: 517 _T = _np.expand_dims(_T, 0) 518 J = _np.zeros((N, Np+N)) 519 520 if (sT is not None) and error_from in ['sT', 'both']: 521 for k in range(N): 522 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 523 524 if error_from in ['calib', 'both']: 525 526 for k in range(Np): 527 528 p1 = {_: self.bfp[_] for _ in self.bfp} 529 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 530 531 p2 = {_: self.bfp[_] for _ in self.bfp} 532 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 533 534 J[:, k] = (self.model_fun(p1, (_T+273.15)**-1) - self.model_fun(p2, (_T+273.15)**-1)) / (0.002 * self.bfp_CM[k,k]**.5) 535 536 ### Error propagation: 537 CM_D47 = J @ CM @ J.T 538 539 if return_covar: 540 return D47, CM_D47 541 else: 542 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 543 544 if D47 is not None: 545 546 T = self._T_from_D47(D47) 547 Np = len(self.degrees) 548 N = T.size 549 550 ### Compute covariance matrix of (*bfp, *T): 551 CM = _np.zeros((Np+N, Np+N)) 552 553 if error_from in ['calib', 'both']: 554 CM[:Np, :Np] = self.bfp_CM[:,:] 555 556 if (sD47 is not None) and error_from in ['sD47', 'both']: 557 _sD47 = _np.asarray(sD47) 558 if _sD47.ndim == 0: 559 for k in range(N): 560 CM[Np+k, Np+k] = _sD47**2 561 elif _sD47.ndim == 1: 562 for k in range(N): 563 CM[Np+k, Np+k] = _sD47[k]**2 564 elif _sD47.ndim == 2: 565 CM[-N:, -N:] = _sD47[:,:] 566 567 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 568 _D47 = _np.asarray(D47) 569 if _D47.ndim == 0: 570 _D47 = _np.expand_dims(_D47, 0) 571 J = _np.zeros((N, Np+N)) 572 if (sD47 is not None) and error_from in ['sD47', 'both']: 573 for k in range(N): 574 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 575 if error_from in ['calib', 'both']: 576 577 xi = _np.linspace(0,200**-1,1001)[1:] 578 for k in range(Np): 579 580 if self.bfp_CM[k,k]: 581 _epsilon_ = self.bfp_CM[k,k]**.5 582 else: 583 _epsilon_ = 1e-6 584 585 p1 = {_: self.bfp[_] for _ in self.bfp} 586 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 587 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 588 589 p2 = {_: self.bfp[_] for _ in self.bfp} 590 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 591 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 592 593 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 594 595 ### Error propagation: 596 CM_T = J @ CM @ J.T 597 598 if return_covar: 599 return T, CM_T 600 else: 601 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5 602 603 604 def plot_T47_errors( 605 self, 606 calibname = None, 607 rD47 = 0.010, 608 Nr = [2,4,8,12,20], 609 Tmin = 0, 610 Tmax = 120, 611 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 612 yscale = 'lin', 613 ): 614 """ 615 Plot SE of T reconstructed using the calibration as a function of T for various 616 combinations of analytical precision and number of analytical replicates. 617 618 **Arguments** 619 620 + **calibname**: 621 Which calibration name to display. By default, use `label` attribute. 622 + **rD47**: 623 Analytical precision of a single analysis. 624 + **Nr**: 625 A list of lines to plot, each corresponding to a given number of replicates. 626 + **Tmin**: 627 Minimum T to plot. 628 + **Tmax**: 629 Maximum T to plot. 630 + **colors**: 631 A list of colors to distinguish the plotted lines. 632 + **yscale**: 633 + If `'lin'`, the Y axis uses a linear scale. 634 + If `'log'`, the Y axis uses a logarithmic scale. 635 636 **Example** 637 638 ````py 639 from matplotlib import pyplot as ppl 640 from D47calib import devils_laghetto_2023 as calib 641 642 fig = ppl.figure(figsize = (3.5,4)) 643 ppl.subplots_adjust(bottom = .2, left = .15) 644 calib.plot_T47_errors( 645 calibname = 'Devils Laghetto calibration', 646 Nr = [1,2,4,16], 647 Tmin =0, 648 Tmax = 40, 649 ) 650 ppl.savefig('example_SE_T.png', dpi = 100) 651 ```` 652 653 This should result in something like this: 654 655 <img src="example_SE_T.png"> 656 """ 657 658 if calibname is None: 659 calibname = self.label 660 661 Nr = _np.array(Nr) 662 if len(colors) < Nr.size: 663 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 664 from colorsys import hsv_to_rgb 665 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 666 colors = [hsv_to_rgb(*x) for x in hsv] 667 668 Ti = _np.linspace(Tmin, Tmax) 669 D47i, _ = self.T47(T = Ti) 670 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 671 672 ymax, ymin = 0, 1e6 673 for N,c in zip(Nr, colors): 674 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 675 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 676 ymin = min(ymin, min(sT)) 677 ymax = max(ymax, max(sT)) 678 679 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 680 681 _ppl.legend(fontsize=9) 682 _ppl.xlabel("T (°C)") 683 684 _ppl.ylabel("Standard error on reconstructed T (°C)") 685 686 # yticks([0,.5,1,1.5,2]) 687 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 688 _ppl.grid( alpha = .25) 689 if yscale == 'lin': 690 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 691 t1, t2 = self.T.min(), self.T.max() 692 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 693 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 694 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 695 elif yscale == 'log': 696 ymin /= 2 697 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 698 _ppl.yscale('log') 699 t1, t2 = self.T.min(), self.T.max() 700 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 701 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 702 703 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 704 """ 705 Write calibration data to a csv file. 706 707 ### Parameters 708 709 + **csvfile**: 710 The filename to write data to. 711 + **sep**: 712 The separator between CSV fields. 713 + **label**: 714 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 715 + If specified as a `str`, include a `Dataset` column with that string. 716 + If specified as `False`, do not include a `Dataset` column. 717 + **T_correl**: 718 + If `True`, include correlations between all `T` values. 719 + **D47_correl**: 720 + If `True`, include correlations between all `D47` values. 721 722 ### Example 723 724 ````py 725 D47calib.huyghe_2022.export_data( 726 csvfile = 'example_export_data.csv', 727 T_correl = True, 728 D47_correl = True, 729 ) 730 ```` 731 732 This should result in something like this ([link](example_export_data.csv)): 733 734 .. include:: ../../docs/example_export_data.md 735 736 """ 737 n = len(str(self.N)) 738 739 with open(csvfile, 'w') as f: 740 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 741 742 if label: 743 f.write(f'{sep}Dataset') 744 745 if T_correl: 746 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 747 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 748 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 749 750 if D47_correl: 751 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 752 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 753 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 754 755 for k, (s, T, sT, D47, sD47) in enumerate(zip( 756 self.samples, 757 self.T, 758 _np.diag(self.sT)**.5, 759 self.D47, 760 _np.diag(self.sD47)**.5, 761 )): 762 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 763 if label: 764 if label is True: 765 f.write(f'{sep}{self.label}') 766 else: 767 f.write(f'{sep}{label}') 768 if T_correl: 769 f.write(sep.join(['']+[ 770 f'{Tcorrel[k,_]:.0f}' 771 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 772 else f'{Tcorrel[k,_]:.6f}' 773 for _ in range(self.N)])) 774 if D47_correl: 775 f.write(sep.join(['']+[ 776 f'{D47correl[k,_]:.0f}' 777 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 778 else f'{D47correl[k,_]:.6f}' 779 for _ in range(self.N)])) 780 781 782 def export(self, name, filename): 783 """ 784 Save `D47calib` object as an importable file. 785 786 ### Parameters 787 788 + **name**: 789 The name of the variable to export. 790 + **filename**: 791 The filename to write to. 792 793 ### Example 794 795 ````py 796 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 797 ```` 798 799 This should result in a `bar.py` file with the following contents: 800 801 ````py 802 foo = D47calib( 803 samples = ['LGB-2', 'DVH-2'], 804 T = [7.9, 33.7], 805 D47 = [0.6485720997671647, 0.5695972909966959], 806 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 807 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 808 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 809 label = 'Slow-growing calcites from Anderson et al. (2021)', 810 color = (0, 0.5, 0), 811 degrees = [0, 2], 812 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 813 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 814 chisq = 6.421311854486162e-27, 815 Nf = 0, 816 ) 817 ```` 818 """ 819 with open(filename, 'w') as f: 820 f.write(f''' 821{name} = D47calib( 822 samples = {self.samples}, 823 T = {list(self.T)}, 824 D47 = {list(self.D47)}, 825 sT = {[list(l) for l in self.sT]}, 826 sD47 = {[list(l) for l in self.sD47]}, 827 degrees = {self.degrees}, 828 description = {repr(self.description)}, 829 name = {repr(self.name)}, 830 label = {repr(self.label)}, 831 bfp = {self.bfp}, 832 bfp_CM = {[list(l) for l in self.bfp_CM]}, 833 chisq = {self.chisq}, 834 cholesky_residuals = {list(self.cholesky_residuals)}, 835 aic = {self.aic}, 836 bic = {self.bic}, 837 ks_pvalue = {self.ks_pvalue}, 838 ) 839''')
Δ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)
46 def __init__(self, 47 samples, T, D47, 48 sT = None, 49 sD47 = None, 50 degrees = [0,2], 51 xpower = 2, 52 name = '', 53 label = '', 54 description = '', 55 **kwargs, 56 ): 57 """ 58 ### Parameters 59 60 + **samples**: a list of N sample names. 61 + **T**: a 1-D array (or array-like) of temperatures values (in degrees C), of size N. 62 + **D47**: a 1-D array (or array-like) of Δ47 values (in permil), of size N. 63 + **sT**: uncertainties on `T`. If specified as: 64 + a scalar: `sT` is treated as the standard error applicable to all `T` values; 65 + a 1-D array-like of size N: `sT` is treated as the standard errors of `T`; 66 + a 2-D array-like of size (N, N): `sT` is treated as the (co)variance matrix of `T`. 67 + **sD47**: uncertainties on `D47`. If specified as: 68 + a scalar: `sD47` is treated as the standard error applicable to all `D47` values; 69 + a 1-D array-like of size N: `sD47` is treated as the standard errors of `D47`; 70 + a 2-D array-like of size (N, N): `sD47` is treated as the (co)variance matrix of `D47`. 71 + **degrees**: degrees of the polynomial regression, e.g., `[0, 2]` or `[0, 1, 2, 3, 4]`. 72 + **name**: a human-readable, short name assigned to the calibration. 73 + **label**: a short description of the calibration, e.g., to be used in legends. 74 + **description**: a longer description, including relevant references/DOIs. 75 This is not necessary when `bfp` and `CM_bfp` are specified at instantiation time. 76 + **kwargs**: keyword arguments passed to the underlying `ogls.InverseTPolynomial()` call. 77 78 ### Notable attributes 79 80 + **N**: 81 The total number of observations (samples) in the calibration data. 82 + **samples**: 83 The list sample names. 84 + **T**: 85 1-D `ndarray` of temperatures in degrees C. 86 + **D47**: 87 1-D `ndarray` of Δ47 values in permil. 88 + **sT**: 89 2-D `ndarray` equal to the full (co)variance matrix for `T`. 90 + **D47**: 91 2-D `ndarray` equal to the full (co)variance matrix for `D47`. 92 + **xpower**: 93 By default, all `D47calib` graphical methods plot Δ47 as a function of 1/T<sup>2</sup>. 94 It is possible to change this behavior to use a different power of 1/T. 95 This is done by redefining the `xpower` attribute to a different, non-zero `int` value 96 (e.g. `foo.xpower = 1` to plot as a function of 1/T instead of 1/T<sup>2</sup>). 97 + **bfp**: 98 The best-fit parameters of the regression. 99 This is a `dict` with keys equal to the polynomial coefficients (see `bff` definition below) 100 + **bff()**: 101 The best-fit polynomial function of inverse T, defined as: 102 `bff(x) = sum(bfp[f'a{k}'] * x**k for k in degrees)` 103 Note that `bff` takes `x = 1/(T+273.15)` (instead of `T`) as input. 104 105 106 ### Examples 107 108 A very simple example: 109 110 ````py 111 .. include:: ../../code_examples/D47calib_init/example.py 112 ```` 113 114 Should yield: 115 116 ```` 117 .. include:: ../../code_examples/D47calib_init/output.txt 118 ```` 119 120 """ 121 122 self.samples = samples[:] 123 self.name = name 124 self.label = label 125 self.description = description 126 self.D47 = _np.asarray(D47, dtype = 'float') 127 self.N = self.D47.size 128 129 if sD47 is None: 130 self.sD47 = _np.zeros((self.N, self.N)) 131 else: 132 self.sD47 = _np.asarray(sD47) 133 if len(self.sD47.shape) == 1: 134 self.sD47 = _np.diag(self.sD47**2) 135 elif len(self.sD47.shape) == 0: 136 self.sD47 = _np.eye(self.D47.size) * self.sD47**2 137 138 _ogls.InverseTPolynomial.__init__(self, T=T, Y=D47, sT=sT, sY=sD47, degrees = degrees, xpower = xpower, **kwargs) 139 140 if self.bfp is None: 141 self.regress() 142 143 self._bff_deriv = lambda x: _np.array([k * self.bfp[f'a{k}'] * x**(k-1) for k in degrees if k > 0]).sum(axis = 0) 144 145 xi = _np.linspace(0,200**-1,1001) 146 self._inv_bff = _interp1d(self.bff(xi), xi) 147 148 self._D47_from_T = lambda T: self.bff((T+273.15)**-1) 149 self._T_from_D47 = lambda D47: self._inv_bff(D47)**-1 - 273.15 150 self._D47_from_T_deriv = lambda T: -(T+273.15)**-2 * self._bff_deriv((T+273.15)**-1) 151 self._T_from_D47_deriv = lambda D47: self._D47_from_T_deriv(self._T_from_D47(D47))**-1
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:- a scalar:
sTis treated as the standard error applicable to allTvalues; - a 1-D array-like of size N:
sTis treated as the standard errors ofT; - a 2-D array-like of size (N, N):
sTis treated as the (co)variance matrix ofT.
- a scalar:
- sD47: uncertainties on
D47. If specified as:- a scalar:
sD47is treated as the standard error applicable to allD47values; - a 1-D array-like of size N:
sD47is treated as the standard errors ofD47; - a 2-D array-like of size (N, N):
sD47is treated as the (co)variance matrix ofD47.
- a scalar:
- 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]):
156 def invT_xaxis(self, 157 xlabel = None, 158 Ti = [0,20,50,100,250,1000], 159 ): 160 """ 161 Create and return an `Axes` object with X values equal to 1/T<sup>2</sup>, 162 but labeled in degrees Celsius. 163 164 ### Parameters 165 166 + **xlabel**: 167 Custom label for X axis (`r'$1\,/\,T^2$'` by default) 168 + **Ti**: 169 Specify tick locations for X axis, in degrees C. 170 171 ### Returns 172 173 + an `matplotlib.axes.Axes` instance 174 175 ### Examples 176 177 ````py 178 .. include:: ../../code_examples/D47calib_invT_xaxis/example_1.py 179 ```` 180 181 This should result in something like this: 182 183 <img align="center" src="example_invT_xaxis_1.png"> 184 185 It is also possible to define the X axis using a different power of 1/T 186 by first redefining the `xpower` attribute: 187 188 ````py 189 .. include:: ../../code_examples/D47calib_invT_xaxis/example_2.py 190 ```` 191 192 This should result in something like this: 193 194 <img align="center" src="example_invT_xaxis_2.png"> 195 """ 196 if xlabel is None: 197 xlabel = f'$1\\,/\\,T^{self.xpower}$' if self.xpower > 1 else '1/T' 198 _ppl.xlabel(xlabel) 199 _ppl.xticks([(273.15 + t) ** -self.xpower for t in sorted(Ti)[::-1]]) 200 ax = _ppl.gca() 201 ax.set_xticklabels([f"${t}\\,$°C" for t in sorted(Ti)[::-1]]) 202 ax.tick_params(which="major") 203 204 return ax
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):
207 def plot_data(self, label = False, **kwargs): 208 """ 209 Plot Δ47 value of each sample as a function of 1/T<sup>2</sup>. 210 211 ### Parameters 212 213 + **label**: 214 + If `label` is a string, use this string as `label` for the underlyig 215 `matplotlib.pyplot.plot()` call. 216 + If `label = True`, use the caller's `label` attribute instead. 217 + If `label = False`, no label is specified (default behavior). 218 + **kwargs**: 219 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 220 221 ### Returns 222 223 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 224 225 ### Example 226 227 ````py 228 from matplotlib import pyplot as ppl 229 from D47calib import huyghe_2022 as calib 230 231 fig = ppl.figure(figsize = (5,3)) 232 ppl.subplots_adjust(bottom = .25, left = .15) 233 calib.invT_xaxis(Ti = [0,10,25]) 234 calib.plot_data(label = True) 235 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 236 ppl.legend() 237 ppl.savefig('example_plot_data.png', dpi = 100) 238 ````` 239 240 This should result in something like this: 241 242 <img align="center" src="example_plot_data.png"> 243 """ 244# if 'mec' not in kwargs: 245# kwargs['mec'] = self.color 246 if label is not False: 247 kwargs['label'] = self.label if label is True else label 248 return _ogls.InverseTPolynomial.plot_data(self, **kwargs)
Plot Δ47 value of each sample as a function of 1/T2.
Parameters
- label:
- If
labelis a string, use this string aslabelfor the underlyigmatplotlib.pyplot.plot()call. - If
label = True, use the caller'slabelattribute instead. - If
label = False, no label is specified (default behavior).
- If
- 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):
251 def plot_error_bars(self, **kwargs): 252 """ 253 Plot Δ47 error bars (±1.96 SE) of each sample as a function of 1/T<sup>2</sup>. 254 255 ### Parameters 256 257 + **kwargs**: 258 keyword arguments passed to the underlying `matplotlib.pyplot.errrobar()` call. 259 260 ### Returns 261 262 + the return value(s) of the underlying `matplotlib.pyplot.errorbar()` call. 263 264 ### Example 265 266 ````py 267 from matplotlib import pyplot as ppl 268 from D47calib import huyghe_2022 as calib 269 270 fig = ppl.figure(figsize = (5,3)) 271 ppl.subplots_adjust(bottom = .25, left = .15) 272 calib.invT_xaxis(Ti = [0,10,25]) 273 calib.plot_error_bars(alpha = .4) 274 calib.plot_data(label = True) 275 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 276 ppl.legend() 277 ppl.savefig('example_plot_error_bars.png', dpi = 100) 278 ````` 279 280 This should result in something like this: 281 282 <img align="center" src="example_plot_error_bars.png"> 283 """ 284# if 'ecolor' not in kwargs: 285# kwargs['ecolor'] = self.color 286 return _ogls.InverseTPolynomial.plot_error_bars(self, **kwargs)
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):
289 def plot_error_ellipses(self, **kwargs): 290 """ 291 Plot Δ47 error ellipses (95 % confidence) of each sample as a function of 1/T<sup>2</sup>. 292 293 ### Parameters 294 295 + **kwargs**: 296 keyword arguments passed to the underlying `matplotlib.patches.Ellipse()` call. 297 298 ### Returns 299 300 + the return value(s) of the underlying `matplotlib.patches.Ellipse()` call. 301 302 ### Example 303 304 ````py 305 from matplotlib import pyplot as ppl 306 from D47calib import huyghe_2022 as calib 307 308 fig = ppl.figure(figsize = (5,3)) 309 ppl.subplots_adjust(bottom = .25, left = .15) 310 calib.invT_xaxis(Ti = [0,10,25]) 311 calib.plot_error_ellipses(alpha = .4) 312 calib.plot_data(label = True) 313 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 314 ppl.legend() 315 ppl.savefig('example_plot_error_ellipses.png', dpi = 100) 316 ````` 317 318 This should result in something like this: 319 320 <img align="center" src="example_plot_error_ellipses.png"> 321 """ 322# if 'ec' not in kwargs: 323# kwargs['ec'] = self.color 324 return _ogls.InverseTPolynomial.plot_error_ellipses(self, **kwargs)
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):
327 def plot_bff(self, label = False, **kwargs): 328 """ 329 Plot best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 330 331 ### Parameters 332 333 + **label**: 334 + If `label` is a string, use this string as `label` for the underlyig 335 `matplotlib.pyplot.plot()` call. 336 + If `label = True`, use the caller's `label` attribute instead. 337 + If `label = False`, no label is specified (default behavior). 338 + **kwargs**: 339 keyword arguments passed to the underlying `matplotlib.pyplot.plot()` call. 340 341 ### Returns 342 343 + the return value(s) of the underlying `matplotlib.pyplot.plot()` call. 344 345 ### Example 346 347 ````py 348 from matplotlib import pyplot as ppl 349 from D47calib import huyghe_2022 as calib 350 351 fig = ppl.figure(figsize = (5,3)) 352 ppl.subplots_adjust(bottom = .25, left = .15) 353 calib.invT_xaxis(Ti = [0,10,25]) 354 calib.plot_bff(label = True, dashes = (8,2,2,2)) 355 calib.plot_data() 356 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 357 ppl.legend() 358 ppl.savefig('example_plot_bff.png', dpi = 100) 359 ````` 360 361 This should result in something like this: 362 363 <img align="center" src="example_plot_bff.png"> 364 """ 365# if 'color' not in kwargs: 366# kwargs['color'] = self.color 367 if label is not False: 368 kwargs['label'] = self.label if label is True else label 369 return _ogls.InverseTPolynomial.plot_bff(self, **kwargs)
Plot best-fit regression of Δ47 as a function of 1/T2.
Parameters
- label:
- If
labelis a string, use this string aslabelfor the underlyigmatplotlib.pyplot.plot()call. - If
label = True, use the caller'slabelattribute instead. - If
label = False, no label is specified (default behavior).
- If
- 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):
372 def plot_bff_ci(self, **kwargs): 373 """ 374 Plot 95 % confidence region for best-fit regression of Δ47 as a function of 1/T<sup>2</sup>. 375 376 ### Parameters 377 378 + **label**: 379 + **kwargs**: 380 keyword arguments passed to the underlying `matplotlib.pyplot.fill_between()` call. 381 382 ### Returns 383 384 + the return value(s) of the underlying `matplotlib.pyplot.fill_between()` call. 385 386 ### Example 387 388 ````py 389 from matplotlib import pyplot as ppl 390 from D47calib import huyghe_2022 as calib 391 392 fig = ppl.figure(figsize = (5,3)) 393 ppl.subplots_adjust(bottom = .25, left = .15) 394 calib.invT_xaxis(Ti = [0,10,25]) 395 calib.plot_bff_ci(alpha = .15) 396 calib.plot_bff(label = True, dashes = (8,2,2,2)) 397 calib.plot_data() 398 ppl.ylabel('$Δ_{47}$ (‰ I-CDES)') 399 ppl.legend() 400 ppl.savefig('example_plot_bff_ci.png', dpi = 100) 401 ````` 402 403 This should result in something like this: 404 405 <img align="center" src="example_plot_bff_ci.png"> 406 """ 407# if 'color' not in kwargs: 408# kwargs['color'] = self.color 409 return _ogls.InverseTPolynomial.plot_bff_ci(self, **kwargs)
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):
411 def T47(self, 412 D47 = None, 413 sD47 = None, 414 T=None, 415 sT = None, 416 error_from = 'both', 417 return_covar = False, 418 ): 419 ''' 420 When `D47` is input, computes corresponding T value(s). 421 `D47` input may be specified as a scalar, or as a 1-D array. 422 `T` output will then have the same type and size as `D47`. 423 424 When `T` is input, computes corresponding Δ47 value(s). 425 `T` input may be specified as a scalar, or as a 1-D array. 426 `D47` output will then have the same type and size as `T`. 427 428 Only one of either `D47` or `T` may be specified as input. 429 430 **Arguments:** 431 432 * `D47`: Δ47 value(s) to convert into temperature (`float` or 1-D array) 433 * `sD47`: Δ47 uncertainties, which may be: 434 - `None` (default) 435 - `float` or `int` (uniform standard error on `D47`) 436 - 1-D array (standard errors on `D47`) 437 - 2-D array (covariance matrix for `D47`) 438 * `T`: T value(s) to convert into Δ47 (`float` or 1-D array), in degrees C 439 * `sT`: T uncertainties, which may be: 440 - `None` (default) 441 - `float` or `int` (uniform standard error on `T`) 442 - 1-D array (standard errors on `T`) 443 - 2-D array (variance-covariance matrix for `T`) 444 * `error_from`: if set to `'both'` (default), returned errors take into account 445 input uncertainties (`sT` or `sD47`) as well as calibration uncertainties; 446 if set to `'calib'`, only calibration uncertainties are accounted for; 447 if set to `'sT'` or `'sD47'`, calibration uncertainties are ignored. 448 * `return_covar`: (False by default) whether to return the full covariance matrix 449 for returned `T` or `D47` values, otherwise return standard errors for the returned 450 `T` or `D47` values instead. 451 452 **Returns (with `D47` input):** 453 454 * `T`: temperature value(s) computed from `D47` 455 * `sT`: uncertainties on `T` value(s), whether as standard error(s) or covariance matrix 456 457 **Returns (with `T` input):** 458 459 * `D47`: Δ47 value(s) computed from `D47` 460 * `sD47`: uncertainties on `D47` value(s), whether as standard error(s) or covariance matrix 461 462 ### Example 463 464 ````py 465 import numpy as np 466 from matplotlib import pyplot as ppl 467 from D47calib import OGLS23 as calib 468 469 X = np.linspace(1473**-2, 270**-2) 470 D47, sD47 = calib.T47(T = X**-0.5 - 273.15) 471 472 fig = ppl.figure(figsize = (5,3)) 473 ppl.subplots_adjust(bottom = .25, left = .15) 474 calib.invT_xaxis() 475 ppl.plot(X, 1000 * sD47, 'r-') 476 ppl.ylabel('Calibration SE on $Δ_{47}$ values (ppm)') 477 ppl.savefig('example_SE47.png', dpi = 100) 478 ````` 479 480 This should result in something like this: 481 482 <img src="example_SE47.png"> 483 ''' 484 485 if D47 is None and T is None: 486 raise ValueError('Either D47 or T must be specified, but both are undefined.') 487 488 if D47 is not None and T is not None: 489 raise ValueError('Either D47 or T must be specified, but not both.') 490 491 if T is not None: 492 493 D47 = self._D47_from_T(T) 494 Np = len(self.degrees) 495 N = D47.size 496 497 ### Compute covariance matrix of (*bfp, *T): 498 CM = _np.zeros((Np+N, Np+N)) 499 500 if error_from in ['calib', 'both']: 501 CM[:Np, :Np] = self.bfp_CM[:,:] 502 503 if (sT is not None) and error_from in ['sT', 'both']: 504 _sT = _np.asarray(sT) 505 if _sT.ndim == 0: 506 for k in range(N): 507 CM[Np+k, Np+k] = _sT**2 508 elif _sT.ndim == 1: 509 for k in range(N): 510 CM[Np+k, Np+k] = _sT[k]**2 511 elif _sT.ndim == 2: 512 CM[-N:, -N:] = _sT[:,:] 513 514 ### Compute Jacobian of D47(T) relative to (*bfp, *T): 515 _T = _np.asarray(T) 516 if _T.ndim == 0: 517 _T = _np.expand_dims(_T, 0) 518 J = _np.zeros((N, Np+N)) 519 520 if (sT is not None) and error_from in ['sT', 'both']: 521 for k in range(N): 522 J[k, Np+k] = self._D47_from_T_deriv(_T[k]) 523 524 if error_from in ['calib', 'both']: 525 526 for k in range(Np): 527 528 p1 = {_: self.bfp[_] for _ in self.bfp} 529 p1[f'a{self.degrees[k]}'] += 0.001 * self.bfp_CM[k,k]**.5 530 531 p2 = {_: self.bfp[_] for _ in self.bfp} 532 p2[f'a{self.degrees[k]}'] -= 0.001 * self.bfp_CM[k,k]**.5 533 534 J[:, k] = (self.model_fun(p1, (_T+273.15)**-1) - self.model_fun(p2, (_T+273.15)**-1)) / (0.002 * self.bfp_CM[k,k]**.5) 535 536 ### Error propagation: 537 CM_D47 = J @ CM @ J.T 538 539 if return_covar: 540 return D47, CM_D47 541 else: 542 return D47, float(_np.diag(CM_D47)**.5) if D47.ndim == 0 else _np.diag(CM_D47)**.5 543 544 if D47 is not None: 545 546 T = self._T_from_D47(D47) 547 Np = len(self.degrees) 548 N = T.size 549 550 ### Compute covariance matrix of (*bfp, *T): 551 CM = _np.zeros((Np+N, Np+N)) 552 553 if error_from in ['calib', 'both']: 554 CM[:Np, :Np] = self.bfp_CM[:,:] 555 556 if (sD47 is not None) and error_from in ['sD47', 'both']: 557 _sD47 = _np.asarray(sD47) 558 if _sD47.ndim == 0: 559 for k in range(N): 560 CM[Np+k, Np+k] = _sD47**2 561 elif _sD47.ndim == 1: 562 for k in range(N): 563 CM[Np+k, Np+k] = _sD47[k]**2 564 elif _sD47.ndim == 2: 565 CM[-N:, -N:] = _sD47[:,:] 566 567 ### Compute Jacobian of T(D47) relative to (*bfp, *D47): 568 _D47 = _np.asarray(D47) 569 if _D47.ndim == 0: 570 _D47 = _np.expand_dims(_D47, 0) 571 J = _np.zeros((N, Np+N)) 572 if (sD47 is not None) and error_from in ['sD47', 'both']: 573 for k in range(N): 574 J[k, Np+k] = self._T_from_D47_deriv(_D47[k]) 575 if error_from in ['calib', 'both']: 576 577 xi = _np.linspace(0,200**-1,1001)[1:] 578 for k in range(Np): 579 580 if self.bfp_CM[k,k]: 581 _epsilon_ = self.bfp_CM[k,k]**.5 582 else: 583 _epsilon_ = 1e-6 584 585 p1 = {_: self.bfp[_] for _ in self.bfp} 586 p1[f'a{self.degrees[k]}'] += 0.001 * _epsilon_ 587 T_from_D47_p1 = _interp1d(self.model_fun(p1, xi), xi**-1 - 273.15) 588 589 p2 = {_: self.bfp[_] for _ in self.bfp} 590 p2[f'a{self.degrees[k]}'] -= 0.001 * _epsilon_ 591 T_from_D47_p2 = _interp1d(self.model_fun(p2, xi), xi**-1 - 273.15) 592 593 J[:, k] = (T_from_D47_p1(_D47) - T_from_D47_p2(_D47)) / (0.002 * _epsilon_) 594 595 ### Error propagation: 596 CM_T = J @ CM @ J.T 597 598 if return_covar: 599 return T, CM_T 600 else: 601 return T, float(_np.diag(CM_T)**.5) if T.ndim == 0 else _np.diag(CM_T)**.5
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:None(default)floatorint(uniform standard error onD47)- 1-D array (standard errors on
D47) - 2-D array (covariance matrix for
D47)
T: T value(s) to convert into Δ47 (floator 1-D array), in degrees CsT: T uncertainties, which may be:None(default)floatorint(uniform standard error onT)- 1-D array (standard errors on
T) - 2-D array (variance-covariance matrix for
T)
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'):
604 def plot_T47_errors( 605 self, 606 calibname = None, 607 rD47 = 0.010, 608 Nr = [2,4,8,12,20], 609 Tmin = 0, 610 Tmax = 120, 611 colors = [(1,0,0),(1,.5,0),(.25,.75,0),(0,.5,1),(0.5,0.5,0.5)], 612 yscale = 'lin', 613 ): 614 """ 615 Plot SE of T reconstructed using the calibration as a function of T for various 616 combinations of analytical precision and number of analytical replicates. 617 618 **Arguments** 619 620 + **calibname**: 621 Which calibration name to display. By default, use `label` attribute. 622 + **rD47**: 623 Analytical precision of a single analysis. 624 + **Nr**: 625 A list of lines to plot, each corresponding to a given number of replicates. 626 + **Tmin**: 627 Minimum T to plot. 628 + **Tmax**: 629 Maximum T to plot. 630 + **colors**: 631 A list of colors to distinguish the plotted lines. 632 + **yscale**: 633 + If `'lin'`, the Y axis uses a linear scale. 634 + If `'log'`, the Y axis uses a logarithmic scale. 635 636 **Example** 637 638 ````py 639 from matplotlib import pyplot as ppl 640 from D47calib import devils_laghetto_2023 as calib 641 642 fig = ppl.figure(figsize = (3.5,4)) 643 ppl.subplots_adjust(bottom = .2, left = .15) 644 calib.plot_T47_errors( 645 calibname = 'Devils Laghetto calibration', 646 Nr = [1,2,4,16], 647 Tmin =0, 648 Tmax = 40, 649 ) 650 ppl.savefig('example_SE_T.png', dpi = 100) 651 ```` 652 653 This should result in something like this: 654 655 <img src="example_SE_T.png"> 656 """ 657 658 if calibname is None: 659 calibname = self.label 660 661 Nr = _np.array(Nr) 662 if len(colors) < Nr.size: 663 print('WARNING: Too few colors to plot different numbers of replicates; generating new colors.') 664 from colorsys import hsv_to_rgb 665 hsv = [(x*1.0/Nr.size, 1, .9) for x in range(Nr.size)] 666 colors = [hsv_to_rgb(*x) for x in hsv] 667 668 Ti = _np.linspace(Tmin, Tmax) 669 D47i, _ = self.T47(T = Ti) 670 _, sT_calib = self.T47(D47 = D47i, error_from = 'calib') 671 672 ymax, ymin = 0, 1e6 673 for N,c in zip(Nr, colors): 674 _, sT = self.T47(D47 = D47i, sD47 = rD47 / N**.5, error_from = 'sD47') 675 _ppl.plot(Ti, sT, '-', color = c, label=f'SE for {N} replicate{"s" if N > 1 else ""}') 676 ymin = min(ymin, min(sT)) 677 ymax = max(ymax, max(sT)) 678 679 _ppl.plot(Ti, sT_calib, 'k--', label='SE from calibration') 680 681 _ppl.legend(fontsize=9) 682 _ppl.xlabel("T (°C)") 683 684 _ppl.ylabel("Standard error on reconstructed T (°C)") 685 686 # yticks([0,.5,1,1.5,2]) 687 _ppl.title(f"{calibname},\nassuming external Δ$_{{47}}$ repeatability of {rD47:.3f} ‰", size = 9) 688 _ppl.grid( alpha = .25) 689 if yscale == 'lin': 690 _ppl.axis([Ti[0], Ti[-1], 0, ymax*1.05]) 691 t1, t2 = self.T.min(), self.T.max() 692 _ppl.plot([t1, t2], [0, 0], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 693 _ppl.text((t1+t2)/2, 0, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic') 694 _ppl.axis([None, None, None, _ppl.axis()[-1]*1.25]) 695 elif yscale == 'log': 696 ymin /= 2 697 _ppl.axis([Ti[0], Ti[-1], ymin, ymax*1.05]) 698 _ppl.yscale('log') 699 t1, t2 = self.T.min(), self.T.max() 700 _ppl.plot([t1, t2], [ymin, ymin], 'k-', alpha = .25, lw = 8, solid_capstyle = 'butt', clip_on = False) 701 _ppl.text((t1+t2)/2, ymin, 'range of observations\n', alpha = .4, size = 7, ha = 'center', va = 'bottom', style = 'italic')
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):
703 def export_data(self, csvfile, sep = ',', label = False, T_correl = False, D47_correl = False): 704 """ 705 Write calibration data to a csv file. 706 707 ### Parameters 708 709 + **csvfile**: 710 The filename to write data to. 711 + **sep**: 712 The separator between CSV fields. 713 + **label**: 714 + If specified as `True`, include a `Dataset` column with the calibration's `label` attribute. 715 + If specified as a `str`, include a `Dataset` column with that string. 716 + If specified as `False`, do not include a `Dataset` column. 717 + **T_correl**: 718 + If `True`, include correlations between all `T` values. 719 + **D47_correl**: 720 + If `True`, include correlations between all `D47` values. 721 722 ### Example 723 724 ````py 725 D47calib.huyghe_2022.export_data( 726 csvfile = 'example_export_data.csv', 727 T_correl = True, 728 D47_correl = True, 729 ) 730 ```` 731 732 This should result in something like this ([link](example_export_data.csv)): 733 734 .. include:: ../../docs/example_export_data.md 735 736 """ 737 n = len(str(self.N)) 738 739 with open(csvfile, 'w') as f: 740 f.write(sep.join(['ID', 'Sample', 'T', 'SE_T', 'D47', 'SE_D47'])) 741 742 if label: 743 f.write(f'{sep}Dataset') 744 745 if T_correl: 746 inv_diag_sT = _np.diag(_np.diag(self.sT)**-.5) 747 Tcorrel = inv_diag_sT @ self.sT @ inv_diag_sT 748 f.write(sep.join(['']+[f'Tcorrel_{k+1:0{n}d}' for k in range(self.N)])) 749 750 if D47_correl: 751 inv_diag_sD47 = _np.diag(_np.diag(self.sD47)**-.5) 752 D47correl = inv_diag_sD47 @ self.sD47 @ inv_diag_sD47 753 f.write(sep.join(['']+[f'D47correl_{k+1:0{n}d}' for k in range(self.N)])) 754 755 for k, (s, T, sT, D47, sD47) in enumerate(zip( 756 self.samples, 757 self.T, 758 _np.diag(self.sT)**.5, 759 self.D47, 760 _np.diag(self.sD47)**.5, 761 )): 762 f.write('\n' + sep.join([f'{k+1:0{n}d}', s, f'{T:.2f}', f'{sT:.2f}', f'{D47:.4f}', f'{sD47:.4f}'])) 763 if label: 764 if label is True: 765 f.write(f'{sep}{self.label}') 766 else: 767 f.write(f'{sep}{label}') 768 if T_correl: 769 f.write(sep.join(['']+[ 770 f'{Tcorrel[k,_]:.0f}' 771 if f'{Tcorrel[k,_]:.6f}'[-6:] == '000000' 772 else f'{Tcorrel[k,_]:.6f}' 773 for _ in range(self.N)])) 774 if D47_correl: 775 f.write(sep.join(['']+[ 776 f'{D47correl[k,_]:.0f}' 777 if f'{D47correl[k,_]:.6f}'[-6:] == '000000' 778 else f'{D47correl[k,_]:.6f}' 779 for _ in range(self.N)]))
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):
782 def export(self, name, filename): 783 """ 784 Save `D47calib` object as an importable file. 785 786 ### Parameters 787 788 + **name**: 789 The name of the variable to export. 790 + **filename**: 791 The filename to write to. 792 793 ### Example 794 795 ````py 796 D47calib.anderson_2021_lsce.export('foo', 'bar.py') 797 ```` 798 799 This should result in a `bar.py` file with the following contents: 800 801 ````py 802 foo = D47calib( 803 samples = ['LGB-2', 'DVH-2'], 804 T = [7.9, 33.7], 805 D47 = [0.6485720997671647, 0.5695972909966959], 806 sT = [[0.04000000000000001, 0.0], [0.0, 0.04000000000000001]], 807 sD47 = [[8.72797097773764e-06, 2.951894073404263e-06], [2.9518940734042614e-06, 7.498611746762038e-06]], 808 description = 'Devils Hole & Laghetto Basso from Anderson et al. (2021), processed in I-CDES', 809 label = 'Slow-growing calcites from Anderson et al. (2021)', 810 color = (0, 0.5, 0), 811 degrees = [0, 2], 812 bfp = {'a0': 0.1583220210575451, 'a2': 38724.41371782721}, 813 bfp_CM = [[0.00035908667755871876, -30.707016431538836], [-30.70701643153884, 2668091.396598919]], 814 chisq = 6.421311854486162e-27, 815 Nf = 0, 816 ) 817 ```` 818 """ 819 with open(filename, 'w') as f: 820 f.write(f''' 821{name} = D47calib( 822 samples = {self.samples}, 823 T = {list(self.T)}, 824 D47 = {list(self.D47)}, 825 sT = {[list(l) for l in self.sT]}, 826 sD47 = {[list(l) for l in self.sD47]}, 827 degrees = {self.degrees}, 828 description = {repr(self.description)}, 829 name = {repr(self.name)}, 830 label = {repr(self.label)}, 831 bfp = {self.bfp}, 832 bfp_CM = {[list(l) for l in self.bfp_CM]}, 833 chisq = {self.chisq}, 834 cholesky_residuals = {list(self.cholesky_residuals)}, 835 aic = {self.aic}, 836 bic = {self.bic}, 837 ks_pvalue = {self.ks_pvalue}, 838 ) 839''')
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=[]):
841def combine_D47calibs(calibs, degrees = [0,2], same_T = [], exclude_samples = []): 842 ''' 843 Combine data from several `D47calib` instances. 844 845 ### Parameters 846 847 + **calibs**: 848 A list of `D47calib` instances 849 + **degrees**: 850 The polynomial degrees of the combined regression. 851 + **same_T**: 852 Use this `list` to specify when samples from different calibrations are known/postulated 853 to have formed at the same temperature (e.g. `DVH-2` and `DHC2-8` from the `fiebig_2021` 854 and `anderson_2021_lsce` data sets). Each element of `same_T` is a `list` with the names 855 of two or more samples formed at the same temperature. 856 + **exclude_samples**: Use this `list` to specify the names of samples to exclude from 857 the combined calibration. 858 859 For example, the `OGLS23` calibration is computed with: 860 861 `same_T = [['DVH-2', DHC-2-8'], ['ETH-1-1100-SAM', 'ETH-1-1100']]` 862 863 Note that when samples from different calibrations have the same name, 864 it is not necessary to explicitly list them in `same_T`. 865 866 Also note that the regression will fail if samples listed together in `same_T` 867 actually have different `T` values specified in the original calibrations. 868 869 ### Example 870 871 The `devils_laghetto_2023` calibration is computed using the following code: 872 873 ````py 874 K = [fiebig_2021.samples.index(_) for _ in ['LGB-2', 'DVH-2', 'DHC2-8']] 875 876 fiebig_temp = D47calib( 877 samples = [fiebig_2021.samples[_] for _ in K], 878 T = fiebig_2021.T[K], 879 D47 = fiebig_2021.D47[K], 880 sT = fiebig_2021.sT[K,:][:,K], 881 sD47 = fiebig_2021.sD47[K,:][:,K], 882 ) 883 884 devils_laghetto_2023 = combine_D47calibs( 885 calibs = [anderson_2021_lsce, fiebig_temp], 886 degrees = [0,2], 887 same_T = [{'DVH-2', 'DHC2-8'}], 888 ) 889 ```` 890 ''' 891 892 samples = [s for c in calibs for s in c.samples] 893 T = [t for c in calibs for t in c.T] 894 D47 = [x for c in calibs for x in c.D47] 895 sD47 = _block_diag(*[c.sD47 for c in calibs]) 896 sT = _block_diag(*[c.sT for c in calibs]) 897 898 for i in range(len(samples)): 899 for j in range(len(samples)): 900 if i != j: 901 if (samples[i] == samples[j] or 902 any([samples[i] in _ and samples[j] in _ for _ in same_T])): 903 904 sT[i,j] = (sT[i,i] * sT[j,j])**.5 905 906 k = [_ for _, s in enumerate(samples) if s not in exclude_samples] 907 908 calib = D47calib( 909 samples = [samples[_] for _ in k], 910 T = [T[_] for _ in k], 911 D47 = [D47[_] for _ in k], 912 sT = sT[k,:][:,k], 913 sD47 = sD47[k,:][:,k], 914 degrees = degrees, 915 ) 916 917 return calib
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'}],
)