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