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