D95eq
Test for clumped isotope equilibrium and estimate carbonate formation temperatures from dual clumped isotope measurements
1. Installation
1.1 Recommended method
The recommended way is to use via uv (https://docs.astral.sh/uv).
If you only want to run the command-line interface (CLI): after installing uv, this should be as simple as uvx D95eq or uv tool install D95eq.
If you want to import D95eq in some Python code, once you are within a uv project (uv init), you can install the module with uv add D95eq.
After installation, open a new shell window and try D95eq --help.
1.2 Other methods
You can of course install globally via pip (pip install D95eq), or only install the CLI using pipx (pipx install D95eq).
2. Command-line interface
D95eq also provides a command-line interface (CLI).
2.1 Simple examples
(work in progress)
1""" 2Test for clumped isotope equilibrium and estimate carbonate formation temperatures from dual clumped isotope measurements 3 4.. include:: ../../docpages/install.md 5.. include:: ../../docpages/cli.md 6 7* * * 8""" 9 10from __future__ import annotations 11from ._metadata import * 12from ._tools import confidence_band 13 14import sys 15import numpy as _np 16import ogls as _ogls 17import uncertainties as _uc 18import lmfit as _lmfit 19import correldata as _cd 20import typer as _typer 21 22from typing import TYPE_CHECKING 23if TYPE_CHECKING: 24 from matplotlib import pyplot as _ppl 25 from matplotlib.patches import Ellipse as _Ellipse 26 from matplotlib.patches import Polygon as _Polygon 27 28from uncertainties import unumpy as _unp 29from scipy.stats import chi2 as _chi2 30from scipy.stats import norm as _norm 31from scipy.linalg import eigh as _eigh 32from scipy.linalg import cholesky as _cholesky 33from scipy.optimize import fsolve as _fsolve 34from numpy.typing import ArrayLike 35from typing_extensions import Annotated as _Annotated 36from typer import rich_utils as _rich_utils 37 38from warnings import filterwarnings as _filterwarnings 39_filterwarnings('ignore', category = FutureWarning, message = 'AffineScalarFunc') 40_filterwarnings('ignore', category = RuntimeWarning, message = 'The iteration is not making good progress') 41 42 43### Mathematical functions ### 44 45 46def ufloat_compatible_interp( 47 xi: (_cd.uarray | ArrayLike), 48 yi: (_cd.uarray | ArrayLike), 49 x: (float | _uc.UFloat | _cd.uarray | ArrayLike), 50): 51 """ 52 Linear interpolation accepting UFloat values for all three input parameters. 53 Only handles one interpolated value. For interpolated arrays, use `uarray_compatible_interp()` 54 55 **Arguments** 56 * `xi`: x-values defining the interpolated function 57 * `yi`: y-values defining the interpolated function 58 * `x`: x-value of the interpolation point 59 60 Returns y-value of the interpolation point, either as a float or a UFloat. 61 """ 62 xn = x.nominal_value if isinstance(x, _uc.UFloat) else float(x) 63 idx = _np.searchsorted(xi, xn) 64 idx = _np.clip(idx, 1, len(xi) - 1) 65 66 x0 = xi[idx-1] 67 x1 = xi[idx] 68 y0 = yi[idx-1] 69 y1 = yi[idx] 70 71 t = (x - x0) / (x1 - x0) 72 return y0 + t * (y1 - y0) 73 74 75def uarray_compatible_interp(xi, yi): 76 """ 77 Linear interpolation accepting UFloat values for all three input parameters. 78 79 **Arguments** 80 * `xi`: x-values defining the interpolated function 81 * `yi`: y-values defining the interpolated function 82 83 Returns an interpolation function which returns arrays or uarrays of y-values. 84 """ 85 return _np.vectorize( 86 lambda x: ufloat_compatible_interp(xi, yi, x) 87 ) 88 89 90def transform_pdf_monotonic(f_inv, df_inv, mu_x, sigma_x, yi): 91 """ 92 Compute probability distribution function of Y = f(X) 93 where X ~ Normal(mu_x, sigma_x) and f is monotonic, 94 based on the change-of-variables formula: 95 96 p[y=f(x)] = p[x=f_inv(y)] * d(f_inv)/dy 97 98 Additionally, if f_inv returns UFloats, the PDF is convolved with that local 99 source of uncertainty (assumed to be Gaussian) at each grid point. 100 101 As currently implemented, requires `yi` to be an equally spaced array-like. 102 103 **Arguments** 104 f_inv: inverse of f, may return UFloats 105 df_inv: derivative of f_inv, should return UFloats if f_inv does 106 mu_x: mean of X PDF 107 sigma_x: std dev of X PDF 108 yi: regularly spaced grid of y values at which to evaluate the PDF 109 110 **Returns:** 111 pdf: normalized PDF evaluated at yi 112 """ 113 114 if not _np.allclose(_np.diff(yi), yi[1] - yi[0]): 115 raise ValueError("yi must be regularly spaced") 116 117 xi = f_inv(yi) # may be floats or ufloats, depending on f_inv 118 119 try: 120 xi_nom = xi.n 121 sigma_xi = xi.s 122 has_ufloats = True 123 except AttributeError: 124 xi_nom = xi 125 has_ufloats = False 126 127 # Jacobian weights (account for irregular xi spacing) 128 try: 129 df_inv_nom = df_inv(yi).n 130 except AttributeError: 131 df_inv_nom = df_inv(yi) 132 133 w_i = _norm.pdf(xi_nom, loc = mu_x, scale = sigma_x) * _np.abs(df_inv_nom) 134 135 if not has_ufloats: 136 return w_i / (_np.trapezoid(w_i, yi)) 137 138 # Propagate sigma from x-space to y-space via Jacobian: sigma_y = sigma_x / abs( dx/dy ) 139 sigma_yi = sigma_xi / _np.abs(df_inv_nom) 140 141 # Convolution of Gaussians: each grid point j contributes N(yi; yj, σ_yj²) scaled by w_j 142 gaussians = _norm.pdf( 143 yi[:, None], 144 loc = yi[None, :], 145 scale = sigma_yi[None, :] 146 ) # NOTE: nice syntax to reshape ndarrays, perhaps use this in D4x_calib_function? 147 148 pdf = (gaussians * w_i[None, :]).sum(axis = 1) 149 150 return pdf / (_np.trapezoid(pdf, yi)) 151 152 153#### Calibration variables and functions #### 154 155 156_D47_approx_calib_coefs = [0.159502986, 38588.1545] # computed from code in comments below 157# from D47calib import OGLS23 as _OGLS23 158# from D47calib import D47calib as _D47calib 159# 160# _D47_approx = _D47calib( 161# samples = _OGLS23.samples, 162# T = _OGLS23.T, 163# sT = _OGLS23.sT, 164# D47 = _OGLS23.D47, 165# sD47 = _OGLS23.sD47, 166# degrees = [0,2], 167# ) 168# _D47_approx_calib_coefs = [_D47_approx.bfp['a0'], _D47_approx.bfp['a2']] 169 170 171def _compute_D48_calib_coefficients(reprocess = False): 172 """ 173 Based on Fiebig et al. (2021, 2024) 174 """ 175 176 # D64 predictions 177 a1 = 6.002 178 a2 = -1.299e4 179 a3 = 8.996e6 180 a4 = -7.423e8 181 182 if reprocess: 183 184 # M. Bernecker, pers. comm. 185 # after Fiebig et al. (2024) 10.1016/j.chemgeo.2024.122382 186 datastr = ''' 187 Sample, D48, SE_D48, T, SE_T, correl_T 188 LGB-2, 0.2606, 0.0103, 7.9, 0.2, 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 189 DHC2-8, 0.2335, 0.0066, 33.7, 0.2, 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0. 190 DVH-2, 0.2484, 0.0105, 33.7, 0.2, 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0. 191 CA120, 0.1715, 0.0154, 120.0, 2., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0. 192 CA170, 0.1621, 0.0142, 170.0, 2., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0. 193 CA200, 0.1561, 0.0134, 200.0, 2., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0. 194 CA250A, 0.1449, 0.0146, 250.0, 2., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0. 195 CA250B, 0.1301, 0.0134, 250.0, 2., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0. 196 CM351, 0.1220, 0.0073, 726.85, 10., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0. 197 ETH-1-1100, 0.1161, 0.0091, 1100.0, 10., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0. 198 ETH-2-1100, 0.1225, 0.0070, 1100.0, 10., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1. 199 '''[1:-2] 200 201 data = _cd.read_str(datastr) 202 T, D48 = data['T'], data['D48'] 203 204 205 D64_predicted = ( 206 a1 / (273.15 + T) 207 + a2 / (273.15 + T)**2 208 + a3 / (273.15 + T)**3 209 + a4 / (273.15 + T)**4 210 ) 211 212 # affine regression of the form D48 = b0 + b1 * D64_theory 213 R = _ogls.Polynomial( 214 X = D64_predicted.n, 215 sX = D64_predicted.covar, 216 Y = D48.n, 217 sY = D48.covar, 218 degrees = [0,1], 219 ) 220 221 R.regress(overdispersion_scaling = True) 222 b0, b1 = _uc.correlated_values(R.bfp.values(), R.bfp_CM) 223# print(_cd.data_string(dict(affine_coefs = _cd.uarray([b0, b1])))) 224 225 else: 226 227 # M. Bernecker, pers. comm. 228 # after Fiebig et al. (2024) 10.1016/j.chemgeo.2024.122382 229 # Caution: because Fiebig et al. ignored T uncertainties, these 230 # coefficeients have smaller uncertainties than those computed above. 231 b0, b1 = _uc.correlated_values( 232 [ 233 0.12135157920099604, 234 1.0379702801201238, 235 ], [ 236 [ 7.39697438e-06, -6.90467053e-05], 237 [-6.90467053e-05, 1.46002771e-03], 238 ], 239 ) 240 241 a0 = b0 242 a1 *= b1 243 a2 *= b1 244 a3 *= b1 245 a4 *= b1 246 247 return _cd.uarray([a0, a1, a2, a3, a4]) 248 249 250def D4x_calib_function( 251 T: (float | _uc.UFloat | _cd.uarray | ArrayLike), 252 coefs: _cd.uarray, 253 return_without_uncertainties: bool = False, 254 ignore_calib_uncertainties: bool = False, 255) -> (float | _uc.UFloat | _cd.uarray | ArrayLike): 256 """ 257 **Arguments** 258 * `T`: temperature(s) for which to compute Δ<sub>4x</sub> 259 * `return_without_uncertainties`: if `True`, returns Δ<sub>4x</sub> values without error propagation of any kind 260 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 261 262 Returns equilibrium Δ<sub>4x</sub> value(s) corresponding to `T` value(s) 263 """ 264 degs = _np.arange(coefs.size) 265 266 D4x = ( 267 _np.expand_dims(_cd.nv(coefs) if ignore_calib_uncertainties else coefs, 1) # shape = (coefs.size, 1) 268 * _np.expand_dims((T+273.15)**-1, 0) # shape = (1, T.size) 269 ** _np.expand_dims(degs, 1) # shape = (coefs.size, 1) 270 ).sum(axis = 0 if isinstance(T, _np.ndarray) else None) 271 272 if D4x.ndim == 0: 273 return D4x.tolist().n if return_without_uncertainties else D4x.tolist() 274 return D4x.n if return_without_uncertainties else D4x 275 276 277def D4x_calib_derivative( 278 T: (float | _uc.UFloat | _cd.uarray | ArrayLike), 279 coefs: _cd.uarray, 280 return_without_uncertainties: bool = False, 281 ignore_calib_uncertainties: bool = False, 282) -> (float | _uc.UFloat | _cd.uarray | ArrayLike): 283 """ 284 **Arguments** 285 * `T`: temperature(s) for which to compute Δ<sub>4x</sub> 286 * `return_without_uncertainties`: if `True`, returns D4x values without error propagation of any kind. 287 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties. 288 289 Returns d(D4x)/dT corresponding to `T` value(s) 290 """ 291 dcoefs = -_np.arange(len(coefs)) * coefs 292 dcoefs = _cd.uarray((dcoefs[0], *dcoefs)) 293 return D4x_calib_function( 294 T, 295 dcoefs, 296 return_without_uncertainties = return_without_uncertainties, 297 ignore_calib_uncertainties = ignore_calib_uncertainties, 298 ) 299 300 301#### Plotting functions #### 302 303 304def conf_ellipse( 305 X: (_cd.uarray | _np.ndarray | _uc.UFloat | float), 306 Y: (_cd.uarray | _np.ndarray | _uc.UFloat | float) = None, 307 p: float = 0.95, 308 CM: (_np.ndarray | None) = None, 309 Xse: (_np.ndarray | float | None) = None, 310 Yse: (_np.ndarray | float | None) = None, 311 plot: bool = True, 312 ax: (_ppl.Axes | None) = None, 313 **kwargs, 314) -> tuple: 315 """ 316 Compute and (optionally) plot the joint *p*-level confidence ellipses for the elements of (X, Y) 317 318 **Arguments** 319 * `X`: x values 320 * `Y`: y values 321 * `p`: confidence level 322 * `CM`: covariance matrix of (X, Y); not needed if X and Y are of type 323 [`uncertainties.UFloat`](https://pythonhosted.org/uncertainties/tech_guide.html). 324 or if (`Xse`, `Yse`) are specified. 325 * `Xse`, `Yse`: SE of X and Y; not needed if X and Y are of type 326 [`uncertainties.UFloat`](https://pythonhosted.org/uncertainties/tech_guide.html) 327 or if `CM` is specified. 328 * `plot`: whether to plot the ellipse or not. If `False`, return a list of 329 `(x_center, y_center, width, height, angle)` elements 330 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 331 * `kwargs`: passed to `matplotlib.patches.Ellipse()` 332 333 Returns a list of the `Ellipse` objects thus created. 334 """ 335 336 r2 = _chi2.ppf(p, 2) 337 kwargs = dict(fc = 'None', ec = 'k', lw = 0.7) | kwargs 338 339 out = [] 340 341 for x, y in zip( 342 *_cd.as_pair_of_uarrays(X, Y, CM = CM, Xse = Xse, Yse = Yse) 343 ): 344 val, vec = _eigh(_uc.covariance_matrix((x, y))) 345 width, height = 2 * (val[:, None] * r2)**0.5 346 angle = _np.degrees(_np.arctan2(*vec[::-1, 0])) 347 348 if plot: 349 from matplotlib import pyplot as _ppl 350 from matplotlib.patches import Ellipse as _Ellipse 351 352 if ax is None: 353 ax = _ppl.gca() 354 355 out.append( 356 ax.add_patch( 357 _Ellipse( 358 xy = (x.n, y.n), 359 width = width.item(), 360 height = height.item(), 361 angle = angle, 362 **kwargs, 363 ) 364 ) 365 ) 366 else: 367 out.append([x.n, y.n, width, height, angle]) 368 369 return (*out,) 370 371 372### D95eq Engine implementation ### 373 374class _Interpolation(): 375 pass 376 377class Engine(): 378 """ 379 Underlying engine to compute and plot nearest equilibrium temperatures and projected 380 temperatures based on a consistent pair of Δ<sub>47</sub>, Δ<sub>48</sub> calibrations. 381 """ 382 383 # D47_calib_coefs from OGLS23 (D47calib v1.3.1) 384 D47_calib_coefs = _cd.read_str(''' 385 coefs, SE, correl, 3860.17437754366432887, 4.911105567257293e-3, 1. , -0.93797005, 0.8865771 387 -18.14215245127414, 5.632326472234856, -0.93797005, 1. , -0.98994249 38842.65722989162373e3, 1.27712751715908e3, 0.8865771 , -0.98994249, 1. 389'''[1:-1])['coefs'] 390 """ 391 Default (OGLS23) Δ<sub>47</sub> calibration coefficients based on [Daëron & Vermeesch (2024)](https://doi.org/10.1016/j.chemgeo.2023.121881) 392 """ 393 394 # D48_calib_coefs reprocessed from Fiebig et al. (2024): 395 # 396 # D48_calib_coefs = _compute_D48_calib_coefficients(reprocess = True) 397 # print(_cd.data_string( 398 # {'coefs': D48_calib_coefs}, 399 # float_format = 'z.12g', 400 # correl_format = 'z.12f', 401 # )) 402 403 D48_calib_coefs = _cd.read_str(''' 404 coefs, SE_coefs, correl_coefs, , , , 4050.121349237888, 0.00390048540724, 1.000000000000, -0.664181963395, 0.664181963395, -0.664181963395, 0.664181963395 406 6.22931985613, 0.32896761459, -0.664181963395, 1.000000000000, -1.000000000000, 1.000000000000, -1.000000000000 407 -13481.983494, 711.977559735, 0.664181963395, -1.000000000000, 1.000000000000, -1.000000000000, 1.000000000000 408 9336714.66607, 493067.754224, -0.664181963395, 1.000000000000, -1.000000000000, 1.000000000000, -1.000000000000 409-770413883.573, 40685214.9801, 0.664181963395, -1.000000000000, 1.000000000000, -1.000000000000, 1.000000000000 410'''[1:-1])['coefs'] 411 """ 412 Default Δ<sub>48</sub> calibration coefficients based on [Fiebig et al. (2024)](https://doi.org/10.1016/j.chemgeo.2024.122382) 413 """ 414 415 def __init__( 416 self, 417 D47_coefs: (_cd.uarray | ArrayLike | None) = None, 418 D48_coefs: (_cd.uarray | ArrayLike | None) = None, 419 Tmin_interp: float = -23.0, 420 Tmax_interp: float = 1277.0, 421 N_interp: float = 201, 422 ): 423 """ 424 **Arguments** 425 * `D47_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 426 * `D48_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 427 * `Tmin_interp`: minimum temperature over which to interpolate for inverse function computations 428 * `Tmax_interp`: maximum temperature over which to interpolate for inverse function computations 429 * `N_interp`: number of points (equally-spaced in 1/T space) over which to interpolate for inverse function computations 430 """ 431 432 self.D47_coefs = Engine.D47_calib_coefs if D47_coefs is None else D47_coefs 433 """The Δ<sub>47</sub> calibration coefficients used by this `Engine` instance""" 434 435 self.D48_coefs = Engine.D48_calib_coefs if D48_coefs is None else D48_coefs 436 """The Δ<sub>48</sub> calibration coefficients used by this `Engine` instance""" 437 438 self.interp = _Interpolation() 439 """ 440 Holds equilibrium Δ<sub>47</sub> and Δ<sub>48</sub> values (ufloats) interpolated 441 along an array of T values (regularly spaced increments of 1/T<sup>2</sup>). 442 443 * `interp.T`: interpolation T values (floats) in regularly spaced increments of 1/T<sup>2</sup> 444 * `interp.D47`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T` 445 * `interp.D48`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T` 446 * `interp.D47_no_calib_errors`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T`, 447 ignoring calibration uncertainties 448 * `interp.D48_no_calib_errors`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T`, 449 ignoring calibration uncertainties 450 """ 451 452 self.interp.T = _np.linspace( 453 (Tmax_interp+273.15)**-2, 454 (Tmin_interp+273.15)**-2, 455 N_interp, 456 )**-0.5 - 273.15 457 458 self.interp.D47 = self.D47_calib_function( 459 self.interp.T, 460 return_without_uncertainties = False, 461 ignore_calib_uncertainties = False, 462 ) 463 464 self.interp.D47_no_calib_errors = self.D47_calib_function( 465 self.interp.T, 466 return_without_uncertainties = False, 467 ignore_calib_uncertainties = True, 468 ) 469 470 self.interp.D48 = self.D48_calib_function( 471 self.interp.T, 472 return_without_uncertainties = False, 473 ignore_calib_uncertainties = False, 474 ) 475 476 self.interp.D48_no_calib_errors = self.D48_calib_function( 477 self.interp.T, 478 return_without_uncertainties = False, 479 ignore_calib_uncertainties = True, 480 ) 481 482 self.interp.D47u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D47) 483 self.interp.D48u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D48) 484 485 #inverse D47 calibration (ignoring calibration errors) 486 self.interp.Teq_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.T) 487 #inverse D47 calibration (including calibration errors) 488 self.interp.Teq_as_function_of_D47u = uarray_compatible_interp(self.interp.D47, self.interp.T) 489 490 def T_as_function_of_D47( 491 self, 492 D47: (_cd.uarray | ArrayLike), 493 ignore_calib_uncertainties: bool = False, 494 ): 495 """ 496 Provided with one or more Δ<sub>47</sub> values (floats or ufloats), return ufloats for the 497 corresponding equilibrium T values (ufloats with or without Δ<sub>47</sub> calibration uncertainties). 498 499 **Arguments** 500 * `D47`: array of Δ<sub>47</sub> values 501 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 502 """ 503 if ignore_calib_uncertainties: 504 return _cd.uarray(self.interp.Teq_as_function_of_D47n(D47)) 505 else: 506 return _cd.uarray(self.interp.Teq_as_function_of_D47u(D47)) 507 508 def D47u_as_function_of_D47n( 509 self, 510 D47: ArrayLike 511 ): 512 """ 513 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 514 equilibrium Δ<sub>47</sub> values (ufloats with Δ<sub>47</sub> calibration uncertainties). 515 """ 516 return _cd.uarray(self.interp.D47u_as_function_of_D47n(D47)) 517 518 def D48u_as_function_of_D47n( 519 self, 520 D47: ArrayLike 521 ): 522 """ 523 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 524 equilibrium Δ<sub>48</sub> values (ufloats with Δ<sub>48</sub> calibration uncertainties). 525 """ 526 return _cd.uarray(self.interp.D48u_as_function_of_D47n(D47)) 527 528 def D47_calib_function( 529 self, 530 T: (float | _uc.UFloat | _cd.uarray), 531 return_without_uncertainties: bool = False, 532 ignore_calib_uncertainties: bool = False, 533 ): 534 return D4x_calib_function( 535 T = T, 536 coefs = self.D47_coefs, 537 return_without_uncertainties = return_without_uncertainties, 538 ignore_calib_uncertainties = ignore_calib_uncertainties, 539 ) 540 541 def D48_calib_function( 542 self, 543 T: (float | _uc.UFloat | _cd.uarray), 544 return_without_uncertainties: bool = False, 545 ignore_calib_uncertainties: bool = False, 546 ): 547 return D4x_calib_function( 548 T = T, 549 coefs = self.D48_coefs, 550 return_without_uncertainties = return_without_uncertainties, 551 ignore_calib_uncertainties = ignore_calib_uncertainties, 552 ) 553 554 D47_calib_function.__doc__ = D4x_calib_function.__doc__.replace('Δ<sub>4x</sub>', 'Δ<sub>47</sub>') 555 D48_calib_function.__doc__ = D4x_calib_function.__doc__.replace('Δ<sub>4x</sub>', 'Δ<sub>48</sub>') 556 557 def T_ellipse( 558 self, 559 T: (_np.ndarray | _cd.uarray), 560 p: float = 0.95, 561 CM: (_np.ndarray | None) = None, 562 Tse: (_np.ndarray | float | None) = None, 563 plot: bool = True, 564 ax: (_ppl.Axes | None) = None, 565 **kwargs, 566 ) -> list: 567 """ 568 Plot the joint `p`-level confidence ellipses in (Δ<sub>47</sub>, Δ<sub>48</sub>) 569 space, for temperatures equal to the elements of `T`, and return a list of the 570 `Ellipse` objects thus created. 571 572 **Arguments** 573 * `T`: `ndarray` or `uarray` of temperatures to plot 574 * `p`: confidence level 575 * `plot`: whether to plot the ellipse or not. If `False`, return a list of 576 `(x_center, y_center, width, height, angle)` elements 577 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 578 * `kwargs`: passed to `matplotlib.patches.Ellipse()` 579 """ 580 _T = _cd.as_uarray(T, CM = CM, Xse = Tse) 581 return conf_ellipse( 582 self.D47_calib_function(_T), 583 self.D48_calib_function(_T), 584 p = p, 585 plot = plot, 586 ax = ax, 587 **kwargs, 588 ) 589 590 def plot_D95_confidence_band( 591 self, 592 p: float = 0.95, 593 Ti: (ArrayLike | None) = None, 594 plot: bool = True, 595 ax: (_ppl.Axes | None) = None, 596 **kwargs, 597 ): 598 """ 599 Plot, for a given p-value, the confidence band of the thermodynamic equilibrium curve 600 in (Δ<sub>47</sub>, Δ<sub>48</sub>) space. 601 602 **Arguments** 603 * `p`: confidence level 604 * `Ti`: array of temperatures over which to evaluate confidence band (default: use `interp.T` attribute instead) 605 * `plot`: whether to plot the confidence band or not. If `False`, return the (N,2) array of polygon nodes 606 * `ax`: `Axes` instance to plot to (default: use current Axes) 607 * `kwargs`: passed to `patches.Polygon()` 608 609 Returns the corresponding `Polygon` instance. 610 """ 611 612 if Ti is None: 613 Ti = self.interp.T 614 615 cb = confidence_band( 616 Ti, 617 self.D47_calib_function, 618 self.D48_calib_function, 619 p, 620 ) 621 622 if plot: 623 from matplotlib import pyplot as _ppl 624 from matplotlib.patches import Polygon as _Polygon 625 626 if ax is None: 627 ax = _ppl.gca() 628 629 polygon = ax.add_patch( 630 _Polygon( 631 cb, 632 closed = True, 633 **kwargs, 634 ) 635 ) 636 return polygon 637 else: 638 return cb 639 640 641 def plot_D95_equilibrium( 642 self, 643 Tmin: float = 0., 644 Tmax: float = 1000., 645 NT: int = 101, 646 Tmarkers: _np.typing.ArrayLike = [0, 25, 100, 250, 1000], 647 kwargs_Tmarkers: dict = {}, 648 show_Tmarker_labels: bool = True, 649 kwargs_Tmarker_labels: dict = {}, 650 show_Tmarker_ellipses: bool = False, 651 kwargs_Tmarker_ellipses: dict = {}, 652 show_eqline: bool = True, 653 kwargs_eqline: dict = {}, 654 show_confidence: bool = True, 655 confidence_pvalue: float = 0.95, 656 kwargs_confidence: dict = {}, 657 ax: (_ppl.Axes | None) = None, 658 xlabel: str = '$Δ_{47}$ [‰]', 659 ylabel: str = '$Δ_{48}$ [‰]', 660 lw: float = 0.7, 661 ) -> (dict, dict): 662 """ 663 Plot a thermodynamic equilibrium curve in (Δ<sub>47</sub>, Δ<sub>48</sub>) space 664 as a function of temperature. 665 666 **Arguments** 667 * `Tmin`: minimum T to plot 668 * `Tmax`: maximum T to plot 669 * `NT`: number of steps in equilibrium curve (interpolated at constant steps in 1/T<sup>2</sup> space) 670 * `Tmarkers`: T markers to add along the curve 671 * `kwargs_Tmarkers`: passed to `plot()` when plotting T markers 672 * `show_Tmarker_labels`: whether to add T labels to T markers 673 * `kwargs_Tmarker_labels`: passed to `text()` when plotting T markers 674 * `show_Tmarker_ellipses`: whether to add confidence ellipses to T markers 675 * `kwargs_Tmarker_ellipses`: passed to `T_ellipses()` when plotting T marker ellipses 676 * `show_eqline`: whether to plot the equilibrium curve itself 677 * `kwargs_eqline`: passed to `plot()` when plotting the equilibrium curve 678 * `show_confidence`: whether to plot the confidence band of the equilibrium curve 679 * `confidence_pvalue`: confidence level for the confidence band 680 * `kwargs_confidence`: passed to `plot_D95_confidence_band()` when plotting the confidence band 681 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 682 * `xlabel`: string to pass to `xlabel()` 683 * `ylabel`: string to pass to `ylabel()` 684 * `lw`: default line width for most plot elements 685 686 **Returns** 687 * `data`: a dict of the T, Δ<sub>47</sub> and Δ<sub>48</sub> values generated for this plot: 688 - `Te` : temperature interpolated along the equilibrium curve 689 - `D47e`: Δ<sub>47</sub> interpolated along the equilibrium curve 690 - `D48e`: Δ<sub>48</sub> interpolated along the equilibrium curve 691 - `Tm` : temperature of T markers 692 - `D47m`: Δ<sub>47</sub> of T markers 693 - `D48m`: Δ<sub>48</sub> of T markers 694 695 * `plot_elements`: a dict of the `Axes` elements generated for this plot: 696 - `eqline`: `Line2D` of the equilibrium curve 697 - `confidence`: `Polygon` object for the confidence band 698 - `Tm`: `Line2D` of the T markers 699 - `Tme`: list of `Ellipse` objects for the T marker ellipses 700 - `Tml`: list of `Text` objects for the T marker labels 701 """ 702 703 from matplotlib import pyplot as _ppl 704 705 default_kwargs_eqline = dict( 706 marker = 'None', 707 ls = '-', 708 color = 'k', 709 lw = lw, 710 ) 711 default_kwargs_confidence = dict( 712 ec = (0,0,0,1), 713 fc = (0,0,0,0.15), 714 lw = 0., 715 ) 716 default_kwargs_Tmarkers = dict( 717 ls = 'None', 718 marker = 'o', 719 ms = 4, 720 mfc = 'w', 721 mec = 'k', 722 mew = lw, 723 ) 724 default_kwargs_Tmarker_ellipses = dict( 725 fc = 'None', 726 ec = 'k', 727 lw = lw, 728 ) 729 default_kwargs_Tmarker_labels = dict( 730 size = 8, 731 va = 'center', 732 ha = 'left', 733 linespacing = 3, 734 ) 735 736 plot_elements = {} 737 738 Ti = _np.linspace( 739 (Tmin + 273.15)**-2, 740 (Tmax + 273.15)**-2, 741 NT 742 )**-0.5 - 273.15 743 744 Tmarkers = _np.array([_ for _ in Tmarkers if _ >= Ti.min() and _ <= Ti.max()]) 745 746 if ax is None: 747 ax = _ppl.gca() 748 ax.set_xlabel(xlabel) 749 ax.set_ylabel(ylabel) 750 751 Xe = self.D47_calib_function(Ti) 752 Ye = self.D48_calib_function(Ti) 753 754 if show_eqline: 755 plot_elements['eqline'], = ax.plot( 756 _unp.nominal_values(Xe), 757 _unp.nominal_values(Ye), 758 **(default_kwargs_eqline | kwargs_eqline), 759 ) 760 761 if show_confidence: 762 plot_elements['confidence'] = self.plot_D95_confidence_band( 763 p = confidence_pvalue, 764 ax = ax, 765 **(default_kwargs_confidence | kwargs_confidence), 766 ) 767 768 Xm = self.D47_calib_function(Tmarkers) 769 Ym = self.D48_calib_function(Tmarkers) 770 if Tmarkers.size > 0: 771 plot_elements['Tm'] = ax.plot( 772 _unp.nominal_values(Xm), 773 _unp.nominal_values(Ym), 774 **(default_kwargs_Tmarkers | kwargs_Tmarkers), 775 ) 776 if show_Tmarker_ellipses: 777 plot_elements['Tme'] = conf_ellipse( 778 Xm, 779 Ym, 780 ax = ax, 781 **(default_kwargs_Tmarker_ellipses | kwargs_Tmarker_ellipses), 782 ) 783 if show_Tmarker_labels: 784 plot_elements['Tml'] = [] 785 for x,y,t in zip(Xm, Ym, Tmarkers): 786 plot_elements['Tml'].append( 787 ax.text( 788 x.n, 789 y.n, 790 f'\n${t:.0f}\\,$°C', 791 **(default_kwargs_Tmarker_labels | kwargs_Tmarker_labels), 792 ) 793 ) 794 795 ax.autoscale_view() 796 797 data = dict( 798 Te = Ti, 799 D47e = Xe, 800 D48e = Ye, 801 Tm = Tmarkers, 802 D47m = Xm, 803 D48m = Ym, 804 ) 805 806 return data, plot_elements 807 808 def _compute_p_and_D48eq_from_D47eq( 809 self, 810 D47, 811 D48, 812 D47eq, 813 ignore_calib_uncertainties = False, 814 ): 815 """ 816 Used by the various `Engine.nearest_D47eq()` methods 817 """ 818 N = D47.size 819 820 # Compute fit residuals for p values 821 if ignore_calib_uncertainties: 822 R = _cd.uarray(_np.concatenate(( 823 D47 - self.D47u_as_function_of_D47n(D47eq.n).n, 824 D48 - self.D48u_as_function_of_D47n(D47eq.n).n, 825 ))) 826 else: 827 R = _cd.uarray(_np.concatenate(( 828 D47 - self.D47u_as_function_of_D47n(D47eq.n), 829 D48 - self.D48u_as_function_of_D47n(D47eq.n), 830 ))) 831 832 # Compute p values 833 p = _np.zeros((N,)) 834 for k in range(N): 835 r = R[k::N] 836 z2 = r.m 837 p[k] = 1-_chi2.cdf(z2, 1) 838 839 # Compute D48eq 840 D48eq = self.D48u_as_function_of_D47n(D47eq) 841 842 return p, D48eq 843 844 def nearest_D47eq( 845 self, 846 D47: _cd.uarray, 847 D48: _cd.uarray, 848 ignore_calib_uncertainties: bool = False, 849 ): 850 """ 851 Computes a `correldata.uarray` of *equilibrium* Δ<sub>47</sub> values, each of which is 852 the closest (in the OGLS sense) to one (Δ<sub>47</sub>, Δ<sub>48</sub>) observation 853 considered independently of the others. 854 855 Also returns an array of corresponding p-values taking into account errors in Δ<sub>47</sub> 856 and Δ<sub>48</sub> (and any covariance between the two) as well as errors in the 857 Δ<sub>47</sub> and Δ<sub>48</sub> calibrations. 858 859 > [!NOTE] 860 > This is both the fastest and the strongly recommended version of this calculation. 861 > It is expected to yield an `uarray` with reasonably accurate covariance between the 862 > `D47eq` values, but also between `D47eq` and all other variables. 863 """ 864 865 N = D47.size 866 N47 = self.D47_coefs.size 867 N48 = self.D48_coefs.size 868 D47eq = D47 * 0 869 870 # _np.set_printoptions(threshold = _np.inf) 871 # _np.set_printoptions(linewidth = _np.inf) 872 873 for i in range(N): 874 def fun(*args): # args = (D47, D48, *D47_calib_coefs, *D48_calib_coefs) 875 876 args = _np.array(args) 877 D47_n = args[0] 878 D48_n = args[1] 879 D47_calib_coefs_n = args[-N48-N47:-N48] 880 D48_calib_coefs_n = args[-N48:] 881 882 params = _lmfit.Parameters() 883 params.add('D47eq', value = D47_n) 884 885 D47_u = _cd.uarray([_uc.ufloat(D47_n, D47.s[i])]) 886 D48_u = _cd.uarray([_uc.ufloat(D48_n, D48.s[i])]) 887 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 888 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 889 890 D47i = D4x_calib_function( 891 self.interp.T, 892 D47_calib_coefs_u, 893 return_without_uncertainties = False, 894 ignore_calib_uncertainties = ignore_calib_uncertainties, 895 ) 896 897 D48i = D4x_calib_function( 898 self.interp.T, 899 D48_calib_coefs_u, 900 return_without_uncertainties = False, 901 ignore_calib_uncertainties = ignore_calib_uncertainties, 902 ) 903 904 D47_interp = uarray_compatible_interp(D47i.n, D47i) 905 D48_interp = uarray_compatible_interp(D47i.n, D48i) 906 907 def cost_fun(p): 908 R = _cd.uarray(_np.concatenate(( 909 D47_u - D47_interp(p['D47eq'].value), 910 D48_u - D48_interp(p['D47eq'].value), 911 ))) 912 913 invS = _np.linalg.inv(R.covar) 914 L = _cholesky(invS) 915 916 return L @ R.n 917 918 minresult = _lmfit.minimize( 919 cost_fun, 920 params, 921 method = 'least_squares', 922 scale_covar = False, 923 jac = '3-point', 924 ) 925 # slower but yields very similar results: 926 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 927 928 return minresult.params['D47eq'].value 929 930 wrapped_fun = _uc.wrap(fun) 931 D47eq[i] = wrapped_fun(D47[i], D48[i], *self.D47_coefs, *self.D48_coefs) 932 933 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 934 935 return D47eq, D48eq, p 936 937 def joint_nearest_D47eq( 938 self, 939 D47: _cd.uarray, 940 D48: _cd.uarray, 941 ignore_calib_uncertainties: bool = False, 942 ): 943 """ 944 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 945 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 946 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 947 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 948 Δ<sub>48</sub> calibrations. 949 950 > [!CAUTION] 951 > Caution: the use of this function is **not generally recommended** except for 952 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 953 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 954 > as done by the recommended function `nearest_D47eq()`. 955 956 This is the most complete but slowest and not recommended version of this calculation. 957 It is expected to yield an `uarray` with reasonably accurate covariance between the 958 `D47eq` values, but also between `D47eq` and all other variables. 959 960 A faster but incomplete and potentially less accurate version of this calculation is 961 provided by `lazy_joint_nearest_D47eq()`. 962 """ 963 964 N = D47.size 965 N47 = self.D47_coefs.size 966 N48 = self.D48_coefs.size 967 968 def fun(j, *args): 969 970 args = _np.array(args) 971 D47_n = args[:N] 972 D48_n = args[N:2*N] 973 D47_calib_coefs_n = args[-N48-N47:-N48] 974 D48_calib_coefs_n = args[-N48:] 975 976 params = _lmfit.Parameters() 977 for k in range(N): 978 params.add(f'D47eq{k}', value = D47_n[k]) 979 980 D47_u = _cd.uarray(_uc.correlated_values(D47_n, D47.covar)) 981 D48_u = _cd.uarray(_uc.correlated_values(D48_n, D48.covar)) 982 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 983 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 984 985 D47i = D4x_calib_function( 986 self.interp.T, 987 D47_calib_coefs_u, 988 return_without_uncertainties = False, 989 ignore_calib_uncertainties = ignore_calib_uncertainties, 990 ) 991 992 D48i = D4x_calib_function( 993 self.interp.T, 994 D48_calib_coefs_u, 995 return_without_uncertainties = False, 996 ignore_calib_uncertainties = ignore_calib_uncertainties, 997 ) 998 999 D47_interp = uarray_compatible_interp(D47i.n, D47i) 1000 D48_interp = uarray_compatible_interp(D47i.n, D48i) 1001 1002 def cost_fun(p): 1003 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1004 R = _cd.uarray(_np.concatenate(( 1005 D47_u - D47_interp(_D47eq), 1006 D48_u - D48_interp(_D47eq), 1007 ))) 1008 1009 invS = _np.linalg.inv(R.covar) 1010 L = _cholesky(invS) 1011 1012 # print(((L @ R.n)**2).sum()) 1013 return L @ R.n 1014 1015 minresult = _lmfit.minimize( 1016 cost_fun, 1017 params, 1018 method = 'least_squares', 1019 scale_covar = False, 1020 jac = '3-point', 1021 ) 1022 # slower but yields very similar results: 1023 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 1024 1025 return minresult.params[f'D47eq{j}'].value 1026 1027 wrapped_fun = _uc.wrap(fun) 1028 1029 D47eq = _cd.uarray([wrapped_fun(j, *D47, *D48, *self.D47_coefs, *self.D48_coefs) for j in range(N)]) 1030 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1031 1032 return D47eq, D48eq, p 1033 1034 def lazy_joint_nearest_D47eq( 1035 self, 1036 D47: _cd.uarray, 1037 D48: _cd.uarray, 1038 ignore_calib_uncertainties: bool = False, 1039 ): 1040 """ 1041 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 1042 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 1043 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 1044 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 1045 Δ<sub>48</sub> calibrations. 1046 1047 > [!CAUTION] 1048 > Caution: the use of this function is **not generally recommended** except for 1049 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 1050 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 1051 > as done by the recommended function `nearest_D47eq()`. 1052 1053 This is a faster but incomplete version of this calculation. It is expected to yield an 1054 `uarray` with roughly accurate covariance between the `Teq` values, but without computing 1055 the covariance with any other variables. 1056 1057 A slower but complete and more accurate version of this calculation is provided by 1058 `joint_nearest_D47eq()`. 1059 """ 1060 1061 N = D47.size 1062 1063 params = _lmfit.Parameters() 1064 for k in range(N): 1065 params.add(f'D47eq{k}', value = D47[k].n) 1066 1067 def cost_fun(p, ignore_calib_uncertainties = ignore_calib_uncertainties): 1068 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1069 1070 if ignore_calib_uncertainties: 1071 R = _cd.uarray(_np.concatenate(( 1072 D47 - self.D47u_as_function_of_D47n(_D47eq).n, 1073 D48 - self.D48u_as_function_of_D47n(_D47eq).n, 1074 ))) 1075 else: 1076 R = _cd.uarray(_np.concatenate(( 1077 D47 - self.D47u_as_function_of_D47n(_D47eq), 1078 D48 - self.D48u_as_function_of_D47n(_D47eq), 1079 ))) 1080 1081 invS = _np.linalg.inv(R.covar) 1082 L = _cholesky(invS) 1083 1084 # print(((L @ R.n)**2).sum()) 1085 return L @ R.n 1086 1087 minresult = _lmfit.minimize( 1088 cost_fun, 1089 params, 1090 method = 'least_squares', 1091 scale_covar = False, 1092 jac = '3-point', 1093 ) 1094 1095 D47eq = _cd.uarray([minresult.uvars[f'D47eq{k}'] for k in range(N)]) 1096 1097 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1098 1099 return D47eq, D48eq, p 1100 1101 def projected_D47eq( 1102 self, 1103 D47: _cd.uarray, 1104 D48: _cd.uarray, 1105 kinetic_slope: (float | _uc.UFloat), 1106 ): 1107 """ 1108 Projects one or more (Δ<sub>47</sub>, Δ<sub>48</sub>) observations onto the equlibrium curve 1109 following a kinetic fractionation vector with a given slope (∂Δ<sub>48</sub>/∂Δ<sub>47</sub>). 1110 1111 **Arguments** 1112 * `D47`: observed Δ<sub>47</sub> value(s) 1113 * `D48`: observed Δ<sub>48</sub> value(s) 1114 * `kinetic_slope`: kinetic fractionation slopw, with or without uncertainty 1115 1116 Returns a tuple of uarrays corresponding to the projected Δ<sub>47</sub> and Δ<sub>48</sub> values. 1117 1118 > [!NOTE] 1119 > This is not a least-squares minimization problem but a direct calculation, and should thus 1120 > be much faster than the various `CorelData.nearestD47eq()` methods. 1121 """ 1122 1123 D47 = _cd.uarray(D47) 1124 D48 = _cd.uarray(D48) 1125 N = D47.size 1126 N47c = self.D47_coefs.size 1127 N48c = self.D48_coefs.size 1128 D47p = D47 * 0 1129 1130 for i in range(N): 1131 1132 # function to solve 1133 def fun(x, *args): # args = (D47, D48, kinetic_slope, *self.D47_coefs, *self.D48_coefs) 1134 1135 args = _np.array(args) 1136 D47_n = args[0] 1137 D48_n = args[1] 1138 kslope_n = args[2] 1139 D47_calib_coefs_n = args[-N48c-N47c:-N48c] 1140 D48_calib_coefs_n = args[-N48c:] 1141 1142 D47i = D4x_calib_function( 1143 self.interp.T, 1144 D47_calib_coefs_n, 1145 return_without_uncertainties = False, 1146 ) 1147 1148 D48i = D4x_calib_function( 1149 self.interp.T, 1150 D48_calib_coefs_n, 1151 return_without_uncertainties = False, 1152 ) 1153 1154 D48_interp = uarray_compatible_interp(D47i, D48i) 1155 1156 return D48_n - D48_interp(x) - kslope_n * (D47_n - x) 1157 1158 def g(*args): 1159 return _fsolve(fun, [100.], args = args)[0] 1160 1161 wg = _uc.wrap(g) 1162 1163 D47p[i] = wg( 1164 D47[i], 1165 D48[i], 1166 kinetic_slope, 1167 *self.D47_coefs, 1168 *self.D48_coefs, 1169 ) 1170 1171 _, D48p = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47p, ignore_calib_uncertainties = False) 1172 1173 return D47p, D48p 1174 1175 def Teq_pdf( 1176 self, 1177 D47: _uc.ufloat, 1178 Tmin: (float | None) = None, 1179 Tmax: (float | None) = None, 1180 Tinc: float = 0.2, 1181 default_D47_sigmas: float = 4.0, 1182 ignore_calib_uncertainties: bool = False, 1183 run_qmc: bool = False, 1184 N_qmc: int = 1024, 1185 ): 1186 """ 1187 Compute the unit-normalized probability distribution function (PDF) of the 1188 equilibrium temperature (`Teq`) for a given (`UFloat`) value of Δ<sub>47</sub>. 1189 1190 **Arguments** 1191 * `D47`: Δ<sub>47</sub> value (with uncertainty) 1192 * `Tmin`: minimum temperature over which to compute the PDF; if not specified, 1193 use temperature corresponding to `D47.n + `default_D47_sigmas` * D47.s` 1194 * `Tmax`: maximum temperature over which to compute the PDF; if not specified, 1195 use temperature corresponding to `D47.n - `default_D47_sigmas` * D47.s` 1196 * `Tinc`: temperature increment over which to compute the PDF 1197 * `default_D47_sigmas`: see `Tmin` and `Tmin` above 1198 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 1199 * `run_qmc`: whether to also run a Quasi Monte carlo simulation to estimate the PDF 1200 * `N_qmc`: number of iterations in the above Quasi Monte Carlo simulation 1201 1202 **Returns** 1203 * `Ti`: Evenly-spaced array of temperature values over which the PDF is computed 1204 * `pdf`: PDF evaluated over `Ti` 1205 * `Tqmc` (only returned if `run_qmc = True`): array of `N_qmc` temperature values 1206 computed in the Quasi Monte Carlo simulation 1207 """ 1208 1209 if Tmin is None: 1210 Tmin = _np.floor(self.T_as_function_of_D47( 1211 D47.n + default_D47_sigmas * D47.s, 1212 ignore_calib_uncertainties = ignore_calib_uncertainties, 1213 ).n) 1214 1215 if Tmax is None: 1216 Tmax = _np.ceil(self.T_as_function_of_D47( 1217 D47.n - default_D47_sigmas * D47.s, 1218 ignore_calib_uncertainties = ignore_calib_uncertainties, 1219 ).n) 1220 1221 assert Tmin < Tmax, "Tmax must be strictly greater than Tmin" 1222 assert Tinc > 0, "Tinc must be strictly greater than zero" 1223 1224 # compute interpolated Ti values 1225 Ti = _np.arange(Tmin, Tmax+Tinc, Tinc) 1226 1227 pdf = transform_pdf_monotonic( 1228 f_inv = lambda T: D4x_calib_function( 1229 T, 1230 self.D47_coefs, 1231 return_without_uncertainties = ignore_calib_uncertainties, 1232 ignore_calib_uncertainties = ignore_calib_uncertainties, 1233 ), 1234 df_inv = lambda T: D4x_calib_derivative( 1235 T, 1236 self.D47_coefs, 1237 return_without_uncertainties = ignore_calib_uncertainties, 1238 ignore_calib_uncertainties = ignore_calib_uncertainties, 1239 ), 1240 mu_x = D47.n, 1241 sigma_x = D47.s, 1242 yi = Ti, 1243 ) 1244 1245 if run_qmc: 1246 1247 from scipy.stats import qmc 1248 from tqdm.rich import tqdm 1249 1250 #parameters to jiggle 1251 input_params = _cd.uarray([D47, *self.D47_coefs]) 1252 1253 # QMC sampler for the correlation matrix of these parameters 1254 qmc_dist = qmc.MultivariateNormalQMC( 1255 mean = input_params.n*0, 1256 cov = input_params.cor, 1257 ) 1258 1259 # QMC samples 1260 qmc_draws = input_params.n + qmc_dist.random(N_qmc) * input_params.s 1261 1262 # initialize T_qmc 1263 Tqmc = _cd.uarray(_np.zeros((N_qmc,))) 1264 1265 for k in tqdm(range(N_qmc)): 1266 # jiggled D47 and D47coefs 1267 _D47 = qmc_draws[k,0] 1268 if ignore_calib_uncertainties: 1269 _coefs = self.D47_coefs 1270 else: 1271 _coefs = _cd.uarray(_uc.correlated_values(qmc_draws[k,1:], self.D47_coefs.covar)) 1272 1273 # jiggled D47 1274 _D47i = D4x_calib_function(self.interp.T, _coefs) 1275 _f = uarray_compatible_interp(_D47i.n, self.interp.T) 1276 Tqmc[k] = _f(_D47) 1277 1278 return Ti, pdf, Tqmc 1279 1280 return Ti, pdf 1281 1282 1283### Utilities and CLI ### 1284 1285 1286def save_Teq_report( 1287 X, 1288 Y, 1289 T, 1290 p, 1291 filename, 1292 Xname = 'D47', 1293 Yname = 'D48', 1294 Tname = 'T95', 1295 labelname = 'Sample', 1296 fmt_Xnv = '.4f', 1297 fmt_Xse = '.4f', 1298 fmt_Ynv = '.4f', 1299 fmt_Yse = '.4f', 1300 fmt_Tnv = '.1f', 1301 fmt_Tse = '.1f', 1302 fmt_cm = '.6f', 1303 fmt_pv = '.2e', 1304 labels = None, 1305 sep = ',', 1306 p_cutoff = 0.05, 1307): 1308 """ 1309 Save a temperature report to a csv file. 1310 Includes observed `D47`, `D48`, p-equilibrium values, and nearest `Teq` with sensible precision defaults. 1311 Alternatively, users may find [`correldata.CorrelData.str()`](https://mdaeron.github.io/correldata/#CorrelData.str) 1312 to be more versatile. 1313 """ 1314 N = T.size 1315 if labels is None: 1316 labels = [str(k+1) for k in range(N)] 1317 1318 with open(filename, 'w') as fid: 1319 fid.write(f'{labelname}{sep}{Xname}{sep}SE{sep}correl{sep*N}{Yname}{sep}SE{sep}correl{sep*N}p-value{sep}{Tname}{sep}SE{sep}correl') 1320 Xnv = _unp.nominal_values(X) 1321 Xse = _unp.std_devs(X) 1322 Xcm = _np.array(_uc.correlation_matrix(X)) 1323 Ynv = _unp.nominal_values(Y) 1324 Yse = _unp.std_devs(Y) 1325 Ycm = _np.array(_uc.correlation_matrix(Y)) 1326 Tnv = _unp.nominal_values(T) 1327 Tse = _unp.std_devs(T) 1328 Tcm = _np.array(_uc.correlation_matrix(T)) 1329 for k in range(X.size): 1330 fid.write(f'\n{labels[k]}{sep}{Xnv[k]:{fmt_Xnv}}{sep}{Xse[k]:{fmt_Xse}}{sep}') 1331 fid.write(sep.join([f'{Xcm[j,k]:{fmt_cm}}' for j in range(N)])) 1332 fid.write(f'{sep}{Ynv[k]:{fmt_Ynv}}{sep}{Yse[k]:{fmt_Yse}}{sep}') 1333 fid.write(sep.join([f'{Ycm[j,k]:{fmt_cm}}' for j in range(N)])) 1334 fid.write(f'{sep}{p[k]:{fmt_pv}}') 1335 if p[k] >= p_cutoff: 1336 fid.write(f'{sep}{Tnv[k]:{fmt_Tnv}}{sep}{Tse[k]:{fmt_Tse}}{sep}') 1337 fid.write(sep.join([f'{Tcm[j,k]:{fmt_cm}}' for j in range(N)])) 1338 1339_rich_utils.STYLE_HELPTEXT = '' 1340 1341__app = _typer.Typer( 1342 add_completion = False, 1343 context_settings={'help_option_names': ['-h', '--help']}, 1344 rich_markup_mode = 'rich', 1345) 1346 1347@__app.command() 1348def _cli( 1349 input: _Annotated[str, _typer.Option('--input', '-i', help = "Input file to read from (otherwise read from stdin).")] = None, 1350 output: _Annotated[str, _typer.Option('--output', '-o', help = "Output file to write to (otherwise write to stdout).")] = None, 1351 kslope: _Annotated[str, _typer.Option('--kslope', '-k', help = "Kinetic fractionation slope, using format [bold]'n(s)'[/bold] (with quotes), where [bold]n[/bold] is the slope and [bold]s[/bold] its standard error.")] = None, 1352 hpoutput: _Annotated[bool, _typer.Option('--high-precision-output', '-p', help = "Generate higher precision output.")] = False, 1353 show_mixed_correl: _Annotated[bool, _typer.Option('--show_mixed_correl', '-m', help = "Show correlations between different fields.")] = False, 1354 version: _Annotated[bool, _typer.Option('--version', '-v', help = 'Show version and exit.')] = False, 1355): 1356 """ 1357[b]Purpose:[/b] 1358 1359Reads data from an input file, computes p-value and T estimates, and print out the results. 1360""" 1361 if version: 1362 print(__version__) 1363 return None 1364 1365 if input is None: 1366 datastring = ''.join(sys.stdin) 1367 elif isinstance(input, str): 1368 with open(input) as fid: 1369 datastring = fid.read() 1370 1371 data = _cd.read_str(datastring) 1372 1373 E = Engine() 1374 1375 D47eq, D48eq, p = E.nearest_D47eq(data['D47'], data['D48']) 1376 Teq = E.T_as_function_of_D47(D47eq) 1377 data['eq_pvalue'] = p 1378 data['Teq'] = Teq 1379 1380 if isinstance(kslope, str): 1381 kslope = kslope.split(')')[0] 1382 kslope = kslope.split('(') 1383 kslope = _uc.ufloat(float(kslope[0]), float(kslope[1])) 1384 1385 D47kp, D48kp = E.projected_D47eq(data['D47'], data['D48'], kinetic_slope = kslope) 1386 Tkp = E.T_as_function_of_D47(D47kp) 1387 1388 data['kslope'] = _cd.uarray([kslope for _ in data['D47']]) 1389 1390 data['Tkp'] = Tkp 1391 1392 ffmt = { 1393 'D47': '.6f', 1394 'D48': '.6f', 1395 'kslope': lambda x: f'{x:z.6f}'.rstrip('0'), 1396 'Teq': 'z.6f', 1397 'Tkp': 'z.6f', 1398 } if hpoutput else { 1399 'D47': '.4f', 1400 'D48': '.4f', 1401 'kslope': lambda x: f'{x:z.6f}'.rstrip('0'), 1402 'Teq': 'z.2f', 1403 'Tkp': 'z.2f', 1404 } 1405 1406 out = data.str( 1407 float_format = ffmt, 1408 show_mixed_correl = show_mixed_correl, 1409 exclude_fields = ['correl_kslope'], 1410 ) 1411 1412 if output is None: 1413 print(out) 1414 elif isinstance(output, str): 1415 with open(output, 'w') as fid: 1416 fid.write(out) 1417 1418def __cli(): __app()
def
ufloat_compatible_interp( xi: correldata.uarray | ArrayLike, yi: correldata.uarray | ArrayLike, x: float | uncertainties.core.AffineScalarFunc | correldata.uarray | ArrayLike):
47def ufloat_compatible_interp( 48 xi: (_cd.uarray | ArrayLike), 49 yi: (_cd.uarray | ArrayLike), 50 x: (float | _uc.UFloat | _cd.uarray | ArrayLike), 51): 52 """ 53 Linear interpolation accepting UFloat values for all three input parameters. 54 Only handles one interpolated value. For interpolated arrays, use `uarray_compatible_interp()` 55 56 **Arguments** 57 * `xi`: x-values defining the interpolated function 58 * `yi`: y-values defining the interpolated function 59 * `x`: x-value of the interpolation point 60 61 Returns y-value of the interpolation point, either as a float or a UFloat. 62 """ 63 xn = x.nominal_value if isinstance(x, _uc.UFloat) else float(x) 64 idx = _np.searchsorted(xi, xn) 65 idx = _np.clip(idx, 1, len(xi) - 1) 66 67 x0 = xi[idx-1] 68 x1 = xi[idx] 69 y0 = yi[idx-1] 70 y1 = yi[idx] 71 72 t = (x - x0) / (x1 - x0) 73 return y0 + t * (y1 - y0)
Linear interpolation accepting UFloat values for all three input parameters.
Only handles one interpolated value. For interpolated arrays, use uarray_compatible_interp()
Arguments
xi: x-values defining the interpolated functionyi: y-values defining the interpolated functionx: x-value of the interpolation point
Returns y-value of the interpolation point, either as a float or a UFloat.
def
uarray_compatible_interp(xi, yi):
76def uarray_compatible_interp(xi, yi): 77 """ 78 Linear interpolation accepting UFloat values for all three input parameters. 79 80 **Arguments** 81 * `xi`: x-values defining the interpolated function 82 * `yi`: y-values defining the interpolated function 83 84 Returns an interpolation function which returns arrays or uarrays of y-values. 85 """ 86 return _np.vectorize( 87 lambda x: ufloat_compatible_interp(xi, yi, x) 88 )
Linear interpolation accepting UFloat values for all three input parameters.
Arguments
xi: x-values defining the interpolated functionyi: y-values defining the interpolated function
Returns an interpolation function which returns arrays or uarrays of y-values.
def
transform_pdf_monotonic(f_inv, df_inv, mu_x, sigma_x, yi):
91def transform_pdf_monotonic(f_inv, df_inv, mu_x, sigma_x, yi): 92 """ 93 Compute probability distribution function of Y = f(X) 94 where X ~ Normal(mu_x, sigma_x) and f is monotonic, 95 based on the change-of-variables formula: 96 97 p[y=f(x)] = p[x=f_inv(y)] * d(f_inv)/dy 98 99 Additionally, if f_inv returns UFloats, the PDF is convolved with that local 100 source of uncertainty (assumed to be Gaussian) at each grid point. 101 102 As currently implemented, requires `yi` to be an equally spaced array-like. 103 104 **Arguments** 105 f_inv: inverse of f, may return UFloats 106 df_inv: derivative of f_inv, should return UFloats if f_inv does 107 mu_x: mean of X PDF 108 sigma_x: std dev of X PDF 109 yi: regularly spaced grid of y values at which to evaluate the PDF 110 111 **Returns:** 112 pdf: normalized PDF evaluated at yi 113 """ 114 115 if not _np.allclose(_np.diff(yi), yi[1] - yi[0]): 116 raise ValueError("yi must be regularly spaced") 117 118 xi = f_inv(yi) # may be floats or ufloats, depending on f_inv 119 120 try: 121 xi_nom = xi.n 122 sigma_xi = xi.s 123 has_ufloats = True 124 except AttributeError: 125 xi_nom = xi 126 has_ufloats = False 127 128 # Jacobian weights (account for irregular xi spacing) 129 try: 130 df_inv_nom = df_inv(yi).n 131 except AttributeError: 132 df_inv_nom = df_inv(yi) 133 134 w_i = _norm.pdf(xi_nom, loc = mu_x, scale = sigma_x) * _np.abs(df_inv_nom) 135 136 if not has_ufloats: 137 return w_i / (_np.trapezoid(w_i, yi)) 138 139 # Propagate sigma from x-space to y-space via Jacobian: sigma_y = sigma_x / abs( dx/dy ) 140 sigma_yi = sigma_xi / _np.abs(df_inv_nom) 141 142 # Convolution of Gaussians: each grid point j contributes N(yi; yj, σ_yj²) scaled by w_j 143 gaussians = _norm.pdf( 144 yi[:, None], 145 loc = yi[None, :], 146 scale = sigma_yi[None, :] 147 ) # NOTE: nice syntax to reshape ndarrays, perhaps use this in D4x_calib_function? 148 149 pdf = (gaussians * w_i[None, :]).sum(axis = 1) 150 151 return pdf / (_np.trapezoid(pdf, yi))
Compute probability distribution function of Y = f(X) where X ~ Normal(mu_x, sigma_x) and f is monotonic, based on the change-of-variables formula:
p[y=f(x)] = p[x=f_inv(y)] * d(f_inv)/dy
Additionally, if f_inv returns UFloats, the PDF is convolved with that local source of uncertainty (assumed to be Gaussian) at each grid point.
As currently implemented, requires yi to be an equally spaced array-like.
Arguments f_inv: inverse of f, may return UFloats df_inv: derivative of f_inv, should return UFloats if f_inv does mu_x: mean of X PDF sigma_x: std dev of X PDF yi: regularly spaced grid of y values at which to evaluate the PDF
Returns: pdf: normalized PDF evaluated at yi
def
D4x_calib_function( T: float | uncertainties.core.AffineScalarFunc | correldata.uarray | ArrayLike, coefs: correldata.uarray, return_without_uncertainties: bool = False, ignore_calib_uncertainties: bool = False) -> float | uncertainties.core.AffineScalarFunc | correldata.uarray | ArrayLike:
251def D4x_calib_function( 252 T: (float | _uc.UFloat | _cd.uarray | ArrayLike), 253 coefs: _cd.uarray, 254 return_without_uncertainties: bool = False, 255 ignore_calib_uncertainties: bool = False, 256) -> (float | _uc.UFloat | _cd.uarray | ArrayLike): 257 """ 258 **Arguments** 259 * `T`: temperature(s) for which to compute Δ<sub>4x</sub> 260 * `return_without_uncertainties`: if `True`, returns Δ<sub>4x</sub> values without error propagation of any kind 261 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 262 263 Returns equilibrium Δ<sub>4x</sub> value(s) corresponding to `T` value(s) 264 """ 265 degs = _np.arange(coefs.size) 266 267 D4x = ( 268 _np.expand_dims(_cd.nv(coefs) if ignore_calib_uncertainties else coefs, 1) # shape = (coefs.size, 1) 269 * _np.expand_dims((T+273.15)**-1, 0) # shape = (1, T.size) 270 ** _np.expand_dims(degs, 1) # shape = (coefs.size, 1) 271 ).sum(axis = 0 if isinstance(T, _np.ndarray) else None) 272 273 if D4x.ndim == 0: 274 return D4x.tolist().n if return_without_uncertainties else D4x.tolist() 275 return D4x.n if return_without_uncertainties else D4x
Arguments
T: temperature(s) for which to compute Δ4xreturn_without_uncertainties: ifTrue, returns Δ4x values without error propagation of any kindignore_calib_uncertainties: whether to propagate calibration uncertainties
Returns equilibrium Δ4x value(s) corresponding to T value(s)
def
D4x_calib_derivative( T: float | uncertainties.core.AffineScalarFunc | correldata.uarray | ArrayLike, coefs: correldata.uarray, return_without_uncertainties: bool = False, ignore_calib_uncertainties: bool = False) -> float | uncertainties.core.AffineScalarFunc | correldata.uarray | ArrayLike:
278def D4x_calib_derivative( 279 T: (float | _uc.UFloat | _cd.uarray | ArrayLike), 280 coefs: _cd.uarray, 281 return_without_uncertainties: bool = False, 282 ignore_calib_uncertainties: bool = False, 283) -> (float | _uc.UFloat | _cd.uarray | ArrayLike): 284 """ 285 **Arguments** 286 * `T`: temperature(s) for which to compute Δ<sub>4x</sub> 287 * `return_without_uncertainties`: if `True`, returns D4x values without error propagation of any kind. 288 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties. 289 290 Returns d(D4x)/dT corresponding to `T` value(s) 291 """ 292 dcoefs = -_np.arange(len(coefs)) * coefs 293 dcoefs = _cd.uarray((dcoefs[0], *dcoefs)) 294 return D4x_calib_function( 295 T, 296 dcoefs, 297 return_without_uncertainties = return_without_uncertainties, 298 ignore_calib_uncertainties = ignore_calib_uncertainties, 299 )
Arguments
T: temperature(s) for which to compute Δ4xreturn_without_uncertainties: ifTrue, returns D4x values without error propagation of any kind.ignore_calib_uncertainties: whether to propagate calibration uncertainties.
Returns d(D4x)/dT corresponding to T value(s)
def
conf_ellipse( X: correldata.uarray | numpy.ndarray | uncertainties.core.AffineScalarFunc | float, Y: correldata.uarray | numpy.ndarray | uncertainties.core.AffineScalarFunc | float = None, p: float = 0.95, CM: numpy.ndarray | None = None, Xse: numpy.ndarray | float | None = None, Yse: numpy.ndarray | float | None = None, plot: bool = True, ax: matplotlib.axes._axes.Axes | None = None, **kwargs) -> tuple:
305def conf_ellipse( 306 X: (_cd.uarray | _np.ndarray | _uc.UFloat | float), 307 Y: (_cd.uarray | _np.ndarray | _uc.UFloat | float) = None, 308 p: float = 0.95, 309 CM: (_np.ndarray | None) = None, 310 Xse: (_np.ndarray | float | None) = None, 311 Yse: (_np.ndarray | float | None) = None, 312 plot: bool = True, 313 ax: (_ppl.Axes | None) = None, 314 **kwargs, 315) -> tuple: 316 """ 317 Compute and (optionally) plot the joint *p*-level confidence ellipses for the elements of (X, Y) 318 319 **Arguments** 320 * `X`: x values 321 * `Y`: y values 322 * `p`: confidence level 323 * `CM`: covariance matrix of (X, Y); not needed if X and Y are of type 324 [`uncertainties.UFloat`](https://pythonhosted.org/uncertainties/tech_guide.html). 325 or if (`Xse`, `Yse`) are specified. 326 * `Xse`, `Yse`: SE of X and Y; not needed if X and Y are of type 327 [`uncertainties.UFloat`](https://pythonhosted.org/uncertainties/tech_guide.html) 328 or if `CM` is specified. 329 * `plot`: whether to plot the ellipse or not. If `False`, return a list of 330 `(x_center, y_center, width, height, angle)` elements 331 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 332 * `kwargs`: passed to `matplotlib.patches.Ellipse()` 333 334 Returns a list of the `Ellipse` objects thus created. 335 """ 336 337 r2 = _chi2.ppf(p, 2) 338 kwargs = dict(fc = 'None', ec = 'k', lw = 0.7) | kwargs 339 340 out = [] 341 342 for x, y in zip( 343 *_cd.as_pair_of_uarrays(X, Y, CM = CM, Xse = Xse, Yse = Yse) 344 ): 345 val, vec = _eigh(_uc.covariance_matrix((x, y))) 346 width, height = 2 * (val[:, None] * r2)**0.5 347 angle = _np.degrees(_np.arctan2(*vec[::-1, 0])) 348 349 if plot: 350 from matplotlib import pyplot as _ppl 351 from matplotlib.patches import Ellipse as _Ellipse 352 353 if ax is None: 354 ax = _ppl.gca() 355 356 out.append( 357 ax.add_patch( 358 _Ellipse( 359 xy = (x.n, y.n), 360 width = width.item(), 361 height = height.item(), 362 angle = angle, 363 **kwargs, 364 ) 365 ) 366 ) 367 else: 368 out.append([x.n, y.n, width, height, angle]) 369 370 return (*out,)
Compute and (optionally) plot the joint p-level confidence ellipses for the elements of (X, Y)
Arguments
X: x valuesY: y valuesp: confidence levelCM: covariance matrix of (X, Y); not needed if X and Y are of typeuncertainties.UFloat. or if (Xse,Yse) are specified.Xse,Yse: SE of X and Y; not needed if X and Y are of typeuncertainties.UFloator ifCMis specified.plot: whether to plot the ellipse or not. IfFalse, return a list of(x_center, y_center, width, height, angle)elementsax: which instance ofmatplotlib.axes.Axesto draw in; use current axes ifax=None.kwargs: passed tomatplotlib.patches.Ellipse()
Returns a list of the Ellipse objects thus created.
class
Engine:
378class Engine(): 379 """ 380 Underlying engine to compute and plot nearest equilibrium temperatures and projected 381 temperatures based on a consistent pair of Δ<sub>47</sub>, Δ<sub>48</sub> calibrations. 382 """ 383 384 # D47_calib_coefs from OGLS23 (D47calib v1.3.1) 385 D47_calib_coefs = _cd.read_str(''' 386 coefs, SE, correl, 3870.17437754366432887, 4.911105567257293e-3, 1. , -0.93797005, 0.8865771 388 -18.14215245127414, 5.632326472234856, -0.93797005, 1. , -0.98994249 38942.65722989162373e3, 1.27712751715908e3, 0.8865771 , -0.98994249, 1. 390'''[1:-1])['coefs'] 391 """ 392 Default (OGLS23) Δ<sub>47</sub> calibration coefficients based on [Daëron & Vermeesch (2024)](https://doi.org/10.1016/j.chemgeo.2023.121881) 393 """ 394 395 # D48_calib_coefs reprocessed from Fiebig et al. (2024): 396 # 397 # D48_calib_coefs = _compute_D48_calib_coefficients(reprocess = True) 398 # print(_cd.data_string( 399 # {'coefs': D48_calib_coefs}, 400 # float_format = 'z.12g', 401 # correl_format = 'z.12f', 402 # )) 403 404 D48_calib_coefs = _cd.read_str(''' 405 coefs, SE_coefs, correl_coefs, , , , 4060.121349237888, 0.00390048540724, 1.000000000000, -0.664181963395, 0.664181963395, -0.664181963395, 0.664181963395 407 6.22931985613, 0.32896761459, -0.664181963395, 1.000000000000, -1.000000000000, 1.000000000000, -1.000000000000 408 -13481.983494, 711.977559735, 0.664181963395, -1.000000000000, 1.000000000000, -1.000000000000, 1.000000000000 409 9336714.66607, 493067.754224, -0.664181963395, 1.000000000000, -1.000000000000, 1.000000000000, -1.000000000000 410-770413883.573, 40685214.9801, 0.664181963395, -1.000000000000, 1.000000000000, -1.000000000000, 1.000000000000 411'''[1:-1])['coefs'] 412 """ 413 Default Δ<sub>48</sub> calibration coefficients based on [Fiebig et al. (2024)](https://doi.org/10.1016/j.chemgeo.2024.122382) 414 """ 415 416 def __init__( 417 self, 418 D47_coefs: (_cd.uarray | ArrayLike | None) = None, 419 D48_coefs: (_cd.uarray | ArrayLike | None) = None, 420 Tmin_interp: float = -23.0, 421 Tmax_interp: float = 1277.0, 422 N_interp: float = 201, 423 ): 424 """ 425 **Arguments** 426 * `D47_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 427 * `D48_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 428 * `Tmin_interp`: minimum temperature over which to interpolate for inverse function computations 429 * `Tmax_interp`: maximum temperature over which to interpolate for inverse function computations 430 * `N_interp`: number of points (equally-spaced in 1/T space) over which to interpolate for inverse function computations 431 """ 432 433 self.D47_coefs = Engine.D47_calib_coefs if D47_coefs is None else D47_coefs 434 """The Δ<sub>47</sub> calibration coefficients used by this `Engine` instance""" 435 436 self.D48_coefs = Engine.D48_calib_coefs if D48_coefs is None else D48_coefs 437 """The Δ<sub>48</sub> calibration coefficients used by this `Engine` instance""" 438 439 self.interp = _Interpolation() 440 """ 441 Holds equilibrium Δ<sub>47</sub> and Δ<sub>48</sub> values (ufloats) interpolated 442 along an array of T values (regularly spaced increments of 1/T<sup>2</sup>). 443 444 * `interp.T`: interpolation T values (floats) in regularly spaced increments of 1/T<sup>2</sup> 445 * `interp.D47`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T` 446 * `interp.D48`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T` 447 * `interp.D47_no_calib_errors`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T`, 448 ignoring calibration uncertainties 449 * `interp.D48_no_calib_errors`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T`, 450 ignoring calibration uncertainties 451 """ 452 453 self.interp.T = _np.linspace( 454 (Tmax_interp+273.15)**-2, 455 (Tmin_interp+273.15)**-2, 456 N_interp, 457 )**-0.5 - 273.15 458 459 self.interp.D47 = self.D47_calib_function( 460 self.interp.T, 461 return_without_uncertainties = False, 462 ignore_calib_uncertainties = False, 463 ) 464 465 self.interp.D47_no_calib_errors = self.D47_calib_function( 466 self.interp.T, 467 return_without_uncertainties = False, 468 ignore_calib_uncertainties = True, 469 ) 470 471 self.interp.D48 = self.D48_calib_function( 472 self.interp.T, 473 return_without_uncertainties = False, 474 ignore_calib_uncertainties = False, 475 ) 476 477 self.interp.D48_no_calib_errors = self.D48_calib_function( 478 self.interp.T, 479 return_without_uncertainties = False, 480 ignore_calib_uncertainties = True, 481 ) 482 483 self.interp.D47u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D47) 484 self.interp.D48u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D48) 485 486 #inverse D47 calibration (ignoring calibration errors) 487 self.interp.Teq_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.T) 488 #inverse D47 calibration (including calibration errors) 489 self.interp.Teq_as_function_of_D47u = uarray_compatible_interp(self.interp.D47, self.interp.T) 490 491 def T_as_function_of_D47( 492 self, 493 D47: (_cd.uarray | ArrayLike), 494 ignore_calib_uncertainties: bool = False, 495 ): 496 """ 497 Provided with one or more Δ<sub>47</sub> values (floats or ufloats), return ufloats for the 498 corresponding equilibrium T values (ufloats with or without Δ<sub>47</sub> calibration uncertainties). 499 500 **Arguments** 501 * `D47`: array of Δ<sub>47</sub> values 502 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 503 """ 504 if ignore_calib_uncertainties: 505 return _cd.uarray(self.interp.Teq_as_function_of_D47n(D47)) 506 else: 507 return _cd.uarray(self.interp.Teq_as_function_of_D47u(D47)) 508 509 def D47u_as_function_of_D47n( 510 self, 511 D47: ArrayLike 512 ): 513 """ 514 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 515 equilibrium Δ<sub>47</sub> values (ufloats with Δ<sub>47</sub> calibration uncertainties). 516 """ 517 return _cd.uarray(self.interp.D47u_as_function_of_D47n(D47)) 518 519 def D48u_as_function_of_D47n( 520 self, 521 D47: ArrayLike 522 ): 523 """ 524 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 525 equilibrium Δ<sub>48</sub> values (ufloats with Δ<sub>48</sub> calibration uncertainties). 526 """ 527 return _cd.uarray(self.interp.D48u_as_function_of_D47n(D47)) 528 529 def D47_calib_function( 530 self, 531 T: (float | _uc.UFloat | _cd.uarray), 532 return_without_uncertainties: bool = False, 533 ignore_calib_uncertainties: bool = False, 534 ): 535 return D4x_calib_function( 536 T = T, 537 coefs = self.D47_coefs, 538 return_without_uncertainties = return_without_uncertainties, 539 ignore_calib_uncertainties = ignore_calib_uncertainties, 540 ) 541 542 def D48_calib_function( 543 self, 544 T: (float | _uc.UFloat | _cd.uarray), 545 return_without_uncertainties: bool = False, 546 ignore_calib_uncertainties: bool = False, 547 ): 548 return D4x_calib_function( 549 T = T, 550 coefs = self.D48_coefs, 551 return_without_uncertainties = return_without_uncertainties, 552 ignore_calib_uncertainties = ignore_calib_uncertainties, 553 ) 554 555 D47_calib_function.__doc__ = D4x_calib_function.__doc__.replace('Δ<sub>4x</sub>', 'Δ<sub>47</sub>') 556 D48_calib_function.__doc__ = D4x_calib_function.__doc__.replace('Δ<sub>4x</sub>', 'Δ<sub>48</sub>') 557 558 def T_ellipse( 559 self, 560 T: (_np.ndarray | _cd.uarray), 561 p: float = 0.95, 562 CM: (_np.ndarray | None) = None, 563 Tse: (_np.ndarray | float | None) = None, 564 plot: bool = True, 565 ax: (_ppl.Axes | None) = None, 566 **kwargs, 567 ) -> list: 568 """ 569 Plot the joint `p`-level confidence ellipses in (Δ<sub>47</sub>, Δ<sub>48</sub>) 570 space, for temperatures equal to the elements of `T`, and return a list of the 571 `Ellipse` objects thus created. 572 573 **Arguments** 574 * `T`: `ndarray` or `uarray` of temperatures to plot 575 * `p`: confidence level 576 * `plot`: whether to plot the ellipse or not. If `False`, return a list of 577 `(x_center, y_center, width, height, angle)` elements 578 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 579 * `kwargs`: passed to `matplotlib.patches.Ellipse()` 580 """ 581 _T = _cd.as_uarray(T, CM = CM, Xse = Tse) 582 return conf_ellipse( 583 self.D47_calib_function(_T), 584 self.D48_calib_function(_T), 585 p = p, 586 plot = plot, 587 ax = ax, 588 **kwargs, 589 ) 590 591 def plot_D95_confidence_band( 592 self, 593 p: float = 0.95, 594 Ti: (ArrayLike | None) = None, 595 plot: bool = True, 596 ax: (_ppl.Axes | None) = None, 597 **kwargs, 598 ): 599 """ 600 Plot, for a given p-value, the confidence band of the thermodynamic equilibrium curve 601 in (Δ<sub>47</sub>, Δ<sub>48</sub>) space. 602 603 **Arguments** 604 * `p`: confidence level 605 * `Ti`: array of temperatures over which to evaluate confidence band (default: use `interp.T` attribute instead) 606 * `plot`: whether to plot the confidence band or not. If `False`, return the (N,2) array of polygon nodes 607 * `ax`: `Axes` instance to plot to (default: use current Axes) 608 * `kwargs`: passed to `patches.Polygon()` 609 610 Returns the corresponding `Polygon` instance. 611 """ 612 613 if Ti is None: 614 Ti = self.interp.T 615 616 cb = confidence_band( 617 Ti, 618 self.D47_calib_function, 619 self.D48_calib_function, 620 p, 621 ) 622 623 if plot: 624 from matplotlib import pyplot as _ppl 625 from matplotlib.patches import Polygon as _Polygon 626 627 if ax is None: 628 ax = _ppl.gca() 629 630 polygon = ax.add_patch( 631 _Polygon( 632 cb, 633 closed = True, 634 **kwargs, 635 ) 636 ) 637 return polygon 638 else: 639 return cb 640 641 642 def plot_D95_equilibrium( 643 self, 644 Tmin: float = 0., 645 Tmax: float = 1000., 646 NT: int = 101, 647 Tmarkers: _np.typing.ArrayLike = [0, 25, 100, 250, 1000], 648 kwargs_Tmarkers: dict = {}, 649 show_Tmarker_labels: bool = True, 650 kwargs_Tmarker_labels: dict = {}, 651 show_Tmarker_ellipses: bool = False, 652 kwargs_Tmarker_ellipses: dict = {}, 653 show_eqline: bool = True, 654 kwargs_eqline: dict = {}, 655 show_confidence: bool = True, 656 confidence_pvalue: float = 0.95, 657 kwargs_confidence: dict = {}, 658 ax: (_ppl.Axes | None) = None, 659 xlabel: str = '$Δ_{47}$ [‰]', 660 ylabel: str = '$Δ_{48}$ [‰]', 661 lw: float = 0.7, 662 ) -> (dict, dict): 663 """ 664 Plot a thermodynamic equilibrium curve in (Δ<sub>47</sub>, Δ<sub>48</sub>) space 665 as a function of temperature. 666 667 **Arguments** 668 * `Tmin`: minimum T to plot 669 * `Tmax`: maximum T to plot 670 * `NT`: number of steps in equilibrium curve (interpolated at constant steps in 1/T<sup>2</sup> space) 671 * `Tmarkers`: T markers to add along the curve 672 * `kwargs_Tmarkers`: passed to `plot()` when plotting T markers 673 * `show_Tmarker_labels`: whether to add T labels to T markers 674 * `kwargs_Tmarker_labels`: passed to `text()` when plotting T markers 675 * `show_Tmarker_ellipses`: whether to add confidence ellipses to T markers 676 * `kwargs_Tmarker_ellipses`: passed to `T_ellipses()` when plotting T marker ellipses 677 * `show_eqline`: whether to plot the equilibrium curve itself 678 * `kwargs_eqline`: passed to `plot()` when plotting the equilibrium curve 679 * `show_confidence`: whether to plot the confidence band of the equilibrium curve 680 * `confidence_pvalue`: confidence level for the confidence band 681 * `kwargs_confidence`: passed to `plot_D95_confidence_band()` when plotting the confidence band 682 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 683 * `xlabel`: string to pass to `xlabel()` 684 * `ylabel`: string to pass to `ylabel()` 685 * `lw`: default line width for most plot elements 686 687 **Returns** 688 * `data`: a dict of the T, Δ<sub>47</sub> and Δ<sub>48</sub> values generated for this plot: 689 - `Te` : temperature interpolated along the equilibrium curve 690 - `D47e`: Δ<sub>47</sub> interpolated along the equilibrium curve 691 - `D48e`: Δ<sub>48</sub> interpolated along the equilibrium curve 692 - `Tm` : temperature of T markers 693 - `D47m`: Δ<sub>47</sub> of T markers 694 - `D48m`: Δ<sub>48</sub> of T markers 695 696 * `plot_elements`: a dict of the `Axes` elements generated for this plot: 697 - `eqline`: `Line2D` of the equilibrium curve 698 - `confidence`: `Polygon` object for the confidence band 699 - `Tm`: `Line2D` of the T markers 700 - `Tme`: list of `Ellipse` objects for the T marker ellipses 701 - `Tml`: list of `Text` objects for the T marker labels 702 """ 703 704 from matplotlib import pyplot as _ppl 705 706 default_kwargs_eqline = dict( 707 marker = 'None', 708 ls = '-', 709 color = 'k', 710 lw = lw, 711 ) 712 default_kwargs_confidence = dict( 713 ec = (0,0,0,1), 714 fc = (0,0,0,0.15), 715 lw = 0., 716 ) 717 default_kwargs_Tmarkers = dict( 718 ls = 'None', 719 marker = 'o', 720 ms = 4, 721 mfc = 'w', 722 mec = 'k', 723 mew = lw, 724 ) 725 default_kwargs_Tmarker_ellipses = dict( 726 fc = 'None', 727 ec = 'k', 728 lw = lw, 729 ) 730 default_kwargs_Tmarker_labels = dict( 731 size = 8, 732 va = 'center', 733 ha = 'left', 734 linespacing = 3, 735 ) 736 737 plot_elements = {} 738 739 Ti = _np.linspace( 740 (Tmin + 273.15)**-2, 741 (Tmax + 273.15)**-2, 742 NT 743 )**-0.5 - 273.15 744 745 Tmarkers = _np.array([_ for _ in Tmarkers if _ >= Ti.min() and _ <= Ti.max()]) 746 747 if ax is None: 748 ax = _ppl.gca() 749 ax.set_xlabel(xlabel) 750 ax.set_ylabel(ylabel) 751 752 Xe = self.D47_calib_function(Ti) 753 Ye = self.D48_calib_function(Ti) 754 755 if show_eqline: 756 plot_elements['eqline'], = ax.plot( 757 _unp.nominal_values(Xe), 758 _unp.nominal_values(Ye), 759 **(default_kwargs_eqline | kwargs_eqline), 760 ) 761 762 if show_confidence: 763 plot_elements['confidence'] = self.plot_D95_confidence_band( 764 p = confidence_pvalue, 765 ax = ax, 766 **(default_kwargs_confidence | kwargs_confidence), 767 ) 768 769 Xm = self.D47_calib_function(Tmarkers) 770 Ym = self.D48_calib_function(Tmarkers) 771 if Tmarkers.size > 0: 772 plot_elements['Tm'] = ax.plot( 773 _unp.nominal_values(Xm), 774 _unp.nominal_values(Ym), 775 **(default_kwargs_Tmarkers | kwargs_Tmarkers), 776 ) 777 if show_Tmarker_ellipses: 778 plot_elements['Tme'] = conf_ellipse( 779 Xm, 780 Ym, 781 ax = ax, 782 **(default_kwargs_Tmarker_ellipses | kwargs_Tmarker_ellipses), 783 ) 784 if show_Tmarker_labels: 785 plot_elements['Tml'] = [] 786 for x,y,t in zip(Xm, Ym, Tmarkers): 787 plot_elements['Tml'].append( 788 ax.text( 789 x.n, 790 y.n, 791 f'\n${t:.0f}\\,$°C', 792 **(default_kwargs_Tmarker_labels | kwargs_Tmarker_labels), 793 ) 794 ) 795 796 ax.autoscale_view() 797 798 data = dict( 799 Te = Ti, 800 D47e = Xe, 801 D48e = Ye, 802 Tm = Tmarkers, 803 D47m = Xm, 804 D48m = Ym, 805 ) 806 807 return data, plot_elements 808 809 def _compute_p_and_D48eq_from_D47eq( 810 self, 811 D47, 812 D48, 813 D47eq, 814 ignore_calib_uncertainties = False, 815 ): 816 """ 817 Used by the various `Engine.nearest_D47eq()` methods 818 """ 819 N = D47.size 820 821 # Compute fit residuals for p values 822 if ignore_calib_uncertainties: 823 R = _cd.uarray(_np.concatenate(( 824 D47 - self.D47u_as_function_of_D47n(D47eq.n).n, 825 D48 - self.D48u_as_function_of_D47n(D47eq.n).n, 826 ))) 827 else: 828 R = _cd.uarray(_np.concatenate(( 829 D47 - self.D47u_as_function_of_D47n(D47eq.n), 830 D48 - self.D48u_as_function_of_D47n(D47eq.n), 831 ))) 832 833 # Compute p values 834 p = _np.zeros((N,)) 835 for k in range(N): 836 r = R[k::N] 837 z2 = r.m 838 p[k] = 1-_chi2.cdf(z2, 1) 839 840 # Compute D48eq 841 D48eq = self.D48u_as_function_of_D47n(D47eq) 842 843 return p, D48eq 844 845 def nearest_D47eq( 846 self, 847 D47: _cd.uarray, 848 D48: _cd.uarray, 849 ignore_calib_uncertainties: bool = False, 850 ): 851 """ 852 Computes a `correldata.uarray` of *equilibrium* Δ<sub>47</sub> values, each of which is 853 the closest (in the OGLS sense) to one (Δ<sub>47</sub>, Δ<sub>48</sub>) observation 854 considered independently of the others. 855 856 Also returns an array of corresponding p-values taking into account errors in Δ<sub>47</sub> 857 and Δ<sub>48</sub> (and any covariance between the two) as well as errors in the 858 Δ<sub>47</sub> and Δ<sub>48</sub> calibrations. 859 860 > [!NOTE] 861 > This is both the fastest and the strongly recommended version of this calculation. 862 > It is expected to yield an `uarray` with reasonably accurate covariance between the 863 > `D47eq` values, but also between `D47eq` and all other variables. 864 """ 865 866 N = D47.size 867 N47 = self.D47_coefs.size 868 N48 = self.D48_coefs.size 869 D47eq = D47 * 0 870 871 # _np.set_printoptions(threshold = _np.inf) 872 # _np.set_printoptions(linewidth = _np.inf) 873 874 for i in range(N): 875 def fun(*args): # args = (D47, D48, *D47_calib_coefs, *D48_calib_coefs) 876 877 args = _np.array(args) 878 D47_n = args[0] 879 D48_n = args[1] 880 D47_calib_coefs_n = args[-N48-N47:-N48] 881 D48_calib_coefs_n = args[-N48:] 882 883 params = _lmfit.Parameters() 884 params.add('D47eq', value = D47_n) 885 886 D47_u = _cd.uarray([_uc.ufloat(D47_n, D47.s[i])]) 887 D48_u = _cd.uarray([_uc.ufloat(D48_n, D48.s[i])]) 888 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 889 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 890 891 D47i = D4x_calib_function( 892 self.interp.T, 893 D47_calib_coefs_u, 894 return_without_uncertainties = False, 895 ignore_calib_uncertainties = ignore_calib_uncertainties, 896 ) 897 898 D48i = D4x_calib_function( 899 self.interp.T, 900 D48_calib_coefs_u, 901 return_without_uncertainties = False, 902 ignore_calib_uncertainties = ignore_calib_uncertainties, 903 ) 904 905 D47_interp = uarray_compatible_interp(D47i.n, D47i) 906 D48_interp = uarray_compatible_interp(D47i.n, D48i) 907 908 def cost_fun(p): 909 R = _cd.uarray(_np.concatenate(( 910 D47_u - D47_interp(p['D47eq'].value), 911 D48_u - D48_interp(p['D47eq'].value), 912 ))) 913 914 invS = _np.linalg.inv(R.covar) 915 L = _cholesky(invS) 916 917 return L @ R.n 918 919 minresult = _lmfit.minimize( 920 cost_fun, 921 params, 922 method = 'least_squares', 923 scale_covar = False, 924 jac = '3-point', 925 ) 926 # slower but yields very similar results: 927 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 928 929 return minresult.params['D47eq'].value 930 931 wrapped_fun = _uc.wrap(fun) 932 D47eq[i] = wrapped_fun(D47[i], D48[i], *self.D47_coefs, *self.D48_coefs) 933 934 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 935 936 return D47eq, D48eq, p 937 938 def joint_nearest_D47eq( 939 self, 940 D47: _cd.uarray, 941 D48: _cd.uarray, 942 ignore_calib_uncertainties: bool = False, 943 ): 944 """ 945 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 946 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 947 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 948 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 949 Δ<sub>48</sub> calibrations. 950 951 > [!CAUTION] 952 > Caution: the use of this function is **not generally recommended** except for 953 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 954 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 955 > as done by the recommended function `nearest_D47eq()`. 956 957 This is the most complete but slowest and not recommended version of this calculation. 958 It is expected to yield an `uarray` with reasonably accurate covariance between the 959 `D47eq` values, but also between `D47eq` and all other variables. 960 961 A faster but incomplete and potentially less accurate version of this calculation is 962 provided by `lazy_joint_nearest_D47eq()`. 963 """ 964 965 N = D47.size 966 N47 = self.D47_coefs.size 967 N48 = self.D48_coefs.size 968 969 def fun(j, *args): 970 971 args = _np.array(args) 972 D47_n = args[:N] 973 D48_n = args[N:2*N] 974 D47_calib_coefs_n = args[-N48-N47:-N48] 975 D48_calib_coefs_n = args[-N48:] 976 977 params = _lmfit.Parameters() 978 for k in range(N): 979 params.add(f'D47eq{k}', value = D47_n[k]) 980 981 D47_u = _cd.uarray(_uc.correlated_values(D47_n, D47.covar)) 982 D48_u = _cd.uarray(_uc.correlated_values(D48_n, D48.covar)) 983 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 984 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 985 986 D47i = D4x_calib_function( 987 self.interp.T, 988 D47_calib_coefs_u, 989 return_without_uncertainties = False, 990 ignore_calib_uncertainties = ignore_calib_uncertainties, 991 ) 992 993 D48i = D4x_calib_function( 994 self.interp.T, 995 D48_calib_coefs_u, 996 return_without_uncertainties = False, 997 ignore_calib_uncertainties = ignore_calib_uncertainties, 998 ) 999 1000 D47_interp = uarray_compatible_interp(D47i.n, D47i) 1001 D48_interp = uarray_compatible_interp(D47i.n, D48i) 1002 1003 def cost_fun(p): 1004 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1005 R = _cd.uarray(_np.concatenate(( 1006 D47_u - D47_interp(_D47eq), 1007 D48_u - D48_interp(_D47eq), 1008 ))) 1009 1010 invS = _np.linalg.inv(R.covar) 1011 L = _cholesky(invS) 1012 1013 # print(((L @ R.n)**2).sum()) 1014 return L @ R.n 1015 1016 minresult = _lmfit.minimize( 1017 cost_fun, 1018 params, 1019 method = 'least_squares', 1020 scale_covar = False, 1021 jac = '3-point', 1022 ) 1023 # slower but yields very similar results: 1024 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 1025 1026 return minresult.params[f'D47eq{j}'].value 1027 1028 wrapped_fun = _uc.wrap(fun) 1029 1030 D47eq = _cd.uarray([wrapped_fun(j, *D47, *D48, *self.D47_coefs, *self.D48_coefs) for j in range(N)]) 1031 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1032 1033 return D47eq, D48eq, p 1034 1035 def lazy_joint_nearest_D47eq( 1036 self, 1037 D47: _cd.uarray, 1038 D48: _cd.uarray, 1039 ignore_calib_uncertainties: bool = False, 1040 ): 1041 """ 1042 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 1043 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 1044 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 1045 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 1046 Δ<sub>48</sub> calibrations. 1047 1048 > [!CAUTION] 1049 > Caution: the use of this function is **not generally recommended** except for 1050 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 1051 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 1052 > as done by the recommended function `nearest_D47eq()`. 1053 1054 This is a faster but incomplete version of this calculation. It is expected to yield an 1055 `uarray` with roughly accurate covariance between the `Teq` values, but without computing 1056 the covariance with any other variables. 1057 1058 A slower but complete and more accurate version of this calculation is provided by 1059 `joint_nearest_D47eq()`. 1060 """ 1061 1062 N = D47.size 1063 1064 params = _lmfit.Parameters() 1065 for k in range(N): 1066 params.add(f'D47eq{k}', value = D47[k].n) 1067 1068 def cost_fun(p, ignore_calib_uncertainties = ignore_calib_uncertainties): 1069 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1070 1071 if ignore_calib_uncertainties: 1072 R = _cd.uarray(_np.concatenate(( 1073 D47 - self.D47u_as_function_of_D47n(_D47eq).n, 1074 D48 - self.D48u_as_function_of_D47n(_D47eq).n, 1075 ))) 1076 else: 1077 R = _cd.uarray(_np.concatenate(( 1078 D47 - self.D47u_as_function_of_D47n(_D47eq), 1079 D48 - self.D48u_as_function_of_D47n(_D47eq), 1080 ))) 1081 1082 invS = _np.linalg.inv(R.covar) 1083 L = _cholesky(invS) 1084 1085 # print(((L @ R.n)**2).sum()) 1086 return L @ R.n 1087 1088 minresult = _lmfit.minimize( 1089 cost_fun, 1090 params, 1091 method = 'least_squares', 1092 scale_covar = False, 1093 jac = '3-point', 1094 ) 1095 1096 D47eq = _cd.uarray([minresult.uvars[f'D47eq{k}'] for k in range(N)]) 1097 1098 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1099 1100 return D47eq, D48eq, p 1101 1102 def projected_D47eq( 1103 self, 1104 D47: _cd.uarray, 1105 D48: _cd.uarray, 1106 kinetic_slope: (float | _uc.UFloat), 1107 ): 1108 """ 1109 Projects one or more (Δ<sub>47</sub>, Δ<sub>48</sub>) observations onto the equlibrium curve 1110 following a kinetic fractionation vector with a given slope (∂Δ<sub>48</sub>/∂Δ<sub>47</sub>). 1111 1112 **Arguments** 1113 * `D47`: observed Δ<sub>47</sub> value(s) 1114 * `D48`: observed Δ<sub>48</sub> value(s) 1115 * `kinetic_slope`: kinetic fractionation slopw, with or without uncertainty 1116 1117 Returns a tuple of uarrays corresponding to the projected Δ<sub>47</sub> and Δ<sub>48</sub> values. 1118 1119 > [!NOTE] 1120 > This is not a least-squares minimization problem but a direct calculation, and should thus 1121 > be much faster than the various `CorelData.nearestD47eq()` methods. 1122 """ 1123 1124 D47 = _cd.uarray(D47) 1125 D48 = _cd.uarray(D48) 1126 N = D47.size 1127 N47c = self.D47_coefs.size 1128 N48c = self.D48_coefs.size 1129 D47p = D47 * 0 1130 1131 for i in range(N): 1132 1133 # function to solve 1134 def fun(x, *args): # args = (D47, D48, kinetic_slope, *self.D47_coefs, *self.D48_coefs) 1135 1136 args = _np.array(args) 1137 D47_n = args[0] 1138 D48_n = args[1] 1139 kslope_n = args[2] 1140 D47_calib_coefs_n = args[-N48c-N47c:-N48c] 1141 D48_calib_coefs_n = args[-N48c:] 1142 1143 D47i = D4x_calib_function( 1144 self.interp.T, 1145 D47_calib_coefs_n, 1146 return_without_uncertainties = False, 1147 ) 1148 1149 D48i = D4x_calib_function( 1150 self.interp.T, 1151 D48_calib_coefs_n, 1152 return_without_uncertainties = False, 1153 ) 1154 1155 D48_interp = uarray_compatible_interp(D47i, D48i) 1156 1157 return D48_n - D48_interp(x) - kslope_n * (D47_n - x) 1158 1159 def g(*args): 1160 return _fsolve(fun, [100.], args = args)[0] 1161 1162 wg = _uc.wrap(g) 1163 1164 D47p[i] = wg( 1165 D47[i], 1166 D48[i], 1167 kinetic_slope, 1168 *self.D47_coefs, 1169 *self.D48_coefs, 1170 ) 1171 1172 _, D48p = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47p, ignore_calib_uncertainties = False) 1173 1174 return D47p, D48p 1175 1176 def Teq_pdf( 1177 self, 1178 D47: _uc.ufloat, 1179 Tmin: (float | None) = None, 1180 Tmax: (float | None) = None, 1181 Tinc: float = 0.2, 1182 default_D47_sigmas: float = 4.0, 1183 ignore_calib_uncertainties: bool = False, 1184 run_qmc: bool = False, 1185 N_qmc: int = 1024, 1186 ): 1187 """ 1188 Compute the unit-normalized probability distribution function (PDF) of the 1189 equilibrium temperature (`Teq`) for a given (`UFloat`) value of Δ<sub>47</sub>. 1190 1191 **Arguments** 1192 * `D47`: Δ<sub>47</sub> value (with uncertainty) 1193 * `Tmin`: minimum temperature over which to compute the PDF; if not specified, 1194 use temperature corresponding to `D47.n + `default_D47_sigmas` * D47.s` 1195 * `Tmax`: maximum temperature over which to compute the PDF; if not specified, 1196 use temperature corresponding to `D47.n - `default_D47_sigmas` * D47.s` 1197 * `Tinc`: temperature increment over which to compute the PDF 1198 * `default_D47_sigmas`: see `Tmin` and `Tmin` above 1199 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 1200 * `run_qmc`: whether to also run a Quasi Monte carlo simulation to estimate the PDF 1201 * `N_qmc`: number of iterations in the above Quasi Monte Carlo simulation 1202 1203 **Returns** 1204 * `Ti`: Evenly-spaced array of temperature values over which the PDF is computed 1205 * `pdf`: PDF evaluated over `Ti` 1206 * `Tqmc` (only returned if `run_qmc = True`): array of `N_qmc` temperature values 1207 computed in the Quasi Monte Carlo simulation 1208 """ 1209 1210 if Tmin is None: 1211 Tmin = _np.floor(self.T_as_function_of_D47( 1212 D47.n + default_D47_sigmas * D47.s, 1213 ignore_calib_uncertainties = ignore_calib_uncertainties, 1214 ).n) 1215 1216 if Tmax is None: 1217 Tmax = _np.ceil(self.T_as_function_of_D47( 1218 D47.n - default_D47_sigmas * D47.s, 1219 ignore_calib_uncertainties = ignore_calib_uncertainties, 1220 ).n) 1221 1222 assert Tmin < Tmax, "Tmax must be strictly greater than Tmin" 1223 assert Tinc > 0, "Tinc must be strictly greater than zero" 1224 1225 # compute interpolated Ti values 1226 Ti = _np.arange(Tmin, Tmax+Tinc, Tinc) 1227 1228 pdf = transform_pdf_monotonic( 1229 f_inv = lambda T: D4x_calib_function( 1230 T, 1231 self.D47_coefs, 1232 return_without_uncertainties = ignore_calib_uncertainties, 1233 ignore_calib_uncertainties = ignore_calib_uncertainties, 1234 ), 1235 df_inv = lambda T: D4x_calib_derivative( 1236 T, 1237 self.D47_coefs, 1238 return_without_uncertainties = ignore_calib_uncertainties, 1239 ignore_calib_uncertainties = ignore_calib_uncertainties, 1240 ), 1241 mu_x = D47.n, 1242 sigma_x = D47.s, 1243 yi = Ti, 1244 ) 1245 1246 if run_qmc: 1247 1248 from scipy.stats import qmc 1249 from tqdm.rich import tqdm 1250 1251 #parameters to jiggle 1252 input_params = _cd.uarray([D47, *self.D47_coefs]) 1253 1254 # QMC sampler for the correlation matrix of these parameters 1255 qmc_dist = qmc.MultivariateNormalQMC( 1256 mean = input_params.n*0, 1257 cov = input_params.cor, 1258 ) 1259 1260 # QMC samples 1261 qmc_draws = input_params.n + qmc_dist.random(N_qmc) * input_params.s 1262 1263 # initialize T_qmc 1264 Tqmc = _cd.uarray(_np.zeros((N_qmc,))) 1265 1266 for k in tqdm(range(N_qmc)): 1267 # jiggled D47 and D47coefs 1268 _D47 = qmc_draws[k,0] 1269 if ignore_calib_uncertainties: 1270 _coefs = self.D47_coefs 1271 else: 1272 _coefs = _cd.uarray(_uc.correlated_values(qmc_draws[k,1:], self.D47_coefs.covar)) 1273 1274 # jiggled D47 1275 _D47i = D4x_calib_function(self.interp.T, _coefs) 1276 _f = uarray_compatible_interp(_D47i.n, self.interp.T) 1277 Tqmc[k] = _f(_D47) 1278 1279 return Ti, pdf, Tqmc 1280 1281 return Ti, pdf
Underlying engine to compute and plot nearest equilibrium temperatures and projected temperatures based on a consistent pair of Δ47, Δ48 calibrations.
Engine( D47_coefs: correldata.uarray | ArrayLike | None = None, D48_coefs: correldata.uarray | ArrayLike | None = None, Tmin_interp: float = -23.0, Tmax_interp: float = 1277.0, N_interp: float = 201)
416 def __init__( 417 self, 418 D47_coefs: (_cd.uarray | ArrayLike | None) = None, 419 D48_coefs: (_cd.uarray | ArrayLike | None) = None, 420 Tmin_interp: float = -23.0, 421 Tmax_interp: float = 1277.0, 422 N_interp: float = 201, 423 ): 424 """ 425 **Arguments** 426 * `D47_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 427 * `D48_coefs`: `ndarray` or `uarray` of coefficients to use instead of default ones, ordered as (a0, a1, a2...) 428 * `Tmin_interp`: minimum temperature over which to interpolate for inverse function computations 429 * `Tmax_interp`: maximum temperature over which to interpolate for inverse function computations 430 * `N_interp`: number of points (equally-spaced in 1/T space) over which to interpolate for inverse function computations 431 """ 432 433 self.D47_coefs = Engine.D47_calib_coefs if D47_coefs is None else D47_coefs 434 """The Δ<sub>47</sub> calibration coefficients used by this `Engine` instance""" 435 436 self.D48_coefs = Engine.D48_calib_coefs if D48_coefs is None else D48_coefs 437 """The Δ<sub>48</sub> calibration coefficients used by this `Engine` instance""" 438 439 self.interp = _Interpolation() 440 """ 441 Holds equilibrium Δ<sub>47</sub> and Δ<sub>48</sub> values (ufloats) interpolated 442 along an array of T values (regularly spaced increments of 1/T<sup>2</sup>). 443 444 * `interp.T`: interpolation T values (floats) in regularly spaced increments of 1/T<sup>2</sup> 445 * `interp.D47`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T` 446 * `interp.D48`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T` 447 * `interp.D47_no_calib_errors`: Equilibrium Δ<sub>47</sub> values (ufloats) interpolated along `interp.T`, 448 ignoring calibration uncertainties 449 * `interp.D48_no_calib_errors`: Equilibrium Δ<sub>48</sub> values (ufloats) interpolated along `interp.T`, 450 ignoring calibration uncertainties 451 """ 452 453 self.interp.T = _np.linspace( 454 (Tmax_interp+273.15)**-2, 455 (Tmin_interp+273.15)**-2, 456 N_interp, 457 )**-0.5 - 273.15 458 459 self.interp.D47 = self.D47_calib_function( 460 self.interp.T, 461 return_without_uncertainties = False, 462 ignore_calib_uncertainties = False, 463 ) 464 465 self.interp.D47_no_calib_errors = self.D47_calib_function( 466 self.interp.T, 467 return_without_uncertainties = False, 468 ignore_calib_uncertainties = True, 469 ) 470 471 self.interp.D48 = self.D48_calib_function( 472 self.interp.T, 473 return_without_uncertainties = False, 474 ignore_calib_uncertainties = False, 475 ) 476 477 self.interp.D48_no_calib_errors = self.D48_calib_function( 478 self.interp.T, 479 return_without_uncertainties = False, 480 ignore_calib_uncertainties = True, 481 ) 482 483 self.interp.D47u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D47) 484 self.interp.D48u_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.D48) 485 486 #inverse D47 calibration (ignoring calibration errors) 487 self.interp.Teq_as_function_of_D47n = uarray_compatible_interp(self.interp.D47.n, self.interp.T) 488 #inverse D47 calibration (including calibration errors) 489 self.interp.Teq_as_function_of_D47u = uarray_compatible_interp(self.interp.D47, self.interp.T)
Arguments
D47_coefs:ndarrayoruarrayof coefficients to use instead of default ones, ordered as (a0, a1, a2...)D48_coefs:ndarrayoruarrayof coefficients to use instead of default ones, ordered as (a0, a1, a2...)Tmin_interp: minimum temperature over which to interpolate for inverse function computationsTmax_interp: maximum temperature over which to interpolate for inverse function computationsN_interp: number of points (equally-spaced in 1/T space) over which to interpolate for inverse function computations
D47_calib_coefs =
uarray([0.17437754366432887+/-0.004911105567257294,
-18.14215245127414+/-5.632326472234855,
42657.22989162373+/-1277.12751715908], dtype=object)
Default (OGLS23) Δ47 calibration coefficients based on Daëron & Vermeesch (2024)
D48_calib_coefs =
uarray([0.121349237888+/-0.0039004854072400012,
6.22931985613+/-0.3289676145899999, -13481.983494+/-711.977559735,
9336714.66607+/-493067.75422400003,
-770413883.573+/-40685214.9801], dtype=object)
Default Δ48 calibration coefficients based on Fiebig et al. (2024)
interp
Holds equilibrium Δ47 and Δ48 values (ufloats) interpolated along an array of T values (regularly spaced increments of 1/T2).
interp.T: interpolation T values (floats) in regularly spaced increments of 1/T2interp.D47: Equilibrium Δ47 values (ufloats) interpolated alonginterp.Tinterp.D48: Equilibrium Δ48 values (ufloats) interpolated alonginterp.Tinterp.D47_no_calib_errors: Equilibrium Δ47 values (ufloats) interpolated alonginterp.T, ignoring calibration uncertaintiesinterp.D48_no_calib_errors: Equilibrium Δ48 values (ufloats) interpolated alonginterp.T, ignoring calibration uncertainties
def
T_as_function_of_D47( self, D47: correldata.uarray | ArrayLike, ignore_calib_uncertainties: bool = False):
491 def T_as_function_of_D47( 492 self, 493 D47: (_cd.uarray | ArrayLike), 494 ignore_calib_uncertainties: bool = False, 495 ): 496 """ 497 Provided with one or more Δ<sub>47</sub> values (floats or ufloats), return ufloats for the 498 corresponding equilibrium T values (ufloats with or without Δ<sub>47</sub> calibration uncertainties). 499 500 **Arguments** 501 * `D47`: array of Δ<sub>47</sub> values 502 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 503 """ 504 if ignore_calib_uncertainties: 505 return _cd.uarray(self.interp.Teq_as_function_of_D47n(D47)) 506 else: 507 return _cd.uarray(self.interp.Teq_as_function_of_D47u(D47))
Provided with one or more Δ47 values (floats or ufloats), return ufloats for the corresponding equilibrium T values (ufloats with or without Δ47 calibration uncertainties).
Arguments
D47: array of Δ47 valuesignore_calib_uncertainties: whether to propagate calibration uncertainties
def
D47u_as_function_of_D47n(self, D47: ArrayLike):
509 def D47u_as_function_of_D47n( 510 self, 511 D47: ArrayLike 512 ): 513 """ 514 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 515 equilibrium Δ<sub>47</sub> values (ufloats with Δ<sub>47</sub> calibration uncertainties). 516 """ 517 return _cd.uarray(self.interp.D47u_as_function_of_D47n(D47))
Provided with one or more Δ47 values (floats), return ufloats for the corresponding equilibrium Δ47 values (ufloats with Δ47 calibration uncertainties).
def
D48u_as_function_of_D47n(self, D47: ArrayLike):
519 def D48u_as_function_of_D47n( 520 self, 521 D47: ArrayLike 522 ): 523 """ 524 Provided with one or more Δ<sub>47</sub> values (floats), return ufloats for the corresponding 525 equilibrium Δ<sub>48</sub> values (ufloats with Δ<sub>48</sub> calibration uncertainties). 526 """ 527 return _cd.uarray(self.interp.D48u_as_function_of_D47n(D47))
Provided with one or more Δ47 values (floats), return ufloats for the corresponding equilibrium Δ48 values (ufloats with Δ48 calibration uncertainties).
def
D47_calib_function( self, T: float | uncertainties.core.AffineScalarFunc | correldata.uarray, return_without_uncertainties: bool = False, ignore_calib_uncertainties: bool = False):
529 def D47_calib_function( 530 self, 531 T: (float | _uc.UFloat | _cd.uarray), 532 return_without_uncertainties: bool = False, 533 ignore_calib_uncertainties: bool = False, 534 ): 535 return D4x_calib_function( 536 T = T, 537 coefs = self.D47_coefs, 538 return_without_uncertainties = return_without_uncertainties, 539 ignore_calib_uncertainties = ignore_calib_uncertainties, 540 )
Arguments
T: temperature(s) for which to compute Δ47return_without_uncertainties: ifTrue, returns Δ47 values without error propagation of any kindignore_calib_uncertainties: whether to propagate calibration uncertainties
Returns equilibrium Δ47 value(s) corresponding to T value(s)
def
D48_calib_function( self, T: float | uncertainties.core.AffineScalarFunc | correldata.uarray, return_without_uncertainties: bool = False, ignore_calib_uncertainties: bool = False):
542 def D48_calib_function( 543 self, 544 T: (float | _uc.UFloat | _cd.uarray), 545 return_without_uncertainties: bool = False, 546 ignore_calib_uncertainties: bool = False, 547 ): 548 return D4x_calib_function( 549 T = T, 550 coefs = self.D48_coefs, 551 return_without_uncertainties = return_without_uncertainties, 552 ignore_calib_uncertainties = ignore_calib_uncertainties, 553 )
Arguments
T: temperature(s) for which to compute Δ48return_without_uncertainties: ifTrue, returns Δ48 values without error propagation of any kindignore_calib_uncertainties: whether to propagate calibration uncertainties
Returns equilibrium Δ48 value(s) corresponding to T value(s)
def
T_ellipse( self, T: numpy.ndarray | correldata.uarray, p: float = 0.95, CM: numpy.ndarray | None = None, Tse: numpy.ndarray | float | None = None, plot: bool = True, ax: matplotlib.axes._axes.Axes | None = None, **kwargs) -> list:
558 def T_ellipse( 559 self, 560 T: (_np.ndarray | _cd.uarray), 561 p: float = 0.95, 562 CM: (_np.ndarray | None) = None, 563 Tse: (_np.ndarray | float | None) = None, 564 plot: bool = True, 565 ax: (_ppl.Axes | None) = None, 566 **kwargs, 567 ) -> list: 568 """ 569 Plot the joint `p`-level confidence ellipses in (Δ<sub>47</sub>, Δ<sub>48</sub>) 570 space, for temperatures equal to the elements of `T`, and return a list of the 571 `Ellipse` objects thus created. 572 573 **Arguments** 574 * `T`: `ndarray` or `uarray` of temperatures to plot 575 * `p`: confidence level 576 * `plot`: whether to plot the ellipse or not. If `False`, return a list of 577 `(x_center, y_center, width, height, angle)` elements 578 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 579 * `kwargs`: passed to `matplotlib.patches.Ellipse()` 580 """ 581 _T = _cd.as_uarray(T, CM = CM, Xse = Tse) 582 return conf_ellipse( 583 self.D47_calib_function(_T), 584 self.D48_calib_function(_T), 585 p = p, 586 plot = plot, 587 ax = ax, 588 **kwargs, 589 )
Plot the joint p-level confidence ellipses in (Δ47, Δ48)
space, for temperatures equal to the elements of T, and return a list of the
Ellipse objects thus created.
Arguments
T:ndarrayoruarrayof temperatures to plotp: confidence levelplot: whether to plot the ellipse or not. IfFalse, return a list of(x_center, y_center, width, height, angle)elementsax: which instance ofmatplotlib.axes.Axesto draw in; use current axes ifax=None.kwargs: passed tomatplotlib.patches.Ellipse()
def
plot_D95_confidence_band( self, p: float = 0.95, Ti: ArrayLike | None = None, plot: bool = True, ax: matplotlib.axes._axes.Axes | None = None, **kwargs):
591 def plot_D95_confidence_band( 592 self, 593 p: float = 0.95, 594 Ti: (ArrayLike | None) = None, 595 plot: bool = True, 596 ax: (_ppl.Axes | None) = None, 597 **kwargs, 598 ): 599 """ 600 Plot, for a given p-value, the confidence band of the thermodynamic equilibrium curve 601 in (Δ<sub>47</sub>, Δ<sub>48</sub>) space. 602 603 **Arguments** 604 * `p`: confidence level 605 * `Ti`: array of temperatures over which to evaluate confidence band (default: use `interp.T` attribute instead) 606 * `plot`: whether to plot the confidence band or not. If `False`, return the (N,2) array of polygon nodes 607 * `ax`: `Axes` instance to plot to (default: use current Axes) 608 * `kwargs`: passed to `patches.Polygon()` 609 610 Returns the corresponding `Polygon` instance. 611 """ 612 613 if Ti is None: 614 Ti = self.interp.T 615 616 cb = confidence_band( 617 Ti, 618 self.D47_calib_function, 619 self.D48_calib_function, 620 p, 621 ) 622 623 if plot: 624 from matplotlib import pyplot as _ppl 625 from matplotlib.patches import Polygon as _Polygon 626 627 if ax is None: 628 ax = _ppl.gca() 629 630 polygon = ax.add_patch( 631 _Polygon( 632 cb, 633 closed = True, 634 **kwargs, 635 ) 636 ) 637 return polygon 638 else: 639 return cb
Plot, for a given p-value, the confidence band of the thermodynamic equilibrium curve in (Δ47, Δ48) space.
Arguments
p: confidence levelTi: array of temperatures over which to evaluate confidence band (default: useinterp.Tattribute instead)plot: whether to plot the confidence band or not. IfFalse, return the (N,2) array of polygon nodesax:Axesinstance to plot to (default: use current Axes)kwargs: passed topatches.Polygon()
Returns the corresponding Polygon instance.
def
plot_D95_equilibrium( self, Tmin: float = 0.0, Tmax: float = 1000.0, NT: int = 101, Tmarkers: ArrayLike = [0, 25, 100, 250, 1000], kwargs_Tmarkers: dict = {}, show_Tmarker_labels: bool = True, kwargs_Tmarker_labels: dict = {}, show_Tmarker_ellipses: bool = False, kwargs_Tmarker_ellipses: dict = {}, show_eqline: bool = True, kwargs_eqline: dict = {}, show_confidence: bool = True, confidence_pvalue: float = 0.95, kwargs_confidence: dict = {}, ax: matplotlib.axes._axes.Axes | None = None, xlabel: str = '$Δ_{47}$ [‰]', ylabel: str = '$Δ_{48}$ [‰]', lw: float = 0.7) -> (<class 'dict'>, <class 'dict'>):
642 def plot_D95_equilibrium( 643 self, 644 Tmin: float = 0., 645 Tmax: float = 1000., 646 NT: int = 101, 647 Tmarkers: _np.typing.ArrayLike = [0, 25, 100, 250, 1000], 648 kwargs_Tmarkers: dict = {}, 649 show_Tmarker_labels: bool = True, 650 kwargs_Tmarker_labels: dict = {}, 651 show_Tmarker_ellipses: bool = False, 652 kwargs_Tmarker_ellipses: dict = {}, 653 show_eqline: bool = True, 654 kwargs_eqline: dict = {}, 655 show_confidence: bool = True, 656 confidence_pvalue: float = 0.95, 657 kwargs_confidence: dict = {}, 658 ax: (_ppl.Axes | None) = None, 659 xlabel: str = '$Δ_{47}$ [‰]', 660 ylabel: str = '$Δ_{48}$ [‰]', 661 lw: float = 0.7, 662 ) -> (dict, dict): 663 """ 664 Plot a thermodynamic equilibrium curve in (Δ<sub>47</sub>, Δ<sub>48</sub>) space 665 as a function of temperature. 666 667 **Arguments** 668 * `Tmin`: minimum T to plot 669 * `Tmax`: maximum T to plot 670 * `NT`: number of steps in equilibrium curve (interpolated at constant steps in 1/T<sup>2</sup> space) 671 * `Tmarkers`: T markers to add along the curve 672 * `kwargs_Tmarkers`: passed to `plot()` when plotting T markers 673 * `show_Tmarker_labels`: whether to add T labels to T markers 674 * `kwargs_Tmarker_labels`: passed to `text()` when plotting T markers 675 * `show_Tmarker_ellipses`: whether to add confidence ellipses to T markers 676 * `kwargs_Tmarker_ellipses`: passed to `T_ellipses()` when plotting T marker ellipses 677 * `show_eqline`: whether to plot the equilibrium curve itself 678 * `kwargs_eqline`: passed to `plot()` when plotting the equilibrium curve 679 * `show_confidence`: whether to plot the confidence band of the equilibrium curve 680 * `confidence_pvalue`: confidence level for the confidence band 681 * `kwargs_confidence`: passed to `plot_D95_confidence_band()` when plotting the confidence band 682 * `ax`: which instance of `matplotlib.axes.Axes` to draw in; use current axes if `ax` = `None`. 683 * `xlabel`: string to pass to `xlabel()` 684 * `ylabel`: string to pass to `ylabel()` 685 * `lw`: default line width for most plot elements 686 687 **Returns** 688 * `data`: a dict of the T, Δ<sub>47</sub> and Δ<sub>48</sub> values generated for this plot: 689 - `Te` : temperature interpolated along the equilibrium curve 690 - `D47e`: Δ<sub>47</sub> interpolated along the equilibrium curve 691 - `D48e`: Δ<sub>48</sub> interpolated along the equilibrium curve 692 - `Tm` : temperature of T markers 693 - `D47m`: Δ<sub>47</sub> of T markers 694 - `D48m`: Δ<sub>48</sub> of T markers 695 696 * `plot_elements`: a dict of the `Axes` elements generated for this plot: 697 - `eqline`: `Line2D` of the equilibrium curve 698 - `confidence`: `Polygon` object for the confidence band 699 - `Tm`: `Line2D` of the T markers 700 - `Tme`: list of `Ellipse` objects for the T marker ellipses 701 - `Tml`: list of `Text` objects for the T marker labels 702 """ 703 704 from matplotlib import pyplot as _ppl 705 706 default_kwargs_eqline = dict( 707 marker = 'None', 708 ls = '-', 709 color = 'k', 710 lw = lw, 711 ) 712 default_kwargs_confidence = dict( 713 ec = (0,0,0,1), 714 fc = (0,0,0,0.15), 715 lw = 0., 716 ) 717 default_kwargs_Tmarkers = dict( 718 ls = 'None', 719 marker = 'o', 720 ms = 4, 721 mfc = 'w', 722 mec = 'k', 723 mew = lw, 724 ) 725 default_kwargs_Tmarker_ellipses = dict( 726 fc = 'None', 727 ec = 'k', 728 lw = lw, 729 ) 730 default_kwargs_Tmarker_labels = dict( 731 size = 8, 732 va = 'center', 733 ha = 'left', 734 linespacing = 3, 735 ) 736 737 plot_elements = {} 738 739 Ti = _np.linspace( 740 (Tmin + 273.15)**-2, 741 (Tmax + 273.15)**-2, 742 NT 743 )**-0.5 - 273.15 744 745 Tmarkers = _np.array([_ for _ in Tmarkers if _ >= Ti.min() and _ <= Ti.max()]) 746 747 if ax is None: 748 ax = _ppl.gca() 749 ax.set_xlabel(xlabel) 750 ax.set_ylabel(ylabel) 751 752 Xe = self.D47_calib_function(Ti) 753 Ye = self.D48_calib_function(Ti) 754 755 if show_eqline: 756 plot_elements['eqline'], = ax.plot( 757 _unp.nominal_values(Xe), 758 _unp.nominal_values(Ye), 759 **(default_kwargs_eqline | kwargs_eqline), 760 ) 761 762 if show_confidence: 763 plot_elements['confidence'] = self.plot_D95_confidence_band( 764 p = confidence_pvalue, 765 ax = ax, 766 **(default_kwargs_confidence | kwargs_confidence), 767 ) 768 769 Xm = self.D47_calib_function(Tmarkers) 770 Ym = self.D48_calib_function(Tmarkers) 771 if Tmarkers.size > 0: 772 plot_elements['Tm'] = ax.plot( 773 _unp.nominal_values(Xm), 774 _unp.nominal_values(Ym), 775 **(default_kwargs_Tmarkers | kwargs_Tmarkers), 776 ) 777 if show_Tmarker_ellipses: 778 plot_elements['Tme'] = conf_ellipse( 779 Xm, 780 Ym, 781 ax = ax, 782 **(default_kwargs_Tmarker_ellipses | kwargs_Tmarker_ellipses), 783 ) 784 if show_Tmarker_labels: 785 plot_elements['Tml'] = [] 786 for x,y,t in zip(Xm, Ym, Tmarkers): 787 plot_elements['Tml'].append( 788 ax.text( 789 x.n, 790 y.n, 791 f'\n${t:.0f}\\,$°C', 792 **(default_kwargs_Tmarker_labels | kwargs_Tmarker_labels), 793 ) 794 ) 795 796 ax.autoscale_view() 797 798 data = dict( 799 Te = Ti, 800 D47e = Xe, 801 D48e = Ye, 802 Tm = Tmarkers, 803 D47m = Xm, 804 D48m = Ym, 805 ) 806 807 return data, plot_elements
Plot a thermodynamic equilibrium curve in (Δ47, Δ48) space as a function of temperature.
Arguments
Tmin: minimum T to plotTmax: maximum T to plotNT: number of steps in equilibrium curve (interpolated at constant steps in 1/T2 space)Tmarkers: T markers to add along the curvekwargs_Tmarkers: passed toplot()when plotting T markersshow_Tmarker_labels: whether to add T labels to T markerskwargs_Tmarker_labels: passed totext()when plotting T markersshow_Tmarker_ellipses: whether to add confidence ellipses to T markerskwargs_Tmarker_ellipses: passed toT_ellipses()when plotting T marker ellipsesshow_eqline: whether to plot the equilibrium curve itselfkwargs_eqline: passed toplot()when plotting the equilibrium curveshow_confidence: whether to plot the confidence band of the equilibrium curveconfidence_pvalue: confidence level for the confidence bandkwargs_confidence: passed toplot_D95_confidence_band()when plotting the confidence bandax: which instance ofmatplotlib.axes.Axesto draw in; use current axes ifax=None.xlabel: string to pass toxlabel()ylabel: string to pass toylabel()lw: default line width for most plot elements
Returns
data: a dict of the T, Δ47 and Δ48 values generated for this plot:Te: temperature interpolated along the equilibrium curveD47e: Δ47 interpolated along the equilibrium curveD48e: Δ48 interpolated along the equilibrium curveTm: temperature of T markersD47m: Δ47 of T markersD48m: Δ48 of T markers
plot_elements: a dict of theAxeselements generated for this plot:eqline:Line2Dof the equilibrium curveconfidence:Polygonobject for the confidence bandTm:Line2Dof the T markersTme: list ofEllipseobjects for the T marker ellipsesTml: list ofTextobjects for the T marker labels
def
nearest_D47eq( self, D47: correldata.uarray, D48: correldata.uarray, ignore_calib_uncertainties: bool = False):
845 def nearest_D47eq( 846 self, 847 D47: _cd.uarray, 848 D48: _cd.uarray, 849 ignore_calib_uncertainties: bool = False, 850 ): 851 """ 852 Computes a `correldata.uarray` of *equilibrium* Δ<sub>47</sub> values, each of which is 853 the closest (in the OGLS sense) to one (Δ<sub>47</sub>, Δ<sub>48</sub>) observation 854 considered independently of the others. 855 856 Also returns an array of corresponding p-values taking into account errors in Δ<sub>47</sub> 857 and Δ<sub>48</sub> (and any covariance between the two) as well as errors in the 858 Δ<sub>47</sub> and Δ<sub>48</sub> calibrations. 859 860 > [!NOTE] 861 > This is both the fastest and the strongly recommended version of this calculation. 862 > It is expected to yield an `uarray` with reasonably accurate covariance between the 863 > `D47eq` values, but also between `D47eq` and all other variables. 864 """ 865 866 N = D47.size 867 N47 = self.D47_coefs.size 868 N48 = self.D48_coefs.size 869 D47eq = D47 * 0 870 871 # _np.set_printoptions(threshold = _np.inf) 872 # _np.set_printoptions(linewidth = _np.inf) 873 874 for i in range(N): 875 def fun(*args): # args = (D47, D48, *D47_calib_coefs, *D48_calib_coefs) 876 877 args = _np.array(args) 878 D47_n = args[0] 879 D48_n = args[1] 880 D47_calib_coefs_n = args[-N48-N47:-N48] 881 D48_calib_coefs_n = args[-N48:] 882 883 params = _lmfit.Parameters() 884 params.add('D47eq', value = D47_n) 885 886 D47_u = _cd.uarray([_uc.ufloat(D47_n, D47.s[i])]) 887 D48_u = _cd.uarray([_uc.ufloat(D48_n, D48.s[i])]) 888 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 889 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 890 891 D47i = D4x_calib_function( 892 self.interp.T, 893 D47_calib_coefs_u, 894 return_without_uncertainties = False, 895 ignore_calib_uncertainties = ignore_calib_uncertainties, 896 ) 897 898 D48i = D4x_calib_function( 899 self.interp.T, 900 D48_calib_coefs_u, 901 return_without_uncertainties = False, 902 ignore_calib_uncertainties = ignore_calib_uncertainties, 903 ) 904 905 D47_interp = uarray_compatible_interp(D47i.n, D47i) 906 D48_interp = uarray_compatible_interp(D47i.n, D48i) 907 908 def cost_fun(p): 909 R = _cd.uarray(_np.concatenate(( 910 D47_u - D47_interp(p['D47eq'].value), 911 D48_u - D48_interp(p['D47eq'].value), 912 ))) 913 914 invS = _np.linalg.inv(R.covar) 915 L = _cholesky(invS) 916 917 return L @ R.n 918 919 minresult = _lmfit.minimize( 920 cost_fun, 921 params, 922 method = 'least_squares', 923 scale_covar = False, 924 jac = '3-point', 925 ) 926 # slower but yields very similar results: 927 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 928 929 return minresult.params['D47eq'].value 930 931 wrapped_fun = _uc.wrap(fun) 932 D47eq[i] = wrapped_fun(D47[i], D48[i], *self.D47_coefs, *self.D48_coefs) 933 934 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 935 936 return D47eq, D48eq, p
Computes a correldata.uarray of equilibrium Δ47 values, each of which is
the closest (in the OGLS sense) to one (Δ47, Δ48) observation
considered independently of the others.
Also returns an array of corresponding p-values taking into account errors in Δ47 and Δ48 (and any covariance between the two) as well as errors in the Δ47 and Δ48 calibrations.
Note
This is both the fastest and the strongly recommended version of this calculation.
It is expected to yield an uarray with reasonably accurate covariance between the
D47eq values, but also between D47eq and all other variables.
def
joint_nearest_D47eq( self, D47: correldata.uarray, D48: correldata.uarray, ignore_calib_uncertainties: bool = False):
938 def joint_nearest_D47eq( 939 self, 940 D47: _cd.uarray, 941 D48: _cd.uarray, 942 ignore_calib_uncertainties: bool = False, 943 ): 944 """ 945 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 946 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 947 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 948 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 949 Δ<sub>48</sub> calibrations. 950 951 > [!CAUTION] 952 > Caution: the use of this function is **not generally recommended** except for 953 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 954 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 955 > as done by the recommended function `nearest_D47eq()`. 956 957 This is the most complete but slowest and not recommended version of this calculation. 958 It is expected to yield an `uarray` with reasonably accurate covariance between the 959 `D47eq` values, but also between `D47eq` and all other variables. 960 961 A faster but incomplete and potentially less accurate version of this calculation is 962 provided by `lazy_joint_nearest_D47eq()`. 963 """ 964 965 N = D47.size 966 N47 = self.D47_coefs.size 967 N48 = self.D48_coefs.size 968 969 def fun(j, *args): 970 971 args = _np.array(args) 972 D47_n = args[:N] 973 D48_n = args[N:2*N] 974 D47_calib_coefs_n = args[-N48-N47:-N48] 975 D48_calib_coefs_n = args[-N48:] 976 977 params = _lmfit.Parameters() 978 for k in range(N): 979 params.add(f'D47eq{k}', value = D47_n[k]) 980 981 D47_u = _cd.uarray(_uc.correlated_values(D47_n, D47.covar)) 982 D48_u = _cd.uarray(_uc.correlated_values(D48_n, D48.covar)) 983 D47_calib_coefs_u = _cd.uarray(_uc.correlated_values(D47_calib_coefs_n, self.D47_coefs.covar)) 984 D48_calib_coefs_u = _cd.uarray(_uc.correlated_values(D48_calib_coefs_n, self.D48_coefs.covar)) 985 986 D47i = D4x_calib_function( 987 self.interp.T, 988 D47_calib_coefs_u, 989 return_without_uncertainties = False, 990 ignore_calib_uncertainties = ignore_calib_uncertainties, 991 ) 992 993 D48i = D4x_calib_function( 994 self.interp.T, 995 D48_calib_coefs_u, 996 return_without_uncertainties = False, 997 ignore_calib_uncertainties = ignore_calib_uncertainties, 998 ) 999 1000 D47_interp = uarray_compatible_interp(D47i.n, D47i) 1001 D48_interp = uarray_compatible_interp(D47i.n, D48i) 1002 1003 def cost_fun(p): 1004 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1005 R = _cd.uarray(_np.concatenate(( 1006 D47_u - D47_interp(_D47eq), 1007 D48_u - D48_interp(_D47eq), 1008 ))) 1009 1010 invS = _np.linalg.inv(R.covar) 1011 L = _cholesky(invS) 1012 1013 # print(((L @ R.n)**2).sum()) 1014 return L @ R.n 1015 1016 minresult = _lmfit.minimize( 1017 cost_fun, 1018 params, 1019 method = 'least_squares', 1020 scale_covar = False, 1021 jac = '3-point', 1022 ) 1023 # slower but yields very similar results: 1024 # minresult = _lmfit.minimize(cost_fun, params, method = 'powell', scale_covar = False) 1025 1026 return minresult.params[f'D47eq{j}'].value 1027 1028 wrapped_fun = _uc.wrap(fun) 1029 1030 D47eq = _cd.uarray([wrapped_fun(j, *D47, *D48, *self.D47_coefs, *self.D48_coefs) for j in range(N)]) 1031 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1032 1033 return D47eq, D48eq, p
Returns a correldata.uarray of equilibrium Δ47 values which are jointly closest (in the OGLS sense)
to a sequence of (Δ47, Δ48) pairs. Also returns an array of
corresponding p-values taking into account errors in Δ47 and Δ48
(and any covariance between the two) as well as errors in the Δ47 and
Δ48 calibrations.
Caution
Caution: the use of this function is not generally recommended except for
experimentation purposes, because it is conceptually and numerically risky to jointly
fit the sequence of Teq values, as opposed to fitting each of them individually,
as done by the recommended function nearest_D47eq().
This is the most complete but slowest and not recommended version of this calculation.
It is expected to yield an uarray with reasonably accurate covariance between the
D47eq values, but also between D47eq and all other variables.
A faster but incomplete and potentially less accurate version of this calculation is
provided by lazy_joint_nearest_D47eq().
def
lazy_joint_nearest_D47eq( self, D47: correldata.uarray, D48: correldata.uarray, ignore_calib_uncertainties: bool = False):
1035 def lazy_joint_nearest_D47eq( 1036 self, 1037 D47: _cd.uarray, 1038 D48: _cd.uarray, 1039 ignore_calib_uncertainties: bool = False, 1040 ): 1041 """ 1042 Returns a `correldata.uarray` of equilibrium Δ<sub>47</sub> values which are *jointly* closest (in the OGLS sense) 1043 to a sequence of (Δ<sub>47</sub>, Δ<sub>48</sub>) pairs. Also returns an array of 1044 corresponding p-values taking into account errors in Δ<sub>47</sub> and Δ<sub>48</sub> 1045 (and any covariance between the two) as well as errors in the Δ<sub>47</sub> and 1046 Δ<sub>48</sub> calibrations. 1047 1048 > [!CAUTION] 1049 > Caution: the use of this function is **not generally recommended** except for 1050 > experimentation purposes, because it is conceptually and numerically risky to *jointly* 1051 > fit the sequence of `Teq` values, as opposed to fitting each of them individually, 1052 > as done by the recommended function `nearest_D47eq()`. 1053 1054 This is a faster but incomplete version of this calculation. It is expected to yield an 1055 `uarray` with roughly accurate covariance between the `Teq` values, but without computing 1056 the covariance with any other variables. 1057 1058 A slower but complete and more accurate version of this calculation is provided by 1059 `joint_nearest_D47eq()`. 1060 """ 1061 1062 N = D47.size 1063 1064 params = _lmfit.Parameters() 1065 for k in range(N): 1066 params.add(f'D47eq{k}', value = D47[k].n) 1067 1068 def cost_fun(p, ignore_calib_uncertainties = ignore_calib_uncertainties): 1069 _D47eq = _np.array([p[f'D47eq{k}'] for k in range(N)]) 1070 1071 if ignore_calib_uncertainties: 1072 R = _cd.uarray(_np.concatenate(( 1073 D47 - self.D47u_as_function_of_D47n(_D47eq).n, 1074 D48 - self.D48u_as_function_of_D47n(_D47eq).n, 1075 ))) 1076 else: 1077 R = _cd.uarray(_np.concatenate(( 1078 D47 - self.D47u_as_function_of_D47n(_D47eq), 1079 D48 - self.D48u_as_function_of_D47n(_D47eq), 1080 ))) 1081 1082 invS = _np.linalg.inv(R.covar) 1083 L = _cholesky(invS) 1084 1085 # print(((L @ R.n)**2).sum()) 1086 return L @ R.n 1087 1088 minresult = _lmfit.minimize( 1089 cost_fun, 1090 params, 1091 method = 'least_squares', 1092 scale_covar = False, 1093 jac = '3-point', 1094 ) 1095 1096 D47eq = _cd.uarray([minresult.uvars[f'D47eq{k}'] for k in range(N)]) 1097 1098 p, D48eq = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47eq, ignore_calib_uncertainties = ignore_calib_uncertainties) 1099 1100 return D47eq, D48eq, p
Returns a correldata.uarray of equilibrium Δ47 values which are jointly closest (in the OGLS sense)
to a sequence of (Δ47, Δ48) pairs. Also returns an array of
corresponding p-values taking into account errors in Δ47 and Δ48
(and any covariance between the two) as well as errors in the Δ47 and
Δ48 calibrations.
Caution
Caution: the use of this function is not generally recommended except for
experimentation purposes, because it is conceptually and numerically risky to jointly
fit the sequence of Teq values, as opposed to fitting each of them individually,
as done by the recommended function nearest_D47eq().
This is a faster but incomplete version of this calculation. It is expected to yield an
uarray with roughly accurate covariance between the Teq values, but without computing
the covariance with any other variables.
A slower but complete and more accurate version of this calculation is provided by
joint_nearest_D47eq().
def
projected_D47eq( self, D47: correldata.uarray, D48: correldata.uarray, kinetic_slope: float | uncertainties.core.AffineScalarFunc):
1102 def projected_D47eq( 1103 self, 1104 D47: _cd.uarray, 1105 D48: _cd.uarray, 1106 kinetic_slope: (float | _uc.UFloat), 1107 ): 1108 """ 1109 Projects one or more (Δ<sub>47</sub>, Δ<sub>48</sub>) observations onto the equlibrium curve 1110 following a kinetic fractionation vector with a given slope (∂Δ<sub>48</sub>/∂Δ<sub>47</sub>). 1111 1112 **Arguments** 1113 * `D47`: observed Δ<sub>47</sub> value(s) 1114 * `D48`: observed Δ<sub>48</sub> value(s) 1115 * `kinetic_slope`: kinetic fractionation slopw, with or without uncertainty 1116 1117 Returns a tuple of uarrays corresponding to the projected Δ<sub>47</sub> and Δ<sub>48</sub> values. 1118 1119 > [!NOTE] 1120 > This is not a least-squares minimization problem but a direct calculation, and should thus 1121 > be much faster than the various `CorelData.nearestD47eq()` methods. 1122 """ 1123 1124 D47 = _cd.uarray(D47) 1125 D48 = _cd.uarray(D48) 1126 N = D47.size 1127 N47c = self.D47_coefs.size 1128 N48c = self.D48_coefs.size 1129 D47p = D47 * 0 1130 1131 for i in range(N): 1132 1133 # function to solve 1134 def fun(x, *args): # args = (D47, D48, kinetic_slope, *self.D47_coefs, *self.D48_coefs) 1135 1136 args = _np.array(args) 1137 D47_n = args[0] 1138 D48_n = args[1] 1139 kslope_n = args[2] 1140 D47_calib_coefs_n = args[-N48c-N47c:-N48c] 1141 D48_calib_coefs_n = args[-N48c:] 1142 1143 D47i = D4x_calib_function( 1144 self.interp.T, 1145 D47_calib_coefs_n, 1146 return_without_uncertainties = False, 1147 ) 1148 1149 D48i = D4x_calib_function( 1150 self.interp.T, 1151 D48_calib_coefs_n, 1152 return_without_uncertainties = False, 1153 ) 1154 1155 D48_interp = uarray_compatible_interp(D47i, D48i) 1156 1157 return D48_n - D48_interp(x) - kslope_n * (D47_n - x) 1158 1159 def g(*args): 1160 return _fsolve(fun, [100.], args = args)[0] 1161 1162 wg = _uc.wrap(g) 1163 1164 D47p[i] = wg( 1165 D47[i], 1166 D48[i], 1167 kinetic_slope, 1168 *self.D47_coefs, 1169 *self.D48_coefs, 1170 ) 1171 1172 _, D48p = self._compute_p_and_D48eq_from_D47eq(D47, D48, D47p, ignore_calib_uncertainties = False) 1173 1174 return D47p, D48p
Projects one or more (Δ47, Δ48) observations onto the equlibrium curve following a kinetic fractionation vector with a given slope (∂Δ48/∂Δ47).
Arguments
D47: observed Δ47 value(s)D48: observed Δ48 value(s)kinetic_slope: kinetic fractionation slopw, with or without uncertainty
Returns a tuple of uarrays corresponding to the projected Δ47 and Δ48 values.
Note
This is not a least-squares minimization problem but a direct calculation, and should thus
be much faster than the various CorelData.nearestD47eq() methods.
def
Teq_pdf( self, D47: <function ufloat>, Tmin: float | None = None, Tmax: float | None = None, Tinc: float = 0.2, default_D47_sigmas: float = 4.0, ignore_calib_uncertainties: bool = False, run_qmc: bool = False, N_qmc: int = 1024):
1176 def Teq_pdf( 1177 self, 1178 D47: _uc.ufloat, 1179 Tmin: (float | None) = None, 1180 Tmax: (float | None) = None, 1181 Tinc: float = 0.2, 1182 default_D47_sigmas: float = 4.0, 1183 ignore_calib_uncertainties: bool = False, 1184 run_qmc: bool = False, 1185 N_qmc: int = 1024, 1186 ): 1187 """ 1188 Compute the unit-normalized probability distribution function (PDF) of the 1189 equilibrium temperature (`Teq`) for a given (`UFloat`) value of Δ<sub>47</sub>. 1190 1191 **Arguments** 1192 * `D47`: Δ<sub>47</sub> value (with uncertainty) 1193 * `Tmin`: minimum temperature over which to compute the PDF; if not specified, 1194 use temperature corresponding to `D47.n + `default_D47_sigmas` * D47.s` 1195 * `Tmax`: maximum temperature over which to compute the PDF; if not specified, 1196 use temperature corresponding to `D47.n - `default_D47_sigmas` * D47.s` 1197 * `Tinc`: temperature increment over which to compute the PDF 1198 * `default_D47_sigmas`: see `Tmin` and `Tmin` above 1199 * `ignore_calib_uncertainties`: whether to propagate calibration uncertainties 1200 * `run_qmc`: whether to also run a Quasi Monte carlo simulation to estimate the PDF 1201 * `N_qmc`: number of iterations in the above Quasi Monte Carlo simulation 1202 1203 **Returns** 1204 * `Ti`: Evenly-spaced array of temperature values over which the PDF is computed 1205 * `pdf`: PDF evaluated over `Ti` 1206 * `Tqmc` (only returned if `run_qmc = True`): array of `N_qmc` temperature values 1207 computed in the Quasi Monte Carlo simulation 1208 """ 1209 1210 if Tmin is None: 1211 Tmin = _np.floor(self.T_as_function_of_D47( 1212 D47.n + default_D47_sigmas * D47.s, 1213 ignore_calib_uncertainties = ignore_calib_uncertainties, 1214 ).n) 1215 1216 if Tmax is None: 1217 Tmax = _np.ceil(self.T_as_function_of_D47( 1218 D47.n - default_D47_sigmas * D47.s, 1219 ignore_calib_uncertainties = ignore_calib_uncertainties, 1220 ).n) 1221 1222 assert Tmin < Tmax, "Tmax must be strictly greater than Tmin" 1223 assert Tinc > 0, "Tinc must be strictly greater than zero" 1224 1225 # compute interpolated Ti values 1226 Ti = _np.arange(Tmin, Tmax+Tinc, Tinc) 1227 1228 pdf = transform_pdf_monotonic( 1229 f_inv = lambda T: D4x_calib_function( 1230 T, 1231 self.D47_coefs, 1232 return_without_uncertainties = ignore_calib_uncertainties, 1233 ignore_calib_uncertainties = ignore_calib_uncertainties, 1234 ), 1235 df_inv = lambda T: D4x_calib_derivative( 1236 T, 1237 self.D47_coefs, 1238 return_without_uncertainties = ignore_calib_uncertainties, 1239 ignore_calib_uncertainties = ignore_calib_uncertainties, 1240 ), 1241 mu_x = D47.n, 1242 sigma_x = D47.s, 1243 yi = Ti, 1244 ) 1245 1246 if run_qmc: 1247 1248 from scipy.stats import qmc 1249 from tqdm.rich import tqdm 1250 1251 #parameters to jiggle 1252 input_params = _cd.uarray([D47, *self.D47_coefs]) 1253 1254 # QMC sampler for the correlation matrix of these parameters 1255 qmc_dist = qmc.MultivariateNormalQMC( 1256 mean = input_params.n*0, 1257 cov = input_params.cor, 1258 ) 1259 1260 # QMC samples 1261 qmc_draws = input_params.n + qmc_dist.random(N_qmc) * input_params.s 1262 1263 # initialize T_qmc 1264 Tqmc = _cd.uarray(_np.zeros((N_qmc,))) 1265 1266 for k in tqdm(range(N_qmc)): 1267 # jiggled D47 and D47coefs 1268 _D47 = qmc_draws[k,0] 1269 if ignore_calib_uncertainties: 1270 _coefs = self.D47_coefs 1271 else: 1272 _coefs = _cd.uarray(_uc.correlated_values(qmc_draws[k,1:], self.D47_coefs.covar)) 1273 1274 # jiggled D47 1275 _D47i = D4x_calib_function(self.interp.T, _coefs) 1276 _f = uarray_compatible_interp(_D47i.n, self.interp.T) 1277 Tqmc[k] = _f(_D47) 1278 1279 return Ti, pdf, Tqmc 1280 1281 return Ti, pdf
Compute the unit-normalized probability distribution function (PDF) of the
equilibrium temperature (Teq) for a given (UFloat) value of Δ47.
Arguments
D47: Δ47 value (with uncertainty)Tmin: minimum temperature over which to compute the PDF; if not specified, use temperature corresponding toD47.n +default_D47_sigmas* D47.sTmax: maximum temperature over which to compute the PDF; if not specified, use temperature corresponding toD47.n -default_D47_sigmas* D47.sTinc: temperature increment over which to compute the PDFdefault_D47_sigmas: seeTminandTminaboveignore_calib_uncertainties: whether to propagate calibration uncertaintiesrun_qmc: whether to also run a Quasi Monte carlo simulation to estimate the PDFN_qmc: number of iterations in the above Quasi Monte Carlo simulation
Returns
Ti: Evenly-spaced array of temperature values over which the PDF is computedpdf: PDF evaluated overTiTqmc(only returned ifrun_qmc = True): array ofN_qmctemperature values computed in the Quasi Monte Carlo simulation
def
save_Teq_report( X, Y, T, p, filename, Xname='D47', Yname='D48', Tname='T95', labelname='Sample', fmt_Xnv='.4f', fmt_Xse='.4f', fmt_Ynv='.4f', fmt_Yse='.4f', fmt_Tnv='.1f', fmt_Tse='.1f', fmt_cm='.6f', fmt_pv='.2e', labels=None, sep=',', p_cutoff=0.05):
1287def save_Teq_report( 1288 X, 1289 Y, 1290 T, 1291 p, 1292 filename, 1293 Xname = 'D47', 1294 Yname = 'D48', 1295 Tname = 'T95', 1296 labelname = 'Sample', 1297 fmt_Xnv = '.4f', 1298 fmt_Xse = '.4f', 1299 fmt_Ynv = '.4f', 1300 fmt_Yse = '.4f', 1301 fmt_Tnv = '.1f', 1302 fmt_Tse = '.1f', 1303 fmt_cm = '.6f', 1304 fmt_pv = '.2e', 1305 labels = None, 1306 sep = ',', 1307 p_cutoff = 0.05, 1308): 1309 """ 1310 Save a temperature report to a csv file. 1311 Includes observed `D47`, `D48`, p-equilibrium values, and nearest `Teq` with sensible precision defaults. 1312 Alternatively, users may find [`correldata.CorrelData.str()`](https://mdaeron.github.io/correldata/#CorrelData.str) 1313 to be more versatile. 1314 """ 1315 N = T.size 1316 if labels is None: 1317 labels = [str(k+1) for k in range(N)] 1318 1319 with open(filename, 'w') as fid: 1320 fid.write(f'{labelname}{sep}{Xname}{sep}SE{sep}correl{sep*N}{Yname}{sep}SE{sep}correl{sep*N}p-value{sep}{Tname}{sep}SE{sep}correl') 1321 Xnv = _unp.nominal_values(X) 1322 Xse = _unp.std_devs(X) 1323 Xcm = _np.array(_uc.correlation_matrix(X)) 1324 Ynv = _unp.nominal_values(Y) 1325 Yse = _unp.std_devs(Y) 1326 Ycm = _np.array(_uc.correlation_matrix(Y)) 1327 Tnv = _unp.nominal_values(T) 1328 Tse = _unp.std_devs(T) 1329 Tcm = _np.array(_uc.correlation_matrix(T)) 1330 for k in range(X.size): 1331 fid.write(f'\n{labels[k]}{sep}{Xnv[k]:{fmt_Xnv}}{sep}{Xse[k]:{fmt_Xse}}{sep}') 1332 fid.write(sep.join([f'{Xcm[j,k]:{fmt_cm}}' for j in range(N)])) 1333 fid.write(f'{sep}{Ynv[k]:{fmt_Ynv}}{sep}{Yse[k]:{fmt_Yse}}{sep}') 1334 fid.write(sep.join([f'{Ycm[j,k]:{fmt_cm}}' for j in range(N)])) 1335 fid.write(f'{sep}{p[k]:{fmt_pv}}') 1336 if p[k] >= p_cutoff: 1337 fid.write(f'{sep}{Tnv[k]:{fmt_Tnv}}{sep}{Tse[k]:{fmt_Tse}}{sep}') 1338 fid.write(sep.join([f'{Tcm[j,k]:{fmt_cm}}' for j in range(N)]))
Save a temperature report to a csv file.
Includes observed D47, D48, p-equilibrium values, and nearest Teq with sensible precision defaults.
Alternatively, users may find correldata.CorrelData.str()
to be more versatile.
def
confidence_band( t: ArrayLike, fx: Callable, fy: Callable, p: float = 0.95, dt: float = 1e-09):
9def confidence_band( 10 t: ArrayLike, 11 fx: Callable, 12 fy: Callable, 13 p: float = 0.95, 14 dt: float = 1e-9, 15): 16 """ 17 Return an (N, 2) array of (x, y) vertices outlining a confidence region, at a given p-value, 18 for the central parametric curve ***C*** defined by `x = fx(t)` and `y = fy(t)`. 19 20 This confidence region is defined as the union of confidence ellipses for all points along ***C***. 21 22 **Arguments** 23 * `t`: array of values over which to sample ***C*** 24 * `fx`: parametric function of `t` yielding x values of ***C*** as 25 [UFloat](https://pythonhosted.org/uncertainties/tech_guide.html) values 26 * `fy`: parametric function of `t` yielding y values of ***C*** as 27 [UFloat](https://pythonhosted.org/uncertainties/tech_guide.html) values 28 * `p`: p-value for the confidence region to return 29 * `dt`: `t` scale at which to evaluate derivatives 30 31 Returns a (N, 2) array of (x, y) vertices. 32 """ 33 34 # curve position & covariance 35 curve = lambda _t: np.array([fx(_t).n, fy(_t).n]) 36 def covariance(_t): 37 return np.array(covariance_matrix((fx(_t), fy(_t)))) 38 # corresponding derivatives 39 def deriv(_f, _t, _dt = dt): 40 return (_f(float(_t) + _dt) - _f(float(_t) - _dt)) / (2 * _dt) 41 mu_dot = lambda _t: deriv(curve, _t) 42 sigma_dot = lambda _t: deriv(covariance, _t) 43 44 # ellipse discretization 45 def ellipse_points(mean, cov, chi2_val, n_pts = 120): 46 phi = np.linspace(0, 2 * np.pi, n_pts, endpoint=False) 47 unit = np.stack([np.cos(phi), np.sin(phi)], axis = 1) 48 L = np.linalg.cholesky(cov) 49 return mean + np.sqrt(chi2_val) * (unit @ L.T) 50 51 # find angular positions where a given ellipse is tangent to the union of ellipses 52 def envelope_contact_angles(t, chi2_val, n_pts = 2000): 53 mu = curve(t) 54 Sigma = covariance(t) 55 L = np.linalg.cholesky(Sigma) 56 s = np.sqrt(chi2_val) 57 58 Lambda = np.linalg.inv(Sigma) 59 Sigma_d = sigma_dot(t) 60 Lambda_dot = -Lambda @ Sigma_d @ Lambda 61 mu_d = mu_dot(t) 62 63 phi = np.linspace(0, 2 * np.pi, n_pts, endpoint = False) 64 u = np.stack([np.cos(phi), np.sin(phi)], axis = 1) 65 delta = s * (u @ L.T) 66 67 term1 = -2.0 * (delta @ (Lambda @ mu_d)) 68 term2 = np.einsum('ni,ij,nj->n', delta, Lambda_dot, delta) 69 dFdt = term1 + term2 70 71 signs = np.sign(dFdt) 72 crossings = np.where(np.diff(signs) != 0)[0] 73 74 contact_pts = [] 75 for idx in crossings: 76 phi0, phi1 = phi[idx], phi[idx + 1] 77 f0, f1 = dFdt[idx], dFdt[idx + 1] 78 phi_c = phi0 - f0 * (phi1 - phi0) / (f1 - f0) 79 u_c = np.array([np.cos(phi_c), np.sin(phi_c)]) 80 pt = mu + s * L @ u_c 81 contact_pts.append(pt) 82 83 return contact_pts 84 85 # build the upper and lower limits of the envelope 86 def build_envelope(ts, chi2_val, means): 87 all_contacts = [] 88 all_t = [] 89 90 for i, t in enumerate(ts): 91 pts = envelope_contact_angles(t, chi2_val) 92 for pt in pts: 93 all_contacts.append(pt) 94 all_t.append(i) 95 96 if not all_contacts: 97 return None, None 98 99 pts = np.array(all_contacts) 100 t_idx = np.array(all_t) 101 102 upper, lower = [], [] 103 104 for i, t in enumerate(ts): 105 mask = t_idx == i 106 pts_t = pts[mask] 107 if len(pts_t) == 0: 108 continue 109 110 i0, i1 = max(0, i - 1), min(len(ts) - 1, i + 1) 111 tangent = means[i1] - means[i0] 112 normal = np.array([-tangent[1], tangent[0]]) 113 114 for pt in pts_t: 115 side = np.dot(pt - means[i], normal) 116 if side >= 0: 117 upper.append((i, pt)) 118 else: 119 lower.append((i, pt)) 120 121 upper.sort(key=lambda x: x[0]) 122 lower.sort(key=lambda x: x[0]) 123 124 upper_pts = np.array([p for _, p in upper]) 125 lower_pts = np.array([p for _, p in lower]) 126 127 return upper_pts, lower_pts 128 129 # Trace the arc of the terminal ellipse that faces outward, running exactly 130 # from upper_end to lower_end along the outward-facing side. 131 # Strategy: parametrise the full ellipse by angle, find the angles 132 # corresponding to upper_end and lower_end, then extract the arc between 133 # them that passes through the outward direction. 134 def terminal_cap(mean, cov, chi2_val, outward_tangent, upper_end, lower_end, n_pts = 200): 135 L = np.linalg.cholesky(cov) 136 Linv = np.linalg.inv(L) 137 s = np.sqrt(chi2_val) 138 139 # Map upper_end and lower_end back to angles in the unit circle 140 def point_to_angle(pt): 141 u = Linv @ (pt - mean) / s 142 return np.arctan2(u[1], u[0]) 143 144 phi_upper = point_to_angle(upper_end) 145 phi_lower = point_to_angle(lower_end) 146 phi_out = np.arctan2(outward_tangent[1], outward_tangent[0]) 147 148 # Normalise all angles relative to phi_upper, on [0, 2π) 149 def normalise(phi, ref): 150 return (phi - ref) % (2 * np.pi) 151 152 phi_lower_n = normalise(phi_lower, phi_upper) 153 phi_out_n = normalise(phi_out, phi_upper) 154 155 # The outward arc from phi_upper to phi_lower passes through phi_out. 156 # Determine direction: if phi_out_n < phi_lower_n, the outward arc goes 157 # forward (increasing angle); otherwise it goes backward. 158 if phi_out_n < phi_lower_n: 159 # Forward arc: phi_upper → phi_upper + phi_lower_n 160 phis = np.linspace(phi_upper, phi_upper + phi_lower_n, n_pts) 161 else: 162 # Backward arc: phi_upper → phi_upper - (2π - phi_lower_n) 163 phis = np.linspace(phi_upper, phi_upper - (2 * np.pi - phi_lower_n), n_pts) 164 165 u = np.stack([np.cos(phis), np.sin(phis)], axis=1) 166 arc = mean + s * (u @ L.T) 167 168 return arc 169 170 chi2_value = chi2.ppf(p, df = 2) 171 means = curve(t).T 172 covs = np.array([covariance(_) for _ in t]) 173 upper, lower = build_envelope(t, chi2_value, means) 174 175 # Outward tangents at each tip: unit vector pointing away from curve interior 176 tangent_start = means[0] - means[1] 177 tangent_end = means[-1] - means[-2] 178 179 cap_start = terminal_cap( 180 means[0], covs[0], chi2_value, tangent_start, 181 upper_end=lower[0], # polygon arrives via lower[::-1], which ends at lower[0] 182 lower_end=upper[0], # polygon departs via upper, which starts at upper[0] 183 ) 184 cap_end = terminal_cap( 185 means[-1], covs[-1], chi2_value, tangent_end, 186 upper_end=upper[-1], # polygon arrives via upper, which ends at upper[-1] 187 lower_end=lower[-1], # polygon departs via lower[::-1], which starts at lower[-1] 188 ) 189 190 band_x = np.concatenate([ 191 upper[:, 0], 192 cap_end[:, 0], 193 lower[::-1, 0], 194 cap_start[:, 0], 195 ]) 196 band_y = np.concatenate([ 197 upper[:, 1], 198 cap_end[:, 1], 199 lower[::-1, 1], 200 cap_start[:, 1], 201 ]) 202 203 return np.array([band_x, band_y]).T
Return an (N, 2) array of (x, y) vertices outlining a confidence region, at a given p-value,
for the central parametric curve C defined by x = fx(t) and y = fy(t).
This confidence region is defined as the union of confidence ellipses for all points along C.
Arguments
t: array of values over which to sample Cfx: parametric function oftyielding x values of C as UFloat valuesfy: parametric function oftyielding y values of C as UFloat valuesp: p-value for the confidence region to returndt:tscale at which to evaluate derivatives
Returns a (N, 2) array of (x, y) vertices.