From 2e8a680cfde2df8800e185cb1feb41a7275ae456 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 14 Aug 2023 18:00:46 +0200 Subject: [PATCH 01/22] User-friendly error when setting manually ds.sections Solves #200 --- src/dtscalibration/datastore.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 01faf0c6..f2ac8215 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -213,6 +213,12 @@ def sections(self): def sections(self): self.attrs["_sections"] = yaml.dump(None) + @sections.setter + def sections(self, value): + msg = "Not possible anymore. Instead, pass the sections as an argument to \n" \ + "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." + raise NotImplementedError(msg) + def check_reference_section_values(self): """ Checks if the values of the used sections are of the right datatype From 07270e371de5ac3bb8555b58dc4b28e12b54b4dd Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 14 Aug 2023 18:04:45 +0200 Subject: [PATCH 02/22] User-friendly error when setting manually ds.sections Solves #200 --- src/dtscalibration/datastore.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index f2ac8215..0d81b2ce 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -219,6 +219,12 @@ def sections(self, value): "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." raise NotImplementedError(msg) + @sections.setter + def sections(self, value): + msg = "Not possible anymore. Instead, pass the sections as an argument to \n" \ + "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." + raise NotImplementedError(msg) + def check_reference_section_values(self): """ Checks if the values of the used sections are of the right datatype From 66aa9f00e335ebe90774211052be3f749498f972 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 14 Aug 2023 18:05:55 +0200 Subject: [PATCH 03/22] Revert accidental double implementation of sections.setter --- src/dtscalibration/datastore.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 0d81b2ce..f2ac8215 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -219,12 +219,6 @@ def sections(self, value): "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." raise NotImplementedError(msg) - @sections.setter - def sections(self, value): - msg = "Not possible anymore. Instead, pass the sections as an argument to \n" \ - "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." - raise NotImplementedError(msg) - def check_reference_section_values(self): """ Checks if the values of the used sections are of the right datatype From 10dfbb42df0d785b39d98d44cd5fa1701bdcf18f Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Thu, 17 Aug 2023 10:12:22 +0200 Subject: [PATCH 04/22] Update documentation of datastore class --- src/dtscalibration/datastore.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index f2ac8215..7033bc68 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -876,9 +876,6 @@ def variance_stokes_constant(self, st_label, sections=None, reshape_residuals=Tr dtscalibration/python-dts-calibration/blob/main/examples/notebooks/\ 04Calculate_variance_Stokes.ipynb>`_ """ - - # var_I, resid = variance_stokes_constant_util(st_label, sections=None, reshape_residuals=True) - # def variance_stokes_constant_util(st_label, sections=None, reshape_residuals=True): def func_fit(p, xs): return p[:xs, None] * p[None, xs:] @@ -1020,7 +1017,12 @@ def variance_stokes_exponential( Parameters ---------- - reshape_residuals + suppress_info : bool, optional + Suppress print statements. + use_statsmodels : bool, optional + Use statsmodels to fit the exponential. If `False`, use scipy. + reshape_residuals : bool, optional + Reshape the residuals to the shape of the Stokes intensity st_label : str label of the Stokes, anti-Stokes measurement. E.g., st, ast, rst, rast @@ -1706,6 +1708,9 @@ def calibration_single_ended( variance of the estimate of dalpha. Covariances between alpha and other parameters are not accounted for. + fix_alpha : Tuple[array-like, array-like], optional + A tuple containing two array-likes. The first array-like is the integrated + differential attenuation of length x, and the second item is its variance. Returns ------- @@ -3035,7 +3040,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): arr = da.concatenate( ( da.zeros( - (1, i_splice, 1), chunks=((1, i_splice, 1)), dtype=bool + (1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool ), da.ones( (1, no - i_splice, 1), @@ -3051,7 +3056,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): da.ones((1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool), da.zeros( (1, no - i_splice, 1), - chunks=((1, no - i_splice, 1)), + chunks=(1, no - i_splice, 1), dtype=bool, ), ), @@ -3832,7 +3837,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): arr = da.concatenate( ( da.zeros( - (1, i_splice, 1), chunks=((1, i_splice, 1)), dtype=bool + (1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool ), da.ones( (1, no - i_splice, 1), @@ -3848,7 +3853,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): da.ones((1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool), da.zeros( (1, no - i_splice, 1), - chunks=((1, no - i_splice, 1)), + chunks=(1, no - i_splice, 1), dtype=bool, ), ), From 77d9ede96020e6347feb2a8626a45d4e529624e3 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 09:16:56 +0200 Subject: [PATCH 05/22] Removed numbagg from dependencies llvmlite is a pain to install with pip. Alternatively, it could be installed with conda --- .github/workflows/python-publish-dry-run.yml | 2 +- pyproject.toml | 18 ++++++++++-------- src/dtscalibration/datastore.py | 15 +++++++-------- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/.github/workflows/python-publish-dry-run.yml b/.github/workflows/python-publish-dry-run.yml index 74f2b3a0..66c2d5de 100644 --- a/.github/workflows/python-publish-dry-run.yml +++ b/.github/workflows/python-publish-dry-run.yml @@ -34,5 +34,5 @@ jobs: pip install build - name: Build package run: python -m build - - name: Test build + - name: Check long description for PyPi run: twine check --strict dist/* diff --git a/pyproject.toml b/pyproject.toml index 315e0bb7..162b2448 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,13 +19,13 @@ disable = true # Requires confirmation when publishing to pypi. [project] name = "dtscalibration" -description = "A Python package to load raw DTS files, perform a calibration, and plot the result." +description = "Load Distributed Temperature Sensing (DTS) files, calibrate the temperature and estimate its uncertainty." readme = "README.rst" license = "BSD-3-Clause" requires-python = ">=3.9, <3.12" authors = [ - {email = "bdestombe@gmail.com"}, - {name = "Bas des Tombe, Bart Schilperoort"} + {name = "Bas des Tombe", email = "bdestombe@gmail.com"}, + {name = "Bart Schilperoort", email = "b.schilperoort@gmail.com"}, ] maintainers = [ {name = "Bas des Tombe", email = "bdestombe@gmail.com"}, @@ -52,10 +52,12 @@ classifiers = [ "Topic :: Utilities", ] dependencies = [ - "numpy>=1.22.4", # xarray: recommended to use >= 1.22 for full quantile method support + "numpy>=1.22.4", # >= 1.22 for full quantile method support in xarray "dask", "pandas", - "xarray[parallel,accel]", + "xarray[parallel]", # numbagg (llvmlite) is a pain to install + "bottleneck", # speed up Xarray + "flox", # speed up Xarray "pyyaml>=6.0.1", "xmltodict", "scipy", @@ -74,13 +76,13 @@ dev = [ "isort", "black[jupyter]", "mypy", - "types-PyYAML", # for pyyaml types + "types-PyYAML", # for pyyaml types "types-xmltodict", # for xmltodict types - "pandas-stubs", # for pandas types + "pandas-stubs", # for pandas types "pytest", "pytest-cov", "jupyter", - "nbformat", # Needed to run the tests + "nbformat", # Needed to run the tests ] docs = [ # Required for ReadTheDocs "IPython", diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 7033bc68..434868af 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -215,8 +215,10 @@ def sections(self): @sections.setter def sections(self, value): - msg = "Not possible anymore. Instead, pass the sections as an argument to \n" \ - "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." + msg = ( + "Not possible anymore. Instead, pass the sections as an argument to \n" + "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." + ) raise NotImplementedError(msg) def check_reference_section_values(self): @@ -876,6 +878,7 @@ def variance_stokes_constant(self, st_label, sections=None, reshape_residuals=Tr dtscalibration/python-dts-calibration/blob/main/examples/notebooks/\ 04Calculate_variance_Stokes.ipynb>`_ """ + def func_fit(p, xs): return p[:xs, None] * p[None, xs:] @@ -3039,9 +3042,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): if direction == "fw": arr = da.concatenate( ( - da.zeros( - (1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool - ), + da.zeros((1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool), da.ones( (1, no - i_splice, 1), chunks=(1, no - i_splice, 1), @@ -3836,9 +3837,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): if direction == "fw": arr = da.concatenate( ( - da.zeros( - (1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool - ), + da.zeros((1, i_splice, 1), chunks=(1, i_splice, 1), dtype=bool), da.ones( (1, no - i_splice, 1), chunks=(1, no - i_splice, 1), From c18f8eab3e9cc46fe6e7d4a54decef250c424a40 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 12:42:07 +0200 Subject: [PATCH 06/22] Explicitly pass sections to more functions --- pyproject.toml | 2 +- src/dtscalibration/calibrate_utils.py | 33 ++++++++++------- .../calibration/section_utils.py | 14 ++++---- src/dtscalibration/datastore.py | 36 ++++++++++++++----- src/dtscalibration/plot.py | 24 ++++++++++--- tests/test_dtscalibration.py | 36 ++++++++++++------- 6 files changed, 100 insertions(+), 45 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 162b2448..ccb756f8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,7 +55,7 @@ dependencies = [ "numpy>=1.22.4", # >= 1.22 for full quantile method support in xarray "dask", "pandas", - "xarray[parallel]", # numbagg (llvmlite) is a pain to install + "xarray[parallel]", # numbagg (llvmlite) is a pain to install with pip "bottleneck", # speed up Xarray "flox", # speed up Xarray "pyyaml>=6.0.1", diff --git a/src/dtscalibration/calibrate_utils.py b/src/dtscalibration/calibrate_utils.py index 6e90e208..5c343ac8 100644 --- a/src/dtscalibration/calibrate_utils.py +++ b/src/dtscalibration/calibrate_utils.py @@ -48,6 +48,7 @@ def parse_st_var(ds, st_var, st_label="st"): def calibration_single_ended_helper( self, + sections, st_var, ast_var, fix_alpha, @@ -67,6 +68,7 @@ def calibration_single_ended_helper( calc_cov = True split = calibration_single_ended_solver( self, + sections, st_var, ast_var, calc_cov=calc_cov, @@ -189,6 +191,7 @@ def calibration_single_ended_helper( def calibration_single_ended_solver( # noqa: MC0001 ds, + sections=None, st_var=None, ast_var=None, calc_cov=True, @@ -235,7 +238,7 @@ def calibration_single_ended_solver( # noqa: MC0001 """ # get ix_sec argsort so the sections are in order of increasing x - ix_sec = ds.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = ds.ufunc_per_section(sections=sections, x_indices=True, calc_per="all") ds_sec = ds.isel(x=ix_sec) x_sec = ds_sec["x"].values @@ -255,7 +258,7 @@ def calibration_single_ended_solver( # noqa: MC0001 # X \gamma # Eq.34 cal_ref = ds.ufunc_per_section( - label="st", ref_temp_broadcasted=True, calc_per="all" + sections=sections, label="st", ref_temp_broadcasted=True, calc_per="all" ) # cal_ref = cal_ref # sort by increasing x data_gamma = 1 / (cal_ref.T.ravel() + 273.15) # gamma @@ -485,6 +488,7 @@ def calibration_single_ended_solver( # noqa: MC0001 def calibration_double_ended_helper( self, + sections, st_var, ast_var, rst_var, @@ -503,6 +507,7 @@ def calibration_double_ended_helper( if fix_alpha or fix_gamma: split = calibration_double_ended_solver( self, + sections, st_var, ast_var, rst_var, @@ -515,6 +520,7 @@ def calibration_double_ended_helper( else: out = calibration_double_ended_solver( self, + sections, st_var, ast_var, rst_var, @@ -887,7 +893,7 @@ def calibration_double_ended_helper( [fix_gamma[0]], out[0][: 2 * nt], E_all_exact, - out[0][2 * nt + n_E_in_cal :], + out[0][2 * nt + n_E_in_cal:], ) ) p_val[1 + 2 * nt + split["ix_from_cal_match_to_glob"]] = out[0][ @@ -987,7 +993,7 @@ def calibration_double_ended_helper( p0_est = np.concatenate( ( split["p0_est"][: 1 + 2 * nt], - split["p0_est"][1 + 2 * nt + n_E_in_cal :], + split["p0_est"][1 + 2 * nt + n_E_in_cal:], ) ) X_E1 = sp.csr_matrix(([], ([], [])), shape=(nt * nx_sec, self.x.size)) @@ -1102,6 +1108,7 @@ def calibration_double_ended_helper( def calibration_double_ended_solver( # noqa: MC0001 ds, + sections=None, st_var=None, ast_var=None, rst_var=None, @@ -1167,7 +1174,7 @@ def calibration_double_ended_solver( # noqa: MC0001 ------- """ - ix_sec = ds.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = ds.ufunc_per_section(sections=sections, x_indices=True, calc_per="all") ds_sec = ds.isel(x=ix_sec) ix_alpha_is_zero = ix_sec[0] # per definition of E @@ -1187,7 +1194,7 @@ def calibration_double_ended_solver( # noqa: MC0001 rast_var=rast_var, ix_alpha_is_zero=ix_alpha_is_zero, ) - df_est, db_est = calc_df_db_double_est(ds, ix_alpha_is_zero, 485.0) + df_est, db_est = calc_df_db_double_est(ds, sections, ix_alpha_is_zero, 485.0) ( E, @@ -1196,7 +1203,7 @@ def calibration_double_ended_solver( # noqa: MC0001 Zero_d, Z_TA_fw, Z_TA_bw, - ) = construct_submatrices(nt, nx_sec, ds, ds.trans_att.values, x_sec) + ) = construct_submatrices(sections, nt, nx_sec, ds, ds.trans_att.values, x_sec) # y # Eq.41--45 y_F = np.log(ds_sec.st / ds_sec.ast).values.ravel() @@ -1822,7 +1829,7 @@ def construct_submatrices_matching_sections(x, ix_sec, hix, tix, nt, trans_att): ) -def construct_submatrices(nt, nx, ds, trans_att, x_sec): +def construct_submatrices(sections, nt, nx, ds, trans_att, x_sec): """Wrapped in a function to reduce memory usage. E is zero at the first index of the reference section (ds_sec) Constructing: @@ -1840,7 +1847,9 @@ def construct_submatrices(nt, nx, ds, trans_att, x_sec): # Z \gamma # Eq.47 cal_ref = np.array( - ds.ufunc_per_section(label="st", ref_temp_broadcasted=True, calc_per="all") + ds.ufunc_per_section( + sections=sections, label="st", ref_temp_broadcasted=True, calc_per="all" + ) ) data_gamma = 1 / (cal_ref.ravel() + 273.15) # gamma coord_gamma_row = np.arange(nt * nx, dtype=int) @@ -2213,7 +2222,7 @@ def calc_alpha_double( return E, E_var -def calc_df_db_double_est(ds, ix_alpha_is_zero, gamma_est): +def calc_df_db_double_est(ds, sections, ix_alpha_is_zero, gamma_est): Ifwx0 = np.log( ds.st.isel(x=ix_alpha_is_zero) / ds.ast.isel(x=ix_alpha_is_zero) ).values @@ -2221,9 +2230,9 @@ def calc_df_db_double_est(ds, ix_alpha_is_zero, gamma_est): ds.rst.isel(x=ix_alpha_is_zero) / ds.rast.isel(x=ix_alpha_is_zero) ).values ref_temps_refs = ds.ufunc_per_section( - label="st", ref_temp_broadcasted=True, calc_per="all" + sections=sections, label="st", ref_temp_broadcasted=True, calc_per="all" ) - ix_sec = ds.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = ds.ufunc_per_section(sections=sections, x_indices=True, calc_per="all") ref_temps_x0 = ( ref_temps_refs[ix_sec == ix_alpha_is_zero].flatten().compute() + 273.15 ) diff --git a/src/dtscalibration/calibration/section_utils.py b/src/dtscalibration/calibration/section_utils.py index 126d575b..63dd1e5a 100644 --- a/src/dtscalibration/calibration/section_utils.py +++ b/src/dtscalibration/calibration/section_utils.py @@ -135,7 +135,7 @@ def ufunc_per_section( 1. Calculate the variance of the residuals in the along ALL the\ reference sections wrt the temperature of the water baths - >>> tmpf_var = d.ufunc_per_section( + >>> tmpf_var = d.ufunc_per_section(sections, >>> func='var', >>> calc_per='all', >>> label='tmpf', @@ -144,7 +144,7 @@ def ufunc_per_section( 2. Calculate the variance of the residuals in the along PER\ reference section wrt the temperature of the water baths - >>> tmpf_var = d.ufunc_per_section( + >>> tmpf_var = d.ufunc_per_section(sections, >>> func='var', >>> calc_per='stretch', >>> label='tmpf', @@ -153,7 +153,7 @@ def ufunc_per_section( 3. Calculate the variance of the residuals in the along PER\ water bath wrt the temperature of the water baths - >>> tmpf_var = d.ufunc_per_section( + >>> tmpf_var = d.ufunc_per_section(sections, >>> func='var', >>> calc_per='section', >>> label='tmpf', @@ -161,7 +161,7 @@ def ufunc_per_section( 4. Obtain the coordinates of the measurements per section - >>> locs = d.ufunc_per_section( + >>> locs = d.ufunc_per_section(sections, >>> func=None, >>> label='x', >>> temp_err=False, @@ -170,7 +170,7 @@ def ufunc_per_section( 5. Number of observations per stretch - >>> nlocs = d.ufunc_per_section( + >>> nlocs = d.ufunc_per_section(sections, >>> func=len, >>> label='x', >>> temp_err=False, @@ -182,14 +182,14 @@ def ufunc_per_section( temperature (a timeseries) is broadcasted to the shape of self[\ label]. The self[label] is not used for anything else. - >>> temp_ref = d.ufunc_per_section( + >>> temp_ref = d.ufunc_per_section(sections, >>> label='st', >>> ref_temp_broadcasted=True, >>> calc_per='all') 7. x-coordinate index - >>> ix_loc = d.ufunc_per_section(x_indices=True) + >>> ix_loc = d.ufunc_per_section(sections, x_indices=True) Note diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 434868af..69eac39c 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -1754,7 +1754,9 @@ def calibration_single_ended( else: matching_indices = None - ix_sec = self.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = self.ufunc_per_section( + sections=sections, x_indices=True, calc_per="all" + ) assert not np.any(self.st.isel(x=ix_sec) <= 0.0), ( "There is uncontrolled noise in the ST signal. Are your sections" "correctly defined?" @@ -1767,6 +1769,7 @@ def calibration_single_ended( if method == "wls": p_cov, p_val, p_var = calibration_single_ended_helper( self, + sections, st_var, ast_var, fix_alpha, @@ -2172,7 +2175,9 @@ def calibration_double_ended( nx = self.x.size nt = self["time"].size nta = self.trans_att.size - ix_sec = self.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = self.ufunc_per_section( + sections=sections, x_indices=True, calc_per="all" + ) nx_sec = ix_sec.size assert self.st.dims[0] == "x", "Stokes are transposed" @@ -2218,6 +2223,7 @@ def calibration_double_ended( if method == "wls": p_cov, p_val, p_var = calibration_double_ended_helper( self, + sections, st_var, ast_var, rst_var, @@ -2899,6 +2905,7 @@ def average_single_ended( def average_double_ended( self, + sections=None, p_val="p_val", p_cov="p_cov", st_var=None, @@ -3083,6 +3090,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): pass self.conf_int_double_ended( + sections=sections, p_val=p_val, p_cov=p_cov, st_var=st_var, @@ -3695,6 +3703,7 @@ def conf_int_single_ended( def conf_int_double_ended( self, + sections=None, p_val="p_val", p_cov="p_cov", st_var=None, @@ -3947,7 +3956,9 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): p_cov = self[p_cov].values assert p_cov.shape == (npar, npar) - ix_sec = self.ufunc_per_section(x_indices=True, calc_per="all") + ix_sec = self.ufunc_per_section( + sections=sections, x_indices=True, calc_per="all" + ) nx_sec = ix_sec.size from_i = np.concatenate( ( @@ -4122,7 +4133,9 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): if var_only_sections: # sets the values outside the reference sections to NaN - xi = self.ufunc_per_section(x_indices=True, calc_per="all") + xi = self.ufunc_per_section( + sections=sections, x_indices=True, calc_per="all" + ) x_mask_ = [True if ix in xi else False for ix in range(self.x.size)] x_mask = np.reshape(x_mask_, (1, -1, 1)) self[label + "_mc_set"] = self[label + "_mc_set"].where(x_mask) @@ -4354,6 +4367,7 @@ def ufunc_per_section( reference sections wrt the temperature of the water baths >>> tmpf_var = d.ufunc_per_section( + >>> sections=sections, >>> func='var', >>> calc_per='all', >>> label='tmpf', @@ -4363,6 +4377,7 @@ def ufunc_per_section( reference section wrt the temperature of the water baths >>> tmpf_var = d.ufunc_per_section( + >>> sections=sections, >>> func='var', >>> calc_per='stretch', >>> label='tmpf', @@ -4372,6 +4387,7 @@ def ufunc_per_section( water bath wrt the temperature of the water baths >>> tmpf_var = d.ufunc_per_section( + >>> sections=sections, >>> func='var', >>> calc_per='section', >>> label='tmpf', @@ -4380,6 +4396,7 @@ def ufunc_per_section( 4. Obtain the coordinates of the measurements per section >>> locs = d.ufunc_per_section( + >>> sections=sections, >>> func=None, >>> label='x', >>> temp_err=False, @@ -4389,6 +4406,7 @@ def ufunc_per_section( 5. Number of observations per stretch >>> nlocs = d.ufunc_per_section( + >>> sections=sections, >>> func=len, >>> label='x', >>> temp_err=False, @@ -4407,7 +4425,7 @@ def ufunc_per_section( 7. x-coordinate index - >>> ix_loc = d.ufunc_per_section(x_indices=True) + >>> ix_loc = d.ufunc_per_section(sections=sections, x_indices=True) Note @@ -4415,9 +4433,8 @@ def ufunc_per_section( If `self[label]` or `self[subtract_from_label]` is a Dask array, a Dask array is returned else a numpy array is returned """ - if sections is None: - sections = self.sections - + # if sections is None: + # sections = self.sections if label is None: dataarray = None else: @@ -4426,7 +4443,10 @@ def ufunc_per_section( if x_indices: x_coords = self.x reference_dataset = None + else: + sections = validate_sections(self, sections) + x_coords = None reference_dataset = {k: self[k] for k in sections} diff --git a/src/dtscalibration/plot.py b/src/dtscalibration/plot.py index 27a747e4..43815b86 100755 --- a/src/dtscalibration/plot.py +++ b/src/dtscalibration/plot.py @@ -509,13 +509,14 @@ def plot_accuracy( def plot_sigma_report( - ds, temp_label, temp_var_acc_label, temp_var_prec_label=None, itimes=None + ds, sections, temp_label, temp_var_acc_label, temp_var_prec_label=None, itimes=None ): """Returns two sub-plots. first a temperature with confidence boundaries. Parameters ---------- ds + sections temp_label temp_var_label itimes @@ -573,11 +574,20 @@ def plot_sigma_report( # temp_err=True, # axis=0) sigma_est = ds.ufunc_per_section( - label=temp_label, func=np.std, temp_err=True, calc_per="stretch", axis=0 + sections=sections, + label=temp_label, + func=np.std, + temp_err=True, + calc_per="stretch", + axis=0, ) else: sigma_est = ds.ufunc_per_section( - label=temp_label, func=np.std, temp_err=True, calc_per="stretch" + sections=sections, + label=temp_label, + func=np.std, + temp_err=True, + calc_per="stretch", ) for (k, v), (k_se, v_se) in zip(ds.sections.items(), sigma_est.items()): @@ -628,9 +638,13 @@ def plot_sigma_report( ax1.set_ylabel(r"Temperature [$^\circ$C]") err_ref = ds.ufunc_per_section( - label=temp_label, func=None, temp_err=True, calc_per="stretch" + sections=sections, + label=temp_label, + func=None, + temp_err=True, + calc_per="stretch", ) - x_ref = ds.ufunc_per_section(label="x", calc_per="stretch") + x_ref = ds.ufunc_per_section(sections=sections, label="x", calc_per="stretch") for (k, v), (k_se, v_se), (kx, vx) in zip( ds.sections.items(), err_ref.items(), x_ref.items() diff --git a/tests/test_dtscalibration.py b/tests/test_dtscalibration.py index 805a900e..c4b674c1 100644 --- a/tests/test_dtscalibration.py +++ b/tests/test_dtscalibration.py @@ -321,6 +321,7 @@ def test_variance_input_types_double(): ) ds.conf_int_double_ended( + sections=sections, st_var=st_var, ast_var=st_var, rst_var=st_var, @@ -353,6 +354,7 @@ def st_var_callable(stokes): ) ds.conf_int_double_ended( + sections=sections, st_var=st_var_callable, ast_var=st_var_callable, rst_var=st_var_callable, @@ -382,6 +384,7 @@ def st_var_callable(stokes): ) ds.conf_int_double_ended( + sections=sections, st_var=st_var, ast_var=st_var, rst_var=st_var, @@ -408,6 +411,7 @@ def st_var_callable(stokes): ) ds.conf_int_double_ended( + sections=sections, st_var=st_var, ast_var=st_var, rst_var=st_var, @@ -437,6 +441,7 @@ def st_var_callable(stokes): ) ds.conf_int_double_ended( + sections=sections, st_var=st_var, ast_var=st_var, rst_var=st_var, @@ -565,6 +570,7 @@ def test_double_ended_variance_estimate_synthetic(): assert_almost_equal_verbose(ds.tmpb.mean(), 12.0, decimal=3) ds.conf_int_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=mst_var, @@ -577,20 +583,20 @@ def test_double_ended_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section( + stdsf1 = ds.ufunc_per_section(sections=sections, label="tmpf", func=np.std, temp_err=True, calc_per="stretch" ) - stdsb1 = ds.ufunc_per_section( + stdsb1 = ds.ufunc_per_section(sections=sections, label="tmpb", func=np.std, temp_err=True, calc_per="stretch" ) # Use a single timestep to better check if the parameter uncertainties propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section( + stdsf2 = ds1.ufunc_per_section(sections=sections, label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) - stdsb2 = ds1.ufunc_per_section( + stdsb2 = ds1.ufunc_per_section(sections=sections, label="tmpb_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) @@ -701,14 +707,14 @@ def test_single_ended_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section( + stdsf1 = ds.ufunc_per_section(sections=sections, label="tmpf", func=np.std, temp_err=True, calc_per="stretch", ddof=1 ) # Use a single timestep to better check if the parameter uncertainties propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section( + stdsf2 = ds1.ufunc_per_section(sections=sections, label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) @@ -2513,6 +2519,7 @@ def test_double_ended_exponential_variance_estimate_synthetic(): ) ds.conf_int_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_label=st_label, @@ -2529,10 +2536,10 @@ def test_double_ended_exponential_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section( + stdsf1 = ds.ufunc_per_section(sections=sections, label="tmpf", func=np.std, temp_err=True, calc_per="stretch" ) - stdsb1 = ds.ufunc_per_section( + stdsb1 = ds.ufunc_per_section(sections=sections, label="tmpb", func=np.std, temp_err=True, calc_per="stretch" ) @@ -2540,10 +2547,10 @@ def test_double_ended_exponential_variance_estimate_synthetic(): ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section( + stdsf2 = ds1.ufunc_per_section(sections=sections, label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) - stdsb2 = ds1.ufunc_per_section( + stdsb2 = ds1.ufunc_per_section(sections=sections, label="tmpb_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) @@ -2662,6 +2669,7 @@ def test_estimate_variance_of_temperature_estimate(): ) ds.conf_int_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=mst_var, @@ -3526,7 +3534,7 @@ def test_single_ended_exponential_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section( + stdsf1 = ds.ufunc_per_section(sections=sections, label="tmpf", func=np.var, temp_err=True, calc_per="stretch", ddof=1 ) @@ -3534,7 +3542,7 @@ def test_single_ended_exponential_variance_estimate_synthetic(): # propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section( + stdsf2 = ds1.ufunc_per_section(sections=sections, label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" ) @@ -3675,6 +3683,7 @@ def test_average_measurements_double_ended(): solver="sparse", ) ds.average_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=st_var, @@ -3688,6 +3697,7 @@ def test_average_measurements_double_ended(): ) ix = ds.get_section_indices(slice(6, 10)) ds.average_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=st_var, @@ -3704,6 +3714,7 @@ def test_average_measurements_double_ended(): np.datetime64("2018-03-28T00:41:12.084000000"), ) ds.average_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=st_var, @@ -3717,6 +3728,7 @@ def test_average_measurements_double_ended(): ci_avg_time_sel=sl, ) ds.average_double_ended( + sections=sections, p_val="p_val", p_cov="p_cov", st_var=st_var, From bb1b477851a02a358067513affdce9af63917c47 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 12:43:23 +0200 Subject: [PATCH 07/22] Format --- src/dtscalibration/calibrate_utils.py | 4 +- src/dtscalibration/datastore.py | 2 +- tests/test_dtscalibration.py | 82 +++++++++++++++++++-------- 3 files changed, 61 insertions(+), 27 deletions(-) diff --git a/src/dtscalibration/calibrate_utils.py b/src/dtscalibration/calibrate_utils.py index 5c343ac8..fd0117b0 100644 --- a/src/dtscalibration/calibrate_utils.py +++ b/src/dtscalibration/calibrate_utils.py @@ -893,7 +893,7 @@ def calibration_double_ended_helper( [fix_gamma[0]], out[0][: 2 * nt], E_all_exact, - out[0][2 * nt + n_E_in_cal:], + out[0][2 * nt + n_E_in_cal :], ) ) p_val[1 + 2 * nt + split["ix_from_cal_match_to_glob"]] = out[0][ @@ -993,7 +993,7 @@ def calibration_double_ended_helper( p0_est = np.concatenate( ( split["p0_est"][: 1 + 2 * nt], - split["p0_est"][1 + 2 * nt + n_E_in_cal:], + split["p0_est"][1 + 2 * nt + n_E_in_cal :], ) ) X_E1 = sp.csr_matrix(([], ([], [])), shape=(nt * nx_sec, self.x.size)) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 69eac39c..fa367dbf 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -4443,7 +4443,7 @@ def ufunc_per_section( if x_indices: x_coords = self.x reference_dataset = None - + else: sections = validate_sections(self, sections) diff --git a/tests/test_dtscalibration.py b/tests/test_dtscalibration.py index c4b674c1..66fd43c2 100644 --- a/tests/test_dtscalibration.py +++ b/tests/test_dtscalibration.py @@ -583,21 +583,29 @@ def test_double_ended_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section(sections=sections, - label="tmpf", func=np.std, temp_err=True, calc_per="stretch" + stdsf1 = ds.ufunc_per_section( + sections=sections, label="tmpf", func=np.std, temp_err=True, calc_per="stretch" ) - stdsb1 = ds.ufunc_per_section(sections=sections, - label="tmpb", func=np.std, temp_err=True, calc_per="stretch" + stdsb1 = ds.ufunc_per_section( + sections=sections, label="tmpb", func=np.std, temp_err=True, calc_per="stretch" ) # Use a single timestep to better check if the parameter uncertainties propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section(sections=sections, - label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsf2 = ds1.ufunc_per_section( + sections=sections, + label="tmpf_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) - stdsb2 = ds1.ufunc_per_section(sections=sections, - label="tmpb_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsb2 = ds1.ufunc_per_section( + sections=sections, + label="tmpb_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) for (_, v1), (_, v2) in zip(stdsf1.items(), stdsf2.items()): @@ -707,15 +715,24 @@ def test_single_ended_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section(sections=sections, - label="tmpf", func=np.std, temp_err=True, calc_per="stretch", ddof=1 + stdsf1 = ds.ufunc_per_section( + sections=sections, + label="tmpf", + func=np.std, + temp_err=True, + calc_per="stretch", + ddof=1, ) # Use a single timestep to better check if the parameter uncertainties propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section(sections=sections, - label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsf2 = ds1.ufunc_per_section( + sections=sections, + label="tmpf_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) for (_, v1), (_, v2) in zip(stdsf1.items(), stdsf2.items()): @@ -2536,22 +2553,30 @@ def test_double_ended_exponential_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section(sections=sections, - label="tmpf", func=np.std, temp_err=True, calc_per="stretch" + stdsf1 = ds.ufunc_per_section( + sections=sections, label="tmpf", func=np.std, temp_err=True, calc_per="stretch" ) - stdsb1 = ds.ufunc_per_section(sections=sections, - label="tmpb", func=np.std, temp_err=True, calc_per="stretch" + stdsb1 = ds.ufunc_per_section( + sections=sections, label="tmpb", func=np.std, temp_err=True, calc_per="stretch" ) # Use a single timestep to better check if the parameter uncertainties propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section(sections=sections, - label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsf2 = ds1.ufunc_per_section( + sections=sections, + label="tmpf_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) - stdsb2 = ds1.ufunc_per_section(sections=sections, - label="tmpb_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsb2 = ds1.ufunc_per_section( + sections=sections, + label="tmpb_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) for (_, v1), (_, v2) in zip(stdsf1.items(), stdsf2.items()): @@ -3534,16 +3559,25 @@ def test_single_ended_exponential_variance_estimate_synthetic(): ) # Calibrated variance - stdsf1 = ds.ufunc_per_section(sections=sections, - label="tmpf", func=np.var, temp_err=True, calc_per="stretch", ddof=1 + stdsf1 = ds.ufunc_per_section( + sections=sections, + label="tmpf", + func=np.var, + temp_err=True, + calc_per="stretch", + ddof=1, ) # Use a single timestep to better check if the parameter uncertainties # propagate ds1 = ds.isel(time=1) # Estimated VAR - stdsf2 = ds1.ufunc_per_section(sections=sections, - label="tmpf_mc_var", func=np.mean, temp_err=False, calc_per="stretch" + stdsf2 = ds1.ufunc_per_section( + sections=sections, + label="tmpf_mc_var", + func=np.mean, + temp_err=False, + calc_per="stretch", ) for (_, v1), (_, v2) in zip(stdsf1.items(), stdsf2.items()): From b453f5965ed54257c2254da10bd8ec2b17677359 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 15:48:15 +0200 Subject: [PATCH 08/22] Moving away from calibrate_double_ended() inplace self.update(out) is still in place at the end --- src/dtscalibration/datastore.py | 454 ++++++++++++++------------ src/dtscalibration/datastore_utils.py | 39 +++ 2 files changed, 293 insertions(+), 200 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index fa367dbf..1c69553d 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -21,6 +21,7 @@ from dtscalibration.datastore_utils import ParameterIndexSingleEnded from dtscalibration.datastore_utils import check_deprecated_kwargs from dtscalibration.datastore_utils import check_timestep_allclose +from dtscalibration.datastore_utils import get_params_from_pval from dtscalibration.datastore_utils import ufunc_per_section_helper from dtscalibration.io_utils import _dim_attrs @@ -2258,215 +2259,243 @@ def calibration_double_ended( assert p_var.size == ip.npar assert p_cov.shape == (ip.npar, ip.npar) - # save estimates and variances to datastore, skip covariances - self["gamma"] = (tuple(), p_val[ip.gamma].item()) - self["alpha"] = (("x",), p_val[ip.alpha]) - self["df"] = (("time",), p_val[ip.df]) - self["db"] = (("time",), p_val[ip.db]) - - if nta: - self["talpha_fw"] = ( - ("time", "trans_att"), - p_val[ip.taf].reshape((nt, nta), order="C"), - ) - self["talpha_bw"] = ( - ("time", "trans_att"), - p_val[ip.tab].reshape((nt, nta), order="C"), - ) - self["talpha_fw_var"] = ( - ("time", "trans_att"), - p_var[ip.taf].reshape((nt, nta), order="C"), - ) - self["talpha_bw_var"] = ( - ("time", "trans_att"), - p_var[ip.tab].reshape((nt, nta), order="C"), - ) - else: - self["talpha_fw"] = (("time", "trans_att"), np.zeros((nt, 0))) - self["talpha_bw"] = (("time", "trans_att"), np.zeros((nt, 0))) - self["talpha_fw_var"] = (("time", "trans_att"), np.zeros((nt, 0))) - self["talpha_bw_var"] = (("time", "trans_att"), np.zeros((nt, 0))) + coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} + params = get_params_from_pval(ip, p_val, coords) + param_covs = get_params_from_pval(ip, p_var, coords) + out = xr.Dataset( + { + "tmpf": params["gamma"] + / ( + np.log(self.st / self.ast) + + params["df"] + + params["alpha"] + + params["talpha_fw_full"] + ) + - 273.15, + "tmpb": params["gamma"] + / ( + np.log(self.rst / self.rast) + + params["db"] + - params["alpha"] + + params["talpha_bw_full"] + ) + - 273.15, + } + ) - self["talpha_fw_full"] = ( - ("x", "time"), - ip.get_taf_values( - pval=p_val, x=self.x.values, trans_att=self.trans_att.values, axis="" - ), + # extract covariances and ensure broadcastable to (nx, nt) + param_covs["gamma_df"] = (("time",), p_cov[np.ix_(ip.gamma, ip.df)][0]) + param_covs["gamma_db"] = (("time",), p_cov[np.ix_(ip.gamma, ip.db)][0]) + param_covs["gamma_alpha"] = (("x",), p_cov[np.ix_(ip.alpha, ip.gamma)][:, 0]) + param_covs["df_db"] = ( + ("time",), + p_cov[ip.df, ip.db], ) - self["talpha_bw_full"] = ( - ("x", "time"), - ip.get_tab_values( - pval=p_val, x=self.x.values, trans_att=self.trans_att.values, axis="" + param_covs["alpha_df"] = ( + ( + "x", + "time", ), + p_cov[np.ix_(ip.alpha, ip.df)], ) - - self["tmpf"] = ( - self.gamma - / ( - np.log(self.st / self.ast) - + self["df"] - + self.alpha - + self["talpha_fw_full"] - ) - - 273.15 - ) - - self["tmpb"] = ( - self.gamma - / ( - np.log(self.rst / self.rast) - + self["db"] - - self.alpha - + self["talpha_bw_full"] - ) - - 273.15 + param_covs["alpha_db"] = ( + ( + "x", + "time", + ), + p_cov[np.ix_(ip.alpha, ip.db)], ) - - self["gamma_var"] = (tuple(), p_var[ip.gamma].item()) - self["alpha_var"] = (("x",), p_var[ip.alpha]) - self["df_var"] = (("time",), p_var[ip.df]) - self["db_var"] = (("time",), p_var[ip.db]) - self["talpha_fw_full_var"] = ( - ("x", "time"), + param_covs["tafw_gamma"] = ( + ( + "x", + "time", + ), ip.get_taf_values( - pval=p_var, x=self.x.values, trans_att=self.trans_att.values, axis="" + pval=p_cov[ip.gamma], + x=self.x.values, + trans_att=self.trans_att.values, + axis="", ), ) - self["talpha_bw_full_var"] = ( - ("x", "time"), + param_covs["tabw_gamma"] = ( + ( + "x", + "time", + ), ip.get_tab_values( - pval=p_var, x=self.x.values, trans_att=self.trans_att.values, axis="" + pval=p_cov[ip.gamma], + x=self.x.values, + trans_att=self.trans_att.values, + axis="", ), ) - - # extract covariances and ensure broadcastable to (nx, nt) - sigma2_gamma_df = p_cov[np.ix_(ip.gamma, ip.df)] - sigma2_gamma_db = p_cov[np.ix_(ip.gamma, ip.db)] - sigma2_gamma_alpha = p_cov[np.ix_(ip.alpha, ip.gamma)] - sigma2_df_db = p_cov[ip.df, ip.db][None] # Shape is [0, nt] ? - sigma2_alpha_df = p_cov[np.ix_(ip.alpha, ip.df)] - sigma2_alpha_db = p_cov[np.ix_(ip.alpha, ip.db)] - sigma2_tafw_gamma = ip.get_taf_values( - pval=p_cov[ip.gamma], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ) - sigma2_tabw_gamma = ip.get_tab_values( - pval=p_cov[ip.gamma], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ) - sigma2_tafw_alpha = ip.get_taf_values( - pval=p_cov[ip.alpha], - x=self.x.values, - trans_att=self.trans_att.values, - axis="x", + param_covs["tafw_alpha"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.alpha], + x=self.x.values, + trans_att=self.trans_att.values, + axis="x", + ), ) - sigma2_tabw_alpha = ip.get_tab_values( - pval=p_cov[ip.alpha], - x=self.x.values, - trans_att=self.trans_att.values, - axis="x", + param_covs["tabw_alpha"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.alpha], + x=self.x.values, + trans_att=self.trans_att.values, + axis="x", + ), ) - sigma2_tafw_df = ip.get_taf_values( - pval=p_cov[ip.df], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", + param_covs["tafw_df"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.df], + x=self.x.values, + trans_att=self.trans_att.values, + axis="time", + ), ) - sigma2_tafw_db = ip.get_taf_values( - pval=p_cov[ip.db], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", + param_covs["tafw_db"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.db], + x=self.x.values, + trans_att=self.trans_att.values, + axis="time", + ), ) - sigma2_tabw_db = ip.get_tab_values( - pval=p_cov[ip.db], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", + param_covs["tabw_db"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.db], + x=self.x.values, + trans_att=self.trans_att.values, + axis="time", + ), ) - sigma2_tabw_df = ip.get_tab_values( - pval=p_cov[ip.df], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", + param_covs["tabw_df"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.df], + x=self.x.values, + trans_att=self.trans_att.values, + axis="time", + ), ) # sigma2_tafw_tabw - tmpf = self["tmpf"] + 273.15 - tmpb = self["tmpb"] + 273.15 + tmpf = out["tmpf"] + 273.15 + tmpb = out["tmpb"] + 273.15 deriv_dict = dict( - T_gamma_fw=tmpf / self.gamma, - T_st_fw=-(tmpf**2) / (self.gamma * self.st), - T_ast_fw=tmpf**2 / (self.gamma * self.ast), - T_df_fw=-(tmpf**2) / self.gamma, - T_alpha_fw=-(tmpf**2) / self.gamma, - T_ta_fw=-(tmpf**2) / self.gamma, - T_gamma_bw=tmpb / self.gamma, - T_rst_bw=-(tmpb**2) / (self.gamma * self.rst), - T_rast_bw=tmpb**2 / (self.gamma * self.rast), - T_db_bw=-(tmpb**2) / self.gamma, - T_alpha_bw=tmpb**2 / self.gamma, - T_ta_bw=-(tmpb**2) / self.gamma, + T_gamma_fw=tmpf / params["gamma"], + T_st_fw=-(tmpf**2) / (params["gamma"] * self.st), + T_ast_fw=tmpf**2 / (params["gamma"] * self.ast), + T_df_fw=-(tmpf**2) / params["gamma"], + T_alpha_fw=-(tmpf**2) / params["gamma"], + T_ta_fw=-(tmpf**2) / params["gamma"], + T_gamma_bw=tmpb / params["gamma"], + T_rst_bw=-(tmpb**2) / (params["gamma"] * self.rst), + T_rast_bw=tmpb**2 / (params["gamma"] * self.rast), + T_db_bw=-(tmpb**2) / params["gamma"], + T_alpha_bw=tmpb**2 / params["gamma"], + T_ta_bw=-(tmpb**2) / params["gamma"], ) deriv_ds = xr.Dataset(deriv_dict) - self["deriv"] = deriv_ds.to_array(dim="com2") + out["deriv"] = deriv_ds.to_array(dim="com2") var_fw_dict = dict( dT_dst=deriv_ds.T_st_fw**2 * parse_st_var(self, st_var, st_label="st"), dT_dast=deriv_ds.T_ast_fw**2 * parse_st_var(self, ast_var, st_label="ast"), - dT_gamma=deriv_ds.T_gamma_fw**2 * self["gamma_var"], - dT_ddf=deriv_ds.T_df_fw**2 * self["df_var"], - dT_dalpha=deriv_ds.T_alpha_fw**2 * self["alpha_var"], - dT_dta=deriv_ds.T_ta_fw**2 * self["talpha_fw_full_var"], - dgamma_ddf=(2 * deriv_ds.T_gamma_fw * deriv_ds.T_df_fw * sigma2_gamma_df), + dT_gamma=deriv_ds.T_gamma_fw**2 * param_covs["gamma"], + dT_ddf=deriv_ds.T_df_fw**2 * param_covs["df"], + dT_dalpha=deriv_ds.T_alpha_fw**2 * param_covs["alpha"], + dT_dta=deriv_ds.T_ta_fw**2 * param_covs["talpha_fw_full"], + dgamma_ddf=( + 2 * deriv_ds.T_gamma_fw * deriv_ds.T_df_fw * param_covs["gamma_df"] + ), dgamma_dalpha=( - 2 * deriv_ds.T_gamma_fw * deriv_ds.T_alpha_fw * sigma2_gamma_alpha + 2 + * deriv_ds.T_gamma_fw + * deriv_ds.T_alpha_fw + * param_covs["gamma_alpha"] + ), + dalpha_ddf=( + 2 * deriv_ds.T_alpha_fw * deriv_ds.T_df_fw * param_covs["alpha_df"] + ), + dta_dgamma=( + 2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * param_covs["tafw_gamma"] + ), + dta_ddf=(2 * deriv_ds.T_ta_fw * deriv_ds.T_df_fw * param_covs["tafw_df"]), + dta_dalpha=( + 2 * deriv_ds.T_ta_fw * deriv_ds.T_alpha_fw * param_covs["tafw_alpha"] ), - dalpha_ddf=(2 * deriv_ds.T_alpha_fw * deriv_ds.T_df_fw * sigma2_alpha_df), - dta_dgamma=(2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * sigma2_tafw_gamma), - dta_ddf=(2 * deriv_ds.T_ta_fw * deriv_ds.T_df_fw * sigma2_tafw_df), - dta_dalpha=(2 * deriv_ds.T_ta_fw * deriv_ds.T_alpha_fw * sigma2_tafw_alpha), ) var_bw_dict = dict( dT_drst=deriv_ds.T_rst_bw**2 * parse_st_var(self, rst_var, st_label="rst"), dT_drast=deriv_ds.T_rast_bw**2 * parse_st_var(self, rast_var, st_label="rast"), - dT_gamma=deriv_ds.T_gamma_bw**2 * self["gamma_var"], - dT_ddb=deriv_ds.T_db_bw**2 * self["db_var"], - dT_dalpha=deriv_ds.T_alpha_bw**2 * self["alpha_var"], - dT_dta=deriv_ds.T_ta_bw**2 * self["talpha_bw_full_var"], - dgamma_ddb=(2 * deriv_ds.T_gamma_bw * deriv_ds.T_db_bw * sigma2_gamma_db), + dT_gamma=deriv_ds.T_gamma_bw**2 * param_covs["gamma"], + dT_ddb=deriv_ds.T_db_bw**2 * param_covs["db"], + dT_dalpha=deriv_ds.T_alpha_bw**2 * param_covs["alpha"], + dT_dta=deriv_ds.T_ta_bw**2 * param_covs["talpha_bw_full"], + dgamma_ddb=( + 2 * deriv_ds.T_gamma_bw * deriv_ds.T_db_bw * param_covs["gamma_db"] + ), dgamma_dalpha=( - 2 * deriv_ds.T_gamma_bw * deriv_ds.T_alpha_bw * sigma2_gamma_alpha + 2 + * deriv_ds.T_gamma_bw + * deriv_ds.T_alpha_bw + * param_covs["gamma_alpha"] + ), + dalpha_ddb=( + 2 * deriv_ds.T_alpha_bw * deriv_ds.T_db_bw * param_covs["alpha_db"] + ), + dta_dgamma=( + 2 * deriv_ds.T_ta_bw * deriv_ds.T_gamma_bw * param_covs["tabw_gamma"] + ), + dta_ddb=(2 * deriv_ds.T_ta_bw * deriv_ds.T_db_bw * param_covs["tabw_db"]), + dta_dalpha=( + 2 * deriv_ds.T_ta_bw * deriv_ds.T_alpha_bw * param_covs["tabw_alpha"] ), - dalpha_ddb=(2 * deriv_ds.T_alpha_bw * deriv_ds.T_db_bw * sigma2_alpha_db), - dta_dgamma=(2 * deriv_ds.T_ta_bw * deriv_ds.T_gamma_bw * sigma2_tabw_gamma), - dta_ddb=(2 * deriv_ds.T_ta_bw * deriv_ds.T_db_bw * sigma2_tabw_db), - dta_dalpha=(2 * deriv_ds.T_ta_bw * deriv_ds.T_alpha_bw * sigma2_tabw_alpha), ) - self["var_fw_da"] = xr.Dataset(var_fw_dict).to_array(dim="comp_fw") - self["var_bw_da"] = xr.Dataset(var_bw_dict).to_array(dim="comp_bw") + out["var_fw_da"] = xr.Dataset(var_fw_dict).to_array(dim="comp_fw") + out["var_bw_da"] = xr.Dataset(var_bw_dict).to_array(dim="comp_bw") - self["tmpf_var"] = self["var_fw_da"].sum(dim="comp_fw") - self["tmpb_var"] = self["var_bw_da"].sum(dim="comp_bw") + out["tmpf_var"] = out["var_fw_da"].sum(dim="comp_fw") + out["tmpb_var"] = out["var_bw_da"].sum(dim="comp_bw") # First estimate of tmpw_var - self["tmpw_var" + "_approx"] = 1 / (1 / self["tmpf_var"] + 1 / self["tmpb_var"]) - self["tmpw"] = ( - (tmpf / self["tmpf_var"] + tmpb / self["tmpb_var"]) - * self["tmpw_var" + "_approx"] + out["tmpw_var" + "_approx"] = 1 / (1 / out["tmpf_var"] + 1 / out["tmpb_var"]) + out["tmpw"] = ( + (tmpf / out["tmpf_var"] + tmpb / out["tmpb_var"]) + * out["tmpw_var" + "_approx"] ) - 273.15 - weightsf = self["tmpw_var" + "_approx"] / self["tmpf_var"] - weightsb = self["tmpw_var" + "_approx"] / self["tmpb_var"] + weightsf = out["tmpw_var" + "_approx"] / out["tmpf_var"] + weightsb = out["tmpw_var" + "_approx"] / out["tmpb_var"] deriv_dict2 = dict( T_gamma_w=weightsf * deriv_dict["T_gamma_fw"] @@ -2493,61 +2522,85 @@ def calibration_double_ended( * parse_st_var(self, rst_var, st_label="rst"), dT_drast=deriv_ds2.T_rast_w**2 * parse_st_var(self, rast_var, st_label="rast"), - dT_gamma=deriv_ds2.T_gamma_w**2 * self["gamma_var"], - dT_ddf=deriv_ds2.T_df_w**2 * self["df_var"], - dT_ddb=deriv_ds2.T_db_w**2 * self["db_var"], - dT_dalpha=deriv_ds2.T_alpha_w**2 * self["alpha_var"], - dT_dtaf=deriv_ds2.T_taf_w**2 * self["talpha_fw_full_var"], - dT_dtab=deriv_ds2.T_tab_w**2 * self["talpha_bw_full_var"], - dgamma_ddf=2 * deriv_ds2.T_gamma_w * deriv_ds2.T_df_w * sigma2_gamma_df, - dgamma_ddb=2 * deriv_ds2.T_gamma_w * deriv_ds2.T_db_w * sigma2_gamma_db, + dT_gamma=deriv_ds2.T_gamma_w**2 * param_covs["gamma"], + dT_ddf=deriv_ds2.T_df_w**2 * param_covs["df"], + dT_ddb=deriv_ds2.T_db_w**2 * param_covs["db"], + dT_dalpha=deriv_ds2.T_alpha_w**2 * param_covs["alpha"], + dT_dtaf=deriv_ds2.T_taf_w**2 * param_covs["talpha_fw_full"], + dT_dtab=deriv_ds2.T_tab_w**2 * param_covs["talpha_bw_full"], + dgamma_ddf=2 + * deriv_ds2.T_gamma_w + * deriv_ds2.T_df_w + * param_covs["gamma_df"], + dgamma_ddb=2 + * deriv_ds2.T_gamma_w + * deriv_ds2.T_db_w + * param_covs["gamma_db"], dgamma_dalpha=2 * deriv_ds2.T_gamma_w * deriv_ds2.T_alpha_w - * sigma2_gamma_alpha, - dgamma_dtaf=2 * deriv_ds2.T_gamma_w * deriv_ds2.T_taf_w * sigma2_tafw_gamma, - dgamma_dtab=2 * deriv_ds2.T_gamma_w * deriv_ds2.T_tab_w * sigma2_tabw_gamma, - ddf_ddb=2 * deriv_ds2.T_df_w * deriv_ds2.T_db_w * sigma2_df_db, - ddf_dalpha=2 * deriv_ds2.T_df_w * deriv_ds2.T_alpha_w * sigma2_alpha_df, - ddf_dtaf=2 * deriv_ds2.T_df_w * deriv_ds2.T_taf_w * sigma2_tafw_df, - ddf_dtab=2 * deriv_ds2.T_df_w * deriv_ds2.T_tab_w * sigma2_tabw_df, - ddb_dalpha=2 * deriv_ds2.T_db_w * deriv_ds2.T_alpha_w * sigma2_alpha_db, - ddb_dtaf=2 * deriv_ds2.T_db_w * deriv_ds2.T_taf_w * sigma2_tafw_db, - ddb_dtab=2 * deriv_ds2.T_db_w * deriv_ds2.T_tab_w * sigma2_tabw_db, - # dtaf_dtab=2 * deriv_ds2.T_tab_w * deriv_ds2.T_tab_w * sigma2_tafw_tabw, + * param_covs["gamma_alpha"], + dgamma_dtaf=2 + * deriv_ds2.T_gamma_w + * deriv_ds2.T_taf_w + * param_covs["tafw_gamma"], + dgamma_dtab=2 + * deriv_ds2.T_gamma_w + * deriv_ds2.T_tab_w + * param_covs["tabw_gamma"], + ddf_ddb=2 * deriv_ds2.T_df_w * deriv_ds2.T_db_w * param_covs["df_db"], + ddf_dalpha=2 + * deriv_ds2.T_df_w + * deriv_ds2.T_alpha_w + * param_covs["alpha_df"], + ddf_dtaf=2 * deriv_ds2.T_df_w * deriv_ds2.T_taf_w * param_covs["tafw_df"], + ddf_dtab=2 * deriv_ds2.T_df_w * deriv_ds2.T_tab_w * param_covs["tabw_df"], + ddb_dalpha=2 + * deriv_ds2.T_db_w + * deriv_ds2.T_alpha_w + * param_covs["alpha_db"], + ddb_dtaf=2 * deriv_ds2.T_db_w * deriv_ds2.T_taf_w * param_covs["tafw_db"], + ddb_dtab=2 * deriv_ds2.T_db_w * deriv_ds2.T_tab_w * param_covs["tabw_db"], + # dtaf_dtab=2 * deriv_ds2.T_tab_w * deriv_ds2.T_tab_w * param_covs["tafw_tabw"], ) - self["var_w_da"] = xr.Dataset(var_w_dict).to_array(dim="comp_w") - self["tmpw_var"] = self["var_w_da"].sum(dim="comp_w") + out["var_w_da"] = xr.Dataset(var_w_dict).to_array(dim="comp_w") + out["tmpw_var"] = out["var_w_da"].sum(dim="comp_w") + # Compute uncertainty solely due to noise in Stokes signal, neglecting parameter uncenrtainty tmpf_var_excl_par = ( - self["var_fw_da"].sel(comp_fw=["dT_dst", "dT_dast"]).sum(dim="comp_fw") + out["var_fw_da"].sel(comp_fw=["dT_dst", "dT_dast"]).sum(dim="comp_fw") ) tmpb_var_excl_par = ( - self["var_bw_da"].sel(comp_bw=["dT_drst", "dT_drast"]).sum(dim="comp_bw") - ) - self["tmpw_var" + "_lower"] = 1 / ( - 1 / tmpf_var_excl_par + 1 / tmpb_var_excl_par + out["var_bw_da"].sel(comp_bw=["dT_drst", "dT_drast"]).sum(dim="comp_bw") ) + out["tmpw_var" + "_lower"] = 1 / (1 / tmpf_var_excl_par + 1 / tmpb_var_excl_par) - self["tmpf"].attrs.update(_dim_attrs[("tmpf",)]) - self["tmpb"].attrs.update(_dim_attrs[("tmpb",)]) - self["tmpw"].attrs.update(_dim_attrs[("tmpw",)]) - self["tmpf_var"].attrs.update(_dim_attrs[("tmpf_var",)]) - self["tmpb_var"].attrs.update(_dim_attrs[("tmpb_var",)]) - self["tmpw_var"].attrs.update(_dim_attrs[("tmpw_var",)]) - self["tmpw_var" + "_approx"].attrs.update(_dim_attrs[("tmpw_var_approx",)]) - self["tmpw_var" + "_lower"].attrs.update(_dim_attrs[("tmpw_var_lower",)]) + out["tmpf"].attrs.update(_dim_attrs[("tmpf",)]) + out["tmpb"].attrs.update(_dim_attrs[("tmpb",)]) + out["tmpw"].attrs.update(_dim_attrs[("tmpw",)]) + out["tmpf_var"].attrs.update(_dim_attrs[("tmpf_var",)]) + out["tmpb_var"].attrs.update(_dim_attrs[("tmpb_var",)]) + out["tmpw_var"].attrs.update(_dim_attrs[("tmpw_var",)]) + out["tmpw_var" + "_approx"].attrs.update(_dim_attrs[("tmpw_var_approx",)]) + out["tmpw_var" + "_lower"].attrs.update(_dim_attrs[("tmpw_var_lower",)]) drop_vars = [ k for k, v in self.items() if {"params1", "params2"}.intersection(v.dims) ] for k in drop_vars: + print(f"removing {k}") del self[k] - self["p_val"] = (("params1",), p_val) - self["p_cov"] = (("params1", "params2"), p_cov) + out["p_val"] = (("params1",), p_val) + out["p_cov"] = (("params1", "params2"), p_cov) + + out.update(params) + for key, da in param_covs.items(): + out[key + "_var"] = da + + self.update(out) pass def average_single_ended( @@ -2879,6 +2932,7 @@ def average_single_ended( ) # The new CI dimension is added as # firsaxis self["tmpf_mc_avg2"] = (("CI", x_dim2), qq) + # Clean up the garbage. All arrays with a Monte Carlo dimension. if mc_remove_set_flag: remove_mc_set = [ diff --git a/src/dtscalibration/datastore_utils.py b/src/dtscalibration/datastore_utils.py index dc2daad3..a4b9f875 100644 --- a/src/dtscalibration/datastore_utils.py +++ b/src/dtscalibration/datastore_utils.py @@ -528,6 +528,45 @@ def get_netcdf_encoding( return encoding +def get_params_from_pval(ip, p_val, coords): + assert len(p_val) == ip.npar, "Length of p_val is incorrect" + + params = xr.Dataset(coords=coords) + + # save estimates and variances to datastore, skip covariances + params["gamma"] = (tuple(), p_val[ip.gamma].item()) + params["alpha"] = (("x",), p_val[ip.alpha]) + params["df"] = (("time",), p_val[ip.df]) + params["db"] = (("time",), p_val[ip.db]) + + if ip.nta: + params["talpha_fw"] = ( + ("time", "trans_att"), + p_val[ip.taf].reshape((ip.nt, ip.nta), order="C"), + ) + params["talpha_bw"] = ( + ("time", "trans_att"), + p_val[ip.tab].reshape((ip.nt, ip.nta), order="C"), + ) + else: + params["talpha_fw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) + params["talpha_bw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) + + params["talpha_fw_full"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_val, x=params.x.values, trans_att=params.trans_att.values, axis="" + ), + ) + params["talpha_bw_full"] = ( + ("x", "time"), + ip.get_tab_values( + pval=p_val, x=params.x.values, trans_att=params.trans_att.values, axis="" + ), + ) + return params + + def merge_double_ended( ds_fw: "DataStore", ds_bw: "DataStore", From e884f945402ea54e5c360bf184477ccacc72d0f0 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 16:03:06 +0200 Subject: [PATCH 09/22] Refactored calibrate_double_ended covariance calculation --- src/dtscalibration/datastore.py | 127 +---------------- src/dtscalibration/datastore_utils.py | 197 +++++++++++++++++++++----- 2 files changed, 167 insertions(+), 157 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 1c69553d..b5af307e 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -2260,8 +2260,9 @@ def calibration_double_ended( assert p_cov.shape == (ip.npar, ip.npar) coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} - params = get_params_from_pval(ip, p_val, coords) - param_covs = get_params_from_pval(ip, p_var, coords) + params = get_params_from_pval(ip, coords, p_val=p_val) + param_covs = get_params_from_pval(ip, coords, p_val=p_var, p_cov=p_cov) + out = xr.Dataset( { "tmpf": params["gamma"] @@ -2283,126 +2284,6 @@ def calibration_double_ended( } ) - # extract covariances and ensure broadcastable to (nx, nt) - param_covs["gamma_df"] = (("time",), p_cov[np.ix_(ip.gamma, ip.df)][0]) - param_covs["gamma_db"] = (("time",), p_cov[np.ix_(ip.gamma, ip.db)][0]) - param_covs["gamma_alpha"] = (("x",), p_cov[np.ix_(ip.alpha, ip.gamma)][:, 0]) - param_covs["df_db"] = ( - ("time",), - p_cov[ip.df, ip.db], - ) - param_covs["alpha_df"] = ( - ( - "x", - "time", - ), - p_cov[np.ix_(ip.alpha, ip.df)], - ) - param_covs["alpha_db"] = ( - ( - "x", - "time", - ), - p_cov[np.ix_(ip.alpha, ip.db)], - ) - param_covs["tafw_gamma"] = ( - ( - "x", - "time", - ), - ip.get_taf_values( - pval=p_cov[ip.gamma], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ), - ) - param_covs["tabw_gamma"] = ( - ( - "x", - "time", - ), - ip.get_tab_values( - pval=p_cov[ip.gamma], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ), - ) - param_covs["tafw_alpha"] = ( - ( - "x", - "time", - ), - ip.get_taf_values( - pval=p_cov[ip.alpha], - x=self.x.values, - trans_att=self.trans_att.values, - axis="x", - ), - ) - param_covs["tabw_alpha"] = ( - ( - "x", - "time", - ), - ip.get_tab_values( - pval=p_cov[ip.alpha], - x=self.x.values, - trans_att=self.trans_att.values, - axis="x", - ), - ) - param_covs["tafw_df"] = ( - ( - "x", - "time", - ), - ip.get_taf_values( - pval=p_cov[ip.df], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", - ), - ) - param_covs["tafw_db"] = ( - ( - "x", - "time", - ), - ip.get_taf_values( - pval=p_cov[ip.db], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", - ), - ) - param_covs["tabw_db"] = ( - ( - "x", - "time", - ), - ip.get_tab_values( - pval=p_cov[ip.db], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", - ), - ) - param_covs["tabw_df"] = ( - ( - "x", - "time", - ), - ip.get_tab_values( - pval=p_cov[ip.df], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", - ), - ) - # sigma2_tafw_tabw - tmpf = out["tmpf"] + 273.15 tmpb = out["tmpb"] + 273.15 @@ -2597,7 +2478,7 @@ def calibration_double_ended( out.update(params) - for key, da in param_covs.items(): + for key, da in param_covs.data_vars.items(): out[key + "_var"] = da self.update(out) diff --git a/src/dtscalibration/datastore_utils.py b/src/dtscalibration/datastore_utils.py index a4b9f875..e12cba3a 100644 --- a/src/dtscalibration/datastore_utils.py +++ b/src/dtscalibration/datastore_utils.py @@ -528,42 +528,171 @@ def get_netcdf_encoding( return encoding -def get_params_from_pval(ip, p_val, coords): - assert len(p_val) == ip.npar, "Length of p_val is incorrect" - - params = xr.Dataset(coords=coords) - - # save estimates and variances to datastore, skip covariances - params["gamma"] = (tuple(), p_val[ip.gamma].item()) - params["alpha"] = (("x",), p_val[ip.alpha]) - params["df"] = (("time",), p_val[ip.df]) - params["db"] = (("time",), p_val[ip.db]) - - if ip.nta: - params["talpha_fw"] = ( - ("time", "trans_att"), - p_val[ip.taf].reshape((ip.nt, ip.nta), order="C"), +def get_params_from_pval(ip, coords, p_val=None, p_cov=None): + if p_val is not None: + assert len(p_val) == ip.npar, "Length of p_val is incorrect" + + params = xr.Dataset(coords=coords) + + # save estimates and variances to datastore, skip covariances + params["gamma"] = (tuple(), p_val[ip.gamma].item()) + params["alpha"] = (("x",), p_val[ip.alpha]) + params["df"] = (("time",), p_val[ip.df]) + params["db"] = (("time",), p_val[ip.db]) + + if ip.nta: + params["talpha_fw"] = ( + ("time", "trans_att"), + p_val[ip.taf].reshape((ip.nt, ip.nta), order="C"), + ) + params["talpha_bw"] = ( + ("time", "trans_att"), + p_val[ip.tab].reshape((ip.nt, ip.nta), order="C"), + ) + else: + params["talpha_fw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) + params["talpha_bw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) + + params["talpha_fw_full"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_val, + x=params.x.values, + trans_att=params.trans_att.values, + axis="", + ), ) - params["talpha_bw"] = ( - ("time", "trans_att"), - p_val[ip.tab].reshape((ip.nt, ip.nta), order="C"), + params["talpha_bw_full"] = ( + ("x", "time"), + ip.get_tab_values( + pval=p_val, + x=params.x.values, + trans_att=params.trans_att.values, + axis="", + ), ) - else: - params["talpha_fw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) - params["talpha_bw"] = (("time", "trans_att"), np.zeros((ip.nt, 0))) - - params["talpha_fw_full"] = ( - ("x", "time"), - ip.get_taf_values( - pval=p_val, x=params.x.values, trans_att=params.trans_att.values, axis="" - ), - ) - params["talpha_bw_full"] = ( - ("x", "time"), - ip.get_tab_values( - pval=p_val, x=params.x.values, trans_att=params.trans_att.values, axis="" - ), - ) + if p_cov is not None: + assert p_cov.shape == (ip.npar, ip.npar), "Shape of p_cov is incorrect" + + # extract covariances and ensure broadcastable to (nx, nt) + params["gamma_df"] = (("time",), p_cov[np.ix_(ip.gamma, ip.df)][0]) + params["gamma_db"] = (("time",), p_cov[np.ix_(ip.gamma, ip.db)][0]) + params["gamma_alpha"] = (("x",), p_cov[np.ix_(ip.alpha, ip.gamma)][:, 0]) + params["df_db"] = ( + ("time",), + p_cov[ip.df, ip.db], + ) + params["alpha_df"] = ( + ( + "x", + "time", + ), + p_cov[np.ix_(ip.alpha, ip.df)], + ) + params["alpha_db"] = ( + ( + "x", + "time", + ), + p_cov[np.ix_(ip.alpha, ip.db)], + ) + params["tafw_gamma"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.gamma], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ), + ) + params["tabw_gamma"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.gamma], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ), + ) + params["tafw_alpha"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.alpha], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="x", + ), + ) + params["tabw_alpha"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.alpha], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="x", + ), + ) + params["tafw_df"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.df], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ), + ) + params["tafw_db"] = ( + ( + "x", + "time", + ), + ip.get_taf_values( + pval=p_cov[ip.db], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ), + ) + params["tabw_db"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.db], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ), + ) + params["tabw_df"] = ( + ( + "x", + "time", + ), + ip.get_tab_values( + pval=p_cov[ip.df], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ), + ) + # sigma2_tafw_tabw return params From a4a1675edb4abb9580a4fa66f2ed4f8dd4ad1cda Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 19:38:55 +0200 Subject: [PATCH 10/22] Parallel testing --- pyproject.toml | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ccb756f8..a3387d12 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,8 +56,8 @@ dependencies = [ "dask", "pandas", "xarray[parallel]", # numbagg (llvmlite) is a pain to install with pip - "bottleneck", # speed up Xarray - "flox", # speed up Xarray + "bottleneck", # optional, speed up Xarray + "flox", # optional, speed up Xarray "pyyaml>=6.0.1", "xmltodict", "scipy", @@ -74,13 +74,14 @@ dev = [ "bump2version", "ruff", "isort", - "black[jupyter]", - "mypy", + "black[jupyter]", # for formatting + "mypy", # for type checking "types-PyYAML", # for pyyaml types "types-xmltodict", # for xmltodict types "pandas-stubs", # for pandas types "pytest", - "pytest-cov", + "pytest-cov", # for coverage + "pytest-xdist", # for parallel testing "jupyter", "nbformat", # Needed to run the tests ] @@ -112,8 +113,8 @@ lint = [ "mypy src/", ] format = ["black .", "isort .", "lint",] -test = ["pytest ./src/ ./tests/",] # --doctest-modules -fast-test = ["pytest ./tests/ -m \"not slow\"",] +test = ["pytest -n auto ./src/ ./tests/",] +fast-test = ["pytest -n auto ./tests/ -m \"not slow\"",] coverage = [ "pytest --cov --cov-report term --cov-report xml --junitxml=xunit-result.xml tests/", ] @@ -149,7 +150,7 @@ select = [ # It would be nice to have the commented out checks working. # "D", # pydocstyle # "C90", # mccabe complexity # "N", # PEP8-naming - # "UP", # pyupgrade (upgrade syntax to current syntax) + "UP", # pyupgrade (upgrade syntax to current syntax) "PLE", # Pylint error https://github.com/charliermarsh/ruff#error-ple # "PLR", # Pylint refactor (e.g. too-many-arguments) # "PLW", # Pylint warning (useless-else-on-loop) From f363413ba6479c2a59db8f9bb58b89f88d1d08e0 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 21:18:04 +0200 Subject: [PATCH 11/22] refactored variance functions --- pyproject.toml | 6 +- src/dtscalibration/__init__.py | 1 - src/dtscalibration/calibrate_utils.py | 1 - src/dtscalibration/datastore.py | 181 +++---------------------- src/dtscalibration/variance_helpers.py | 169 +++++++++++++++++++++++ 5 files changed, 193 insertions(+), 165 deletions(-) create mode 100644 src/dtscalibration/variance_helpers.py diff --git a/pyproject.toml b/pyproject.toml index a3387d12..959c779e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -113,8 +113,8 @@ lint = [ "mypy src/", ] format = ["black .", "isort .", "lint",] -test = ["pytest -n auto ./src/ ./tests/",] -fast-test = ["pytest -n auto ./tests/ -m \"not slow\"",] +test = ["pytest -n auto --dist worksteal ./src/ ./tests/",] +fast-test = ["pytest -n auto --dist worksteal ./tests/ -m \"not slow\"",] coverage = [ "pytest --cov --cov-report term --cov-report xml --junitxml=xunit-result.xml tests/", ] @@ -150,7 +150,7 @@ select = [ # It would be nice to have the commented out checks working. # "D", # pydocstyle # "C90", # mccabe complexity # "N", # PEP8-naming - "UP", # pyupgrade (upgrade syntax to current syntax) + # "UP", # pyupgrade (upgrade syntax to current syntax) "PLE", # Pylint error https://github.com/charliermarsh/ruff#error-ple # "PLR", # Pylint refactor (e.g. too-many-arguments) # "PLW", # Pylint warning (useless-else-on-loop) diff --git a/src/dtscalibration/__init__.py b/src/dtscalibration/__init__.py index 3836a779..4ecabb02 100644 --- a/src/dtscalibration/__init__.py +++ b/src/dtscalibration/__init__.py @@ -1,4 +1,3 @@ -# coding=utf-8 from dtscalibration.datastore import DataStore from dtscalibration.datastore_utils import check_dims from dtscalibration.datastore_utils import check_timestep_allclose diff --git a/src/dtscalibration/calibrate_utils.py b/src/dtscalibration/calibrate_utils.py index fd0117b0..960fa220 100644 --- a/src/dtscalibration/calibrate_utils.py +++ b/src/dtscalibration/calibrate_utils.py @@ -1,4 +1,3 @@ -# coding=utf-8 import numpy as np import scipy.sparse as sp import statsmodels.api as sm diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index b5af307e..6bb774f8 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -4,12 +4,9 @@ import dask import dask.array as da import numpy as np -import scipy.sparse as sp import scipy.stats as sst import xarray as xr import yaml -from scipy.optimize import minimize -from scipy.sparse import linalg as ln from dtscalibration.calibrate_utils import calibration_double_ended_helper from dtscalibration.calibrate_utils import calibration_single_ended_helper @@ -24,6 +21,9 @@ from dtscalibration.datastore_utils import get_params_from_pval from dtscalibration.datastore_utils import ufunc_per_section_helper from dtscalibration.io_utils import _dim_attrs +from dtscalibration.variance_helpers import variance_stokes_constant_helper +from dtscalibration.variance_helpers import variance_stokes_exponential_helper +from dtscalibration.variance_helpers import variance_stokes_linear_helper dtsattr_namelist = ["double_ended_flag"] dim_attrs = {k: v for kl, v in _dim_attrs.items() for k in kl} @@ -879,21 +879,12 @@ def variance_stokes_constant(self, st_label, sections=None, reshape_residuals=Tr dtscalibration/python-dts-calibration/blob/main/examples/notebooks/\ 04Calculate_variance_Stokes.ipynb>`_ """ - - def func_fit(p, xs): - return p[:xs, None] * p[None, xs:] - - def func_cost(p, data, xs): - fit = func_fit(p, xs) - return np.sum((fit - data) ** 2) - if sections is None: sections = self.sections else: sections = validate_sections(self, sections) - assert self[st_label].dims[0] == "x", "Stokes are transposed" - + assert self[st_label].dims[0] == "x", f"{st_label} are transposed" check_timestep_allclose(self, eps=0.01) # should maybe be per section. But then residuals @@ -904,25 +895,7 @@ def func_cost(p, data, xs): ) )[0] - resid_list = [] - - for k, v in data_dict.items(): - for vi in v: - nxs, nt = vi.shape - npar = nt + nxs - - p1 = np.ones(npar) * vi.mean() ** 0.5 - - res = minimize(func_cost, p1, args=(vi, nxs), method="Powell") - assert res.success, "Unable to fit. Try variance_stokes_exponential" - - fit = func_fit(res.x, nxs) - resid_list.append(fit - vi) - - resid = np.concatenate(resid_list) - - # unbiased estimater ddof=1, originally thought it was npar - var_I = resid.var(ddof=1) + var_I, resid = variance_stokes_constant_helper(data_dict) if not reshape_residuals: return var_I, resid @@ -1096,90 +1069,13 @@ def variance_stokes_exponential( x_list.append(da.tile(_x, nt)) len_stretch_list.append(_x.size) - n_sections = len(len_stretch_list) # number of sections - n_locs = sum(len_stretch_list) # total number of locations along cable used - # for reference. - x = np.concatenate(x_list) # coordinates are already in memory y = np.concatenate(y_list) - data1 = x - data2 = np.ones(sum(len_stretch_list) * nt) - data = np.concatenate([data1, data2]) - - # alpha is NOT the same for all -> one column per section - coords1row = np.arange(nt * n_locs) - coords1col = np.hstack( - [np.ones(in_locs * nt) * i for i, in_locs in enumerate(len_stretch_list)] - ) # C for - - # second calibration parameter is different per section and per timestep - coords2row = np.arange(nt * n_locs) - coords2col = np.hstack( - [ - np.repeat( - np.arange(i * nt + n_sections, (i + 1) * nt + n_sections), in_locs - ) - for i, in_locs in enumerate(len_stretch_list) - ] - ) # C for - coords = ( - np.concatenate([coords1row, coords2row]), - np.concatenate([coords1col, coords2col]), + var_I, resid = variance_stokes_exponential_helper( + nt, x, y, len_stretch_list, use_statsmodels, suppress_info ) - lny = np.log(y) - w = y.copy() # 1/std. - - ddof = n_sections + nt * n_sections # see numpy documentation on ddof - - if use_statsmodels: - # returns the same answer with statsmodel - import statsmodels.api as sm - - X = sp.coo_matrix( - (data, coords), shape=(nt * n_locs, ddof), dtype=float, copy=False - ) - - mod_wls = sm.WLS(lny, X.toarray(), weights=w**2) - res_wls = mod_wls.fit() - # print(res_wls.summary()) - a = res_wls.params - - else: - wdata = data * np.hstack((w, w)) - wX = sp.coo_matrix( - (wdata, coords), - shape=(nt * n_locs, n_sections + nt * n_sections), - dtype=float, - copy=False, - ) - - wlny = lny * w - - p0_est = np.asarray(n_sections * [0.0] + nt * n_sections * [8]) - # noinspection PyTypeChecker - a = ln.lsqr(wX, wlny, x0=p0_est, show=not suppress_info, calc_var=False)[0] - - beta = a[:n_sections] - beta_expand_to_sec = np.hstack( - [ - np.repeat(float(beta[i]), leni * nt) - for i, leni in enumerate(len_stretch_list) - ] - ) - G = np.asarray(a[n_sections:]) - G_expand_to_sec = np.hstack( - [ - np.repeat(G[i * nt : (i + 1) * nt], leni) - for i, leni in enumerate(len_stretch_list) - ] - ) - - I_est = np.exp(G_expand_to_sec) * np.exp(x * beta_expand_to_sec) - resid = I_est - y - var_I = resid.var(ddof=1) - if not reshape_residuals: return var_I, resid @@ -1336,60 +1232,25 @@ def variance_stokes_linear( sections = validate_sections(self, sections) assert self[st_label].dims[0] == "x", "Stokes are transposed" - _, resid = self.variance_stokes(sections=sections, st_label=st_label) + _, resid = self.variance_stokes( + sections=sections, st_label=st_label, reshape_residuals=False + ) ix_sec = self.ufunc_per_section( sections=sections, x_indices=True, calc_per="all" ) - st = self.isel(x=ix_sec)[st_label].values.ravel() - diff_st = resid.isel(x=ix_sec).values.ravel() - - # Adjust nbin silently to fit residuals in - # rectangular matrix and use numpy for computation - nbin_ = nbin - while st.size % nbin_: - nbin_ -= 1 - - if nbin_ != nbin: - print( - "Estimation of linear variance of", - st_label, - "Adjusting nbin to:", - nbin_, - ) - nbin = nbin_ - - isort = np.argsort(st) - st_sort_mean = st[isort].reshape((nbin, -1)).mean(axis=1) - st_sort_var = diff_st[isort].reshape((nbin, -1)).var(axis=1) - - if through_zero: - # VAR(Stokes) = slope * Stokes - offset = 0.0 - slope = np.linalg.lstsq(st_sort_mean[:, None], st_sort_var, rcond=None)[0] - - else: - # VAR(Stokes) = slope * Stokes + offset - slope, offset = np.linalg.lstsq( - np.hstack((st_sort_mean[:, None], np.ones((nbin, 1)))), - st_sort_var, - rcond=None, - )[0] - - if offset < 0: - warnings.warn( - f"Warning! Offset of variance_stokes_linear() " - f"of {st_label} is negative. This is phisically " - f"not possible. Most likely, your {st_label} do " - f"not vary enough to fit a linear curve. Either " - f"use `through_zero` option or use " - f"`ds.variance_stokes_constant()`. Another reason " - f"could be that your sections are defined to be " - f"wider than they actually are." - ) - def var_fun(stokes): - return slope * stokes + offset + st = self[st_label].isel(x=ix_sec).values.ravel() + diff_st = resid.ravel() + + ( + slope, + offset, + st_sort_mean, + st_sort_var, + resid, + var_fun, + ) = variance_stokes_linear_helper(st, diff_st, nbin, through_zero) if plot_fit: plt.figure() diff --git a/src/dtscalibration/variance_helpers.py b/src/dtscalibration/variance_helpers.py new file mode 100644 index 00000000..e1b930c3 --- /dev/null +++ b/src/dtscalibration/variance_helpers.py @@ -0,0 +1,169 @@ +import warnings + +import numpy as np +import scipy.sparse as sp +from scipy.optimize import minimize +from scipy.sparse import linalg as ln + + +def variance_stokes_constant_helper(data_dict): + def func_fit(p, xs): + return p[:xs, None] * p[None, xs:] + + def func_cost(p, data, xs): + fit = func_fit(p, xs) + return np.sum((fit - data) ** 2) + + resid_list = [] + + for k, v in data_dict.items(): + for vi in v: + nxs, nt = vi.shape + npar = nt + nxs + + p1 = np.ones(npar) * vi.mean() ** 0.5 + + res = minimize(func_cost, p1, args=(vi, nxs), method="Powell") + assert res.success, "Unable to fit. Try variance_stokes_exponential" + + fit = func_fit(res.x, nxs) + resid_list.append(fit - vi) + + resid = np.concatenate(resid_list) + + # unbiased estimater ddof=1, originally thought it was npar + var_I = resid.var(ddof=1) + + return var_I, resid + + +def variance_stokes_exponential_helper( + nt, x, y, len_stretch_list, use_statsmodels, suppress_info +): + n_sections = len(len_stretch_list) # number of sections + n_locs = sum(len_stretch_list) # total number of locations along cable used + # for reference. + + data1 = x + data2 = np.ones(sum(len_stretch_list) * nt) + data = np.concatenate([data1, data2]) + + # alpha is NOT the same for all -> one column per section + coords1row = np.arange(nt * n_locs) + coords1col = np.hstack( + [np.ones(in_locs * nt) * i for i, in_locs in enumerate(len_stretch_list)] + ) # C for + + # second calibration parameter is different per section and per timestep + coords2row = np.arange(nt * n_locs) + coords2col = np.hstack( + [ + np.repeat( + np.arange(i * nt + n_sections, (i + 1) * nt + n_sections), in_locs + ) + for i, in_locs in enumerate(len_stretch_list) + ] + ) # C for + coords = ( + np.concatenate([coords1row, coords2row]), + np.concatenate([coords1col, coords2col]), + ) + + lny = np.log(y) + w = y.copy() # 1/std. + + ddof = n_sections + nt * n_sections # see numpy documentation on ddof + + if use_statsmodels: + # returns the same answer with statsmodel + import statsmodels.api as sm + + X = sp.coo_matrix( + (data, coords), shape=(nt * n_locs, ddof), dtype=float, copy=False + ) + + mod_wls = sm.WLS(lny, X.toarray(), weights=w**2) + res_wls = mod_wls.fit() + # print(res_wls.summary()) + a = res_wls.params + + else: + wdata = data * np.hstack((w, w)) + wX = sp.coo_matrix( + (wdata, coords), + shape=(nt * n_locs, n_sections + nt * n_sections), + dtype=float, + copy=False, + ) + + wlny = lny * w + + p0_est = np.asarray(n_sections * [0.0] + nt * n_sections * [8]) + # noinspection PyTypeChecker + a = ln.lsqr(wX, wlny, x0=p0_est, show=not suppress_info, calc_var=False)[0] + + beta = a[:n_sections] + beta_expand_to_sec = np.hstack( + [ + np.repeat(float(beta[i]), leni * nt) + for i, leni in enumerate(len_stretch_list) + ] + ) + G = np.asarray(a[n_sections:]) + G_expand_to_sec = np.hstack( + [ + np.repeat(G[i * nt : (i + 1) * nt], leni) + for i, leni in enumerate(len_stretch_list) + ] + ) + + I_est = np.exp(G_expand_to_sec) * np.exp(x * beta_expand_to_sec) + resid = I_est - y + var_I = resid.var(ddof=1) + return var_I, resid + + +def variance_stokes_linear_helper(st_sec, resid_sec, nbin, through_zero): + # Adjust nbin silently to fit residuals in + # rectangular matrix and use numpy for computation + nbin_ = nbin + while st_sec.size % nbin_: + nbin_ -= 1 + + if nbin_ != nbin: + print(f"Adjusting nbin to: {nbin_} to fit residuals in ") + nbin = nbin_ + + isort = np.argsort(st_sec) + st_sort_mean = st_sec[isort].reshape((nbin, -1)).mean(axis=1) + st_sort_var = resid_sec[isort].reshape((nbin, -1)).var(axis=1) + + if through_zero: + # VAR(Stokes) = slope * Stokes + offset = 0.0 + slope = np.linalg.lstsq(st_sort_mean[:, None], st_sort_var, rcond=None)[0] + + else: + # VAR(Stokes) = slope * Stokes + offset + slope, offset = np.linalg.lstsq( + np.hstack((st_sort_mean[:, None], np.ones((nbin, 1)))), + st_sort_var, + rcond=None, + )[0] + + if offset < 0: + warnings.warn( + "Warning! Offset of variance_stokes_linear() " + "is negative. This is phisically " + "not possible. Most likely, your Stokes intensities do " + "not vary enough to fit a linear curve. Either " + "use `through_zero` option or use " + "`ds.variance_stokes_constant()`. Another reason " + "could be that your sections are defined to be " + "wider than they actually are." + ) + + def var_fun(stokes): + return slope * stokes + offset + + return slope, offset, st_sort_mean, st_sort_var, resid_sec, var_fun From 7f795b92698fa75ae8dedb0c7eeb2657207b0a87 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Fri, 18 Aug 2023 21:26:46 +0200 Subject: [PATCH 12/22] Include pyupgrade flag of ruff --- pyproject.toml | 4 ++-- src/dtscalibration/calibration/section_utils.py | 6 ++---- src/dtscalibration/datastore.py | 10 +++++----- src/dtscalibration/datastore_utils.py | 1 - src/dtscalibration/io.py | 6 +++--- src/dtscalibration/io_utils.py | 11 +++++------ src/dtscalibration/plot.py | 7 +++---- tests/test_datastore.py | 1 - tests/test_dtscalibration.py | 1 - 9 files changed, 20 insertions(+), 27 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 959c779e..3597d280 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -150,7 +150,7 @@ select = [ # It would be nice to have the commented out checks working. # "D", # pydocstyle # "C90", # mccabe complexity # "N", # PEP8-naming - # "UP", # pyupgrade (upgrade syntax to current syntax) + "UP", # pyupgrade (upgrade syntax to current syntax) "PLE", # Pylint error https://github.com/charliermarsh/ruff#error-ple # "PLR", # Pylint refactor (e.g. too-many-arguments) # "PLW", # Pylint warning (useless-else-on-loop) @@ -166,7 +166,7 @@ ignore = [ "E501", # Line too long (want to have fixed ] # Allow autofix for all enabled rules (when `--fix`) is provided. -fixable = ["E", "F"] +fixable = ["E", "F", "UP", "PLE"] unfixable = [] line-length = 88 exclude = ["docs", "build"] diff --git a/src/dtscalibration/calibration/section_utils.py b/src/dtscalibration/calibration/section_utils.py index 63dd1e5a..648642be 100644 --- a/src/dtscalibration/calibration/section_utils.py +++ b/src/dtscalibration/calibration/section_utils.py @@ -1,5 +1,3 @@ -from typing import Dict -from typing import List import numpy as np import xarray as xr @@ -8,7 +6,7 @@ from dtscalibration.datastore_utils import ufunc_per_section_helper -def set_sections(ds: xr.Dataset, sections: Dict[str, List[slice]]) -> xr.Dataset: +def set_sections(ds: xr.Dataset, sections: dict[str, list[slice]]) -> xr.Dataset: sections_validated = None if sections is not None: @@ -18,7 +16,7 @@ def set_sections(ds: xr.Dataset, sections: Dict[str, List[slice]]) -> xr.Dataset return ds -def validate_sections(ds: xr.Dataset, sections: Dict[str, List[slice]]): +def validate_sections(ds: xr.Dataset, sections: dict[str, list[slice]]): assert isinstance(sections, dict) # be less restrictive for capitalized labels diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 6bb774f8..64f7bd35 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -144,11 +144,11 @@ def __repr__(self): unit = "" for k, v in self.sections.items(): - preamble_new += " {0: <23}".format(k) + preamble_new += f" {k: <23}" # Compute statistics reference section timeseries - sec_stat = "({0:6.2f}".format(float(self[k].mean())) - sec_stat += " +/-{0:5.2f}".format(float(self[k].std())) + sec_stat = f"({float(self[k].mean()):6.2f}" + sec_stat += f" +/-{float(self[k].std()):5.2f}" sec_stat += "\N{DEGREE SIGN}C)\t" preamble_new += sec_stat @@ -425,7 +425,7 @@ def to_netcdf( if value is None: self.attrs[attribute] = "" - return super(DataStore, self).to_netcdf( + return super().to_netcdf( path, mode, format=format, @@ -552,7 +552,7 @@ def to_mf_netcdf( paths = [ os.path.join( folder_path, - filename_preamble + "{:04d}".format(ix) + filename_extension, + filename_preamble + f"{ix:04d}" + filename_extension, ) for ix in range(len(x)) ] diff --git a/src/dtscalibration/datastore_utils.py b/src/dtscalibration/datastore_utils.py index e12cba3a..339a45a7 100644 --- a/src/dtscalibration/datastore_utils.py +++ b/src/dtscalibration/datastore_utils.py @@ -1,4 +1,3 @@ -# coding=utf-8 from typing import TYPE_CHECKING from typing import Optional from typing import Union diff --git a/src/dtscalibration/io.py b/src/dtscalibration/io.py index be7d0dd6..edbeac63 100644 --- a/src/dtscalibration/io.py +++ b/src/dtscalibration/io.py @@ -282,7 +282,7 @@ def read_silixa_files( else: raise NotImplementedError( - "Silixa xml version " + "{0} not implemented".format(xml_version) + "Silixa xml version " + f"{xml_version} not implemented" ) ds = DataStore(data_vars=data_vars, coords=coords, attrs=attrs, **kwargs) @@ -332,7 +332,7 @@ def read_sensortran_files( # Check if corresponding temperature file exists if not os.path.isfile(filepathlist_temp[ii]): raise FileNotFoundError( - "Could not find BinaryTemp " + "file corresponding to {}".format(fname) + "Could not find BinaryTemp " + f"file corresponding to {fname}" ) version = sensortran_binary_version_check(filepathlist_dts) @@ -347,7 +347,7 @@ def read_sensortran_files( ) else: raise NotImplementedError( - "Sensortran binary version " + "{0} not implemented".format(version) + "Sensortran binary version " + f"{version} not implemented" ) ds = DataStore(data_vars=data_vars, coords=coords, attrs=attrs, **kwargs) diff --git a/src/dtscalibration/io_utils.py b/src/dtscalibration/io_utils.py index e880a6ca..4b480742 100644 --- a/src/dtscalibration/io_utils.py +++ b/src/dtscalibration/io_utils.py @@ -1,4 +1,3 @@ -# coding=utf-8 import os import struct from contextlib import contextmanager @@ -1263,8 +1262,8 @@ def read_sensortran_files_routine( AST = np.zeros((nx, ntime), dtype=np.int32) TMP = np.zeros((nx, ntime)) - ST_zero = np.zeros((ntime)) - AST_zero = np.zeros((ntime)) + ST_zero = np.zeros(ntime) + AST_zero = np.zeros(ntime) for ii in range(ntime): data_dts, meta_dts = read_sensortran_single(filepathlist_dts[ii]) @@ -1468,7 +1467,7 @@ def read_apsensing_files_routine( + "/{0}wellLogSet/{0}wellLog" ).format(namespace) ) - logdata_tree = logtree.find("./{0}logData".format(namespace)) + logdata_tree = logtree.find(f"./{namespace}logData") # Amount of datapoints is the size of the logdata tree nx = len(logdata_tree) @@ -1490,7 +1489,7 @@ def read_apsensing_files_routine( attrs["backwardMeasurementChannel"] = "N/A" data_item_names = [ - attrs["wellbore:wellLogSet:wellLog:logCurveInfo_{0}:mnemonic".format(x)] + attrs[f"wellbore:wellLogSet:wellLog:logCurveInfo_{x}:mnemonic"] for x in range(0, 4) ] @@ -1591,7 +1590,7 @@ def grab_data_per_file(file_handle): data_vars[tld[name]] = (["x", "time"], data_arri, dim_attrs_apsensing[name]) else: raise ValueError( - "Dont know what to do with the" + " {} data column".format(name) + "Dont know what to do with the" + f" {name} data column" ) _time_dtype = [("filename_tstamp", np.int64), ("acquisitionTime", " Date: Sun, 20 Aug 2023 18:07:33 +0200 Subject: [PATCH 13/22] Refactored single ended routine p_val to params and moved averaging routines --- src/dtscalibration/averaging_utils.py | 75 +++++++++++ src/dtscalibration/datastore.py | 173 +++----------------------- src/dtscalibration/datastore_utils.py | 76 ++++++++++- 3 files changed, 169 insertions(+), 155 deletions(-) create mode 100644 src/dtscalibration/averaging_utils.py diff --git a/src/dtscalibration/averaging_utils.py b/src/dtscalibration/averaging_utils.py new file mode 100644 index 00000000..9169c36e --- /dev/null +++ b/src/dtscalibration/averaging_utils.py @@ -0,0 +1,75 @@ + + +def inverse_variance_weighted_mean( + self, + tmp1="tmpf", + tmp2="tmpb", + tmp1_var="tmpf_mc_var", + tmp2_var="tmpb_mc_var", + tmpw_store="tmpw", + tmpw_var_store="tmpw_var", +): + """ + Average two temperature datasets with the inverse of the variance as + weights. The two + temperature datasets `tmp1` and `tmp2` with their variances + `tmp1_var` and `tmp2_var`, + respectively. Are averaged and stored in the DataStore. + + Parameters + ---------- + tmp1 : str + The label of the first temperature dataset that is averaged + tmp2 : str + The label of the second temperature dataset that is averaged + tmp1_var : str + The variance of tmp1 + tmp2_var : str + The variance of tmp2 + tmpw_store : str + The label of the averaged temperature dataset + tmpw_var_store : str + The label of the variance of the averaged temperature dataset + + Returns + ------- + + """ + + self[tmpw_var_store] = 1 / (1 / self[tmp1_var] + 1 / self[tmp2_var]) + + self[tmpw_store] = ( + self[tmp1] / self[tmp1_var] + self[tmp2] / self[tmp2_var] + ) * self[tmpw_var_store] + + pass + +def inverse_variance_weighted_mean_array( + self, + tmp_label="tmpf", + tmp_var_label="tmpf_mc_var", + tmpw_store="tmpw", + tmpw_var_store="tmpw_var", + dim="time", +): + """ + Calculates the weighted average across a dimension. + + Parameters + ---------- + + Returns + ------- + + See Also + -------- + - https://en.wikipedia.org/wiki/Inverse-variance_weighting + + """ + self[tmpw_var_store] = 1 / (1 / self[tmp_var_label]).sum(dim=dim) + + self[tmpw_store] = (self[tmp_label] / self[tmp_var_label]).sum(dim=dim) / ( + 1 / self[tmp_var_label] + ).sum(dim=dim) + + pass \ No newline at end of file diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 64f7bd35..3b6e9b49 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -18,7 +18,8 @@ from dtscalibration.datastore_utils import ParameterIndexSingleEnded from dtscalibration.datastore_utils import check_deprecated_kwargs from dtscalibration.datastore_utils import check_timestep_allclose -from dtscalibration.datastore_utils import get_params_from_pval +from dtscalibration.datastore_utils import get_params_from_pval_double_ended +from dtscalibration.datastore_utils import get_params_from_pval_single_ended from dtscalibration.datastore_utils import ufunc_per_section_helper from dtscalibration.io_utils import _dim_attrs from dtscalibration.variance_helpers import variance_stokes_constant_helper @@ -1331,80 +1332,6 @@ def i_var(self, st_var, ast_var, st_label="st", ast_label="ast"): return st**-2 * st_var + ast**-2 * ast_var - def inverse_variance_weighted_mean( - self, - tmp1="tmpf", - tmp2="tmpb", - tmp1_var="tmpf_mc_var", - tmp2_var="tmpb_mc_var", - tmpw_store="tmpw", - tmpw_var_store="tmpw_var", - ): - """ - Average two temperature datasets with the inverse of the variance as - weights. The two - temperature datasets `tmp1` and `tmp2` with their variances - `tmp1_var` and `tmp2_var`, - respectively. Are averaged and stored in the DataStore. - - Parameters - ---------- - tmp1 : str - The label of the first temperature dataset that is averaged - tmp2 : str - The label of the second temperature dataset that is averaged - tmp1_var : str - The variance of tmp1 - tmp2_var : str - The variance of tmp2 - tmpw_store : str - The label of the averaged temperature dataset - tmpw_var_store : str - The label of the variance of the averaged temperature dataset - - Returns - ------- - - """ - - self[tmpw_var_store] = 1 / (1 / self[tmp1_var] + 1 / self[tmp2_var]) - - self[tmpw_store] = ( - self[tmp1] / self[tmp1_var] + self[tmp2] / self[tmp2_var] - ) * self[tmpw_var_store] - - pass - - def inverse_variance_weighted_mean_array( - self, - tmp_label="tmpf", - tmp_var_label="tmpf_mc_var", - tmpw_store="tmpw", - tmpw_var_store="tmpw_var", - dim="time", - ): - """ - Calculates the weighted average across a dimension. - - Parameters - ---------- - - Returns - ------- - - See Also - -------- - - https://en.wikipedia.org/wiki/Inverse-variance_weighting - - """ - self[tmpw_var_store] = 1 / (1 / self[tmp_var_label]).sum(dim=dim) - - self[tmpw_store] = (self[tmp_label] / self[tmp_var_label]).sum(dim=dim) / ( - 1 / self[tmp_var_label] - ).sum(dim=dim) - - pass - def set_trans_att(self, trans_att=None): """Gracefully set the locations that introduce directional differential attenuation @@ -1669,52 +1596,12 @@ def calibration_single_ended( assert p_cov.shape == (ip.npar, ip.npar) # store calibration parameters in DataStore - self["gamma"] = (tuple(), p_val[ip.gamma].item()) - self["gamma_var"] = (tuple(), p_var[ip.gamma].item()) - - if nta > 0: - self["talpha_fw"] = ( - ("trans_att", "time"), - ip.get_taf_pars(p_val), - ) - self["talpha_fw_var"] = ( - ("trans_att", "time"), - ip.get_taf_pars(p_var), - ) - else: - self["talpha_fw"] = (("trans_att", "time"), np.zeros((0, nt))) - self["talpha_fw_var"] = (("trans_att", "time"), np.zeros((0, nt))) - - self["c"] = (("time",), p_val[ip.c]) - self["c_var"] = (("time",), p_var[ip.c]) - - if fix_alpha: - self["alpha"] = (("x",), fix_alpha[0]) - self["alpha_var"] = (("x",), fix_alpha[1]) - - else: - self["dalpha"] = (tuple(), p_val[ip.dalpha].item()) - self["dalpha_var"] = (tuple(), p_var[ip.dalpha].item()) - - self["alpha"] = self.dalpha * self.x - self["alpha_var"] = self["dalpha_var"] * self.x**2 - - self["talpha_fw_full"] = ( - ("x", "time"), - ip.get_taf_values( - pval=p_val, x=self.x.values, trans_att=self.trans_att.values, axis="" - ), - ) - self["talpha_fw_full_var"] = ( - ("x", "time"), - ip.get_taf_values( - pval=p_var, x=self.x.values, trans_att=self.trans_att.values, axis="" - ), - ) + coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} + params, param_covs = get_params_from_pval_single_ended(ip, coords, p_val=p_val, p_var=p_var, p_cov=p_cov, fix_alpha=fix_alpha) tmpf = self.gamma / ( - (np.log(self.st / self.ast) + (self.c + self["talpha_fw_full"])) - + self.alpha + (np.log(self.st / self.ast) + (params["c"] + params["talpha_fw_full"])) + + params["alpha"] ) self["tmpf"] = tmpf - 273.15 @@ -1729,58 +1616,36 @@ def calibration_single_ended( T_ta_fw=-(tmpf**2) / self.gamma, ) deriv_ds = xr.Dataset(deriv_dict) - # self["deriv"] = deriv_ds.to_array(dim="com2") - - sigma2_gamma_c = p_cov[np.ix_(ip.gamma, ip.c)] - sigma2_tafw_gamma = ip.get_taf_values( - pval=p_cov[ip.gamma], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ) - sigma2_tafw_c = ip.get_taf_values( - pval=p_cov[ip.c], - x=self.x.values, - trans_att=self.trans_att.values, - axis="time", - ) + var_fw_dict = dict( dT_dst=deriv_ds.T_st_fw**2 * parse_st_var(self, st_var, st_label="st"), dT_dast=deriv_ds.T_ast_fw**2 * parse_st_var(self, ast_var, st_label="ast"), - dT_gamma=deriv_ds.T_gamma_fw**2 * self["gamma_var"], - dT_dc=deriv_ds.T_c_fw**2 * self["c_var"], + dT_gamma=deriv_ds.T_gamma_fw**2 * param_covs["gamma"], + dT_dc=deriv_ds.T_c_fw**2 * param_covs["c"], dT_ddalpha=deriv_ds.T_alpha_fw**2 - * self["alpha_var"], # same as dT_dalpha - dT_dta=deriv_ds.T_ta_fw**2 * self["talpha_fw_full_var"], - dgamma_dc=(2 * deriv_ds.T_gamma_fw * deriv_ds.T_c_fw * sigma2_gamma_c), - dta_dgamma=(2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * sigma2_tafw_gamma), - dta_dc=(2 * deriv_ds.T_ta_fw * deriv_ds.T_c_fw * sigma2_tafw_c), + * param_covs["alpha"], # same as dT_dalpha + dT_dta=deriv_ds.T_ta_fw**2 * param_covs["talpha_fw_full"], + dgamma_dc=(2 * deriv_ds.T_gamma_fw * deriv_ds.T_c_fw * param_covs["gamma_c"]), + dta_dgamma=(2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * param_covs["tafw_gamma"]), + dta_dc=(2 * deriv_ds.T_ta_fw * deriv_ds.T_c_fw * param_covs["tafw_c"]), ) if not fix_alpha: # These correlations don't exist in case of fix_alpha. Including them reduces tmpf_var. - sigma2_gamma_dalpha = p_cov[np.ix_(ip.dalpha, ip.gamma)] - sigma2_dalpha_c = p_cov[np.ix_(ip.dalpha, ip.c)] - sigma2_tafw_dalpha = ip.get_taf_values( - pval=p_cov[ip.dalpha], - x=self.x.values, - trans_att=self.trans_att.values, - axis="", - ) var_fw_dict.update( dict( dgamma_ddalpha=( 2 * deriv_ds.T_gamma_fw * deriv_ds.T_dalpha_fw - * sigma2_gamma_dalpha + * param_covs["gamma_dalpha"] ), ddalpha_dc=( - 2 * deriv_ds.T_dalpha_fw * deriv_ds.T_c_fw * sigma2_dalpha_c + 2 * deriv_ds.T_dalpha_fw * deriv_ds.T_c_fw * param_covs["dalpha_c"] ), dta_ddalpha=( - 2 * deriv_ds.T_ta_fw * deriv_ds.T_dalpha_fw * sigma2_tafw_dalpha + 2 * deriv_ds.T_ta_fw * deriv_ds.T_dalpha_fw * param_covs["tafw_dalpha"] ), ) ) @@ -2121,8 +1986,8 @@ def calibration_double_ended( assert p_cov.shape == (ip.npar, ip.npar) coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} - params = get_params_from_pval(ip, coords, p_val=p_val) - param_covs = get_params_from_pval(ip, coords, p_val=p_var, p_cov=p_cov) + params = get_params_from_pval_double_ended(ip, coords, p_val=p_val) + param_covs = get_params_from_pval_double_ended(ip, coords, p_val=p_var, p_cov=p_cov) out = xr.Dataset( { diff --git a/src/dtscalibration/datastore_utils.py b/src/dtscalibration/datastore_utils.py index 339a45a7..c36c1c9b 100644 --- a/src/dtscalibration/datastore_utils.py +++ b/src/dtscalibration/datastore_utils.py @@ -527,7 +527,7 @@ def get_netcdf_encoding( return encoding -def get_params_from_pval(ip, coords, p_val=None, p_cov=None): +def get_params_from_pval_double_ended(ip, coords, p_val=None, p_cov=None): if p_val is not None: assert len(p_val) == ip.npar, "Length of p_val is incorrect" @@ -693,6 +693,80 @@ def get_params_from_pval(ip, coords, p_val=None, p_cov=None): ) # sigma2_tafw_tabw return params + + +def get_params_from_pval_single_ended(ip, coords, p_val=None, p_var=None, p_cov=None, fix_alpha=None): + assert len(p_val) == ip.npar, "Length of p_val is incorrect" + + params = xr.Dataset(coords=coords) + param_covs = xr.Dataset(coords=coords) + + params["gamma"] = (tuple(), p_val[ip.gamma].item()) + param_covs["gamma"] = (tuple(), p_var[ip.gamma].item()) + + if ip.nta > 0: + params["talpha_fw"] = ( + ("trans_att", "time"), + ip.get_taf_pars(p_val), + ) + param_covs["talpha_fw"] = ( + ("trans_att", "time"), + ip.get_taf_pars(p_var), + ) + else: + params["talpha_fw"] = (("trans_att", "time"), np.zeros((0, ip.nt))) + param_covs["talpha_fw"] = (("trans_att", "time"), np.zeros((0, ip.nt))) + + params["c"] = (("time",), p_val[ip.c]) + param_covs["c"] = (("time",), p_var[ip.c]) + + if fix_alpha is not None: + params["alpha"] = (("x",), fix_alpha[0]) + param_covs["alpha"] = (("x",), fix_alpha[1]) + + else: + params["dalpha"] = (tuple(), p_val[ip.dalpha].item()) + param_covs["dalpha"] = (tuple(), p_var[ip.dalpha].item()) + + params["alpha"] = params["dalpha"] * params["x"] + param_covs["alpha"] = param_covs["dalpha"] * params["x"]**2 + + param_covs["gamma_dalpha"] = p_cov[np.ix_(ip.dalpha, ip.gamma)] + param_covs["dalpha_c"] = p_cov[np.ix_(ip.dalpha, ip.c)] + param_covs["tafw_dalpha"] = ip.get_taf_values( + pval=p_cov[ip.dalpha], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ) + + params["talpha_fw_full"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_val, x=params["x"].values, trans_att=params["trans_att"].values, axis="" + ), + ) + param_covs["talpha_fw_full"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_var, x=params["x"].values, trans_att=params["trans_att"].values, axis="" + ), + ) + + param_covs["gamma_c"] = p_cov[np.ix_(ip.gamma, ip.c)] + param_covs["tafw_gamma"] = ip.get_taf_values( + pval=p_cov[ip.gamma], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ) + param_covs["tafw_c"] = ip.get_taf_values( + pval=p_cov[ip.c], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ) + return params, param_covs def merge_double_ended( From 0c427d805eeb2fc79c94a0a7bf5c609a4fd8d420 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 21 Aug 2023 09:47:03 +0200 Subject: [PATCH 14/22] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index b8f681da..f06bd184 100644 --- a/.gitignore +++ b/.gitignore @@ -72,3 +72,4 @@ docs/_build .appveyor.token *.bak .vscode/settings.json +*.code-workspace From 43f42c8849807e77a1c68fb081b6cfe406ab3af7 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Tue, 22 Aug 2023 21:18:41 +0200 Subject: [PATCH 15/22] Refactored calibration routines In preparation of the accessor --- pyproject.toml | 2 +- src/dtscalibration/averaging_utils.py | 5 +- .../calibration/section_utils.py | 1 - src/dtscalibration/datastore.py | 67 ++++++++++++------- src/dtscalibration/datastore_utils.py | 64 +++++++++++------- src/dtscalibration/io.py | 14 ++-- src/dtscalibration/io_utils.py | 4 +- 7 files changed, 92 insertions(+), 65 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 3597d280..af0bd1fd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["hatchling>=1.8.0", "hatch-vcs", "hatch-fancy-pypi-readme"] +requires = ["hatchling>=1.8.0", "hatch-vcs"] build-backend = "hatchling.build" [tool.hatch.version] diff --git a/src/dtscalibration/averaging_utils.py b/src/dtscalibration/averaging_utils.py index 9169c36e..c94194ba 100644 --- a/src/dtscalibration/averaging_utils.py +++ b/src/dtscalibration/averaging_utils.py @@ -1,5 +1,3 @@ - - def inverse_variance_weighted_mean( self, tmp1="tmpf", @@ -44,6 +42,7 @@ def inverse_variance_weighted_mean( pass + def inverse_variance_weighted_mean_array( self, tmp_label="tmpf", @@ -72,4 +71,4 @@ def inverse_variance_weighted_mean_array( 1 / self[tmp_var_label] ).sum(dim=dim) - pass \ No newline at end of file + pass diff --git a/src/dtscalibration/calibration/section_utils.py b/src/dtscalibration/calibration/section_utils.py index 648642be..2d1ef100 100644 --- a/src/dtscalibration/calibration/section_utils.py +++ b/src/dtscalibration/calibration/section_utils.py @@ -1,4 +1,3 @@ - import numpy as np import xarray as xr import yaml diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 3b6e9b49..9e1449c2 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -1597,23 +1597,26 @@ def calibration_single_ended( # store calibration parameters in DataStore coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} - params, param_covs = get_params_from_pval_single_ended(ip, coords, p_val=p_val, p_var=p_var, p_cov=p_cov, fix_alpha=fix_alpha) + params, param_covs = get_params_from_pval_single_ended( + ip, coords, p_val=p_val, p_var=p_var, p_cov=p_cov, fix_alpha=fix_alpha + ) - tmpf = self.gamma / ( + tmpf = params["gamma"] / ( (np.log(self.st / self.ast) + (params["c"] + params["talpha_fw_full"])) + params["alpha"] ) - self["tmpf"] = tmpf - 273.15 + out = xr.Dataset({"tmpf": tmpf - 273.15}) + out["tmpf"].attrs.update(_dim_attrs[("tmpf",)]) # tmpf_var deriv_dict = dict( - T_gamma_fw=tmpf / self.gamma, - T_st_fw=-(tmpf**2) / (self.gamma * self.st), - T_ast_fw=tmpf**2 / (self.gamma * self.ast), - T_c_fw=-(tmpf**2) / self.gamma, - T_dalpha_fw=-self.x * (tmpf**2) / self.gamma, - T_alpha_fw=-(tmpf**2) / self.gamma, - T_ta_fw=-(tmpf**2) / self.gamma, + T_gamma_fw=tmpf / params["gamma"], + T_st_fw=-(tmpf**2) / (params["gamma"] * self.st), + T_ast_fw=tmpf**2 / (params["gamma"] * self.ast), + T_c_fw=-(tmpf**2) / params["gamma"], + T_dalpha_fw=-self.x * (tmpf**2) / params["gamma"], + T_alpha_fw=-(tmpf**2) / params["gamma"], + T_ta_fw=-(tmpf**2) / params["gamma"], ) deriv_ds = xr.Dataset(deriv_dict) @@ -1626,8 +1629,12 @@ def calibration_single_ended( dT_ddalpha=deriv_ds.T_alpha_fw**2 * param_covs["alpha"], # same as dT_dalpha dT_dta=deriv_ds.T_ta_fw**2 * param_covs["talpha_fw_full"], - dgamma_dc=(2 * deriv_ds.T_gamma_fw * deriv_ds.T_c_fw * param_covs["gamma_c"]), - dta_dgamma=(2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * param_covs["tafw_gamma"]), + dgamma_dc=( + 2 * deriv_ds.T_gamma_fw * deriv_ds.T_c_fw * param_covs["gamma_c"] + ), + dta_dgamma=( + 2 * deriv_ds.T_ta_fw * deriv_ds.T_gamma_fw * param_covs["tafw_gamma"] + ), dta_dc=(2 * deriv_ds.T_ta_fw * deriv_ds.T_c_fw * param_covs["tafw_c"]), ) @@ -1642,19 +1649,23 @@ def calibration_single_ended( * param_covs["gamma_dalpha"] ), ddalpha_dc=( - 2 * deriv_ds.T_dalpha_fw * deriv_ds.T_c_fw * param_covs["dalpha_c"] + 2 + * deriv_ds.T_dalpha_fw + * deriv_ds.T_c_fw + * param_covs["dalpha_c"] ), dta_ddalpha=( - 2 * deriv_ds.T_ta_fw * deriv_ds.T_dalpha_fw * param_covs["tafw_dalpha"] + 2 + * deriv_ds.T_ta_fw + * deriv_ds.T_dalpha_fw + * param_covs["tafw_dalpha"] ), ) ) - self["var_fw_da"] = xr.Dataset(var_fw_dict).to_array(dim="comp_fw") - self["tmpf_var"] = self["var_fw_da"].sum(dim="comp_fw") - - self["tmpf"].attrs.update(_dim_attrs[("tmpf",)]) - self["tmpf_var"].attrs.update(_dim_attrs[("tmpf_var",)]) + out["var_fw_da"] = xr.Dataset(var_fw_dict).to_array(dim="comp_fw") + out["tmpf_var"] = out["var_fw_da"].sum(dim="comp_fw") + out["tmpf_var"].attrs.update(_dim_attrs[("tmpf_var",)]) drop_vars = [ k for k, v in self.items() if {"params1", "params2"}.intersection(v.dims) @@ -1663,9 +1674,14 @@ def calibration_single_ended( for k in drop_vars: del self[k] - self["p_val"] = (("params1",), p_val) - self["p_cov"] = (("params1", "params2"), p_cov) + out["p_val"] = (("params1",), p_val) + out["p_cov"] = (("params1", "params2"), p_cov) + + out.update(params) + for key, dataarray in param_covs.data_vars.items(): + out[key + "_var"] = dataarray + self.update(out) pass def calibration_double_ended( @@ -1987,7 +2003,9 @@ def calibration_double_ended( coords = {"x": self["x"], "time": self["time"], "trans_att": self["trans_att"]} params = get_params_from_pval_double_ended(ip, coords, p_val=p_val) - param_covs = get_params_from_pval_double_ended(ip, coords, p_val=p_var, p_cov=p_cov) + param_covs = get_params_from_pval_double_ended( + ip, coords, p_val=p_var, p_cov=p_cov + ) out = xr.Dataset( { @@ -2203,9 +2221,8 @@ def calibration_double_ended( out["p_cov"] = (("params1", "params2"), p_cov) out.update(params) - - for key, da in param_covs.data_vars.items(): - out[key + "_var"] = da + for key, dataarray in param_covs.data_vars.items(): + out[key + "_var"] = dataarray self.update(out) pass diff --git a/src/dtscalibration/datastore_utils.py b/src/dtscalibration/datastore_utils.py index c36c1c9b..0d54dd31 100644 --- a/src/dtscalibration/datastore_utils.py +++ b/src/dtscalibration/datastore_utils.py @@ -693,9 +693,11 @@ def get_params_from_pval_double_ended(ip, coords, p_val=None, p_cov=None): ) # sigma2_tafw_tabw return params - -def get_params_from_pval_single_ended(ip, coords, p_val=None, p_var=None, p_cov=None, fix_alpha=None): + +def get_params_from_pval_single_ended( + ip, coords, p_val=None, p_var=None, p_cov=None, fix_alpha=None +): assert len(p_val) == ip.npar, "Length of p_val is incorrect" params = xr.Dataset(coords=coords) @@ -729,42 +731,56 @@ def get_params_from_pval_single_ended(ip, coords, p_val=None, p_var=None, p_cov= param_covs["dalpha"] = (tuple(), p_var[ip.dalpha].item()) params["alpha"] = params["dalpha"] * params["x"] - param_covs["alpha"] = param_covs["dalpha"] * params["x"]**2 + param_covs["alpha"] = param_covs["dalpha"] * params["x"] ** 2 - param_covs["gamma_dalpha"] = p_cov[np.ix_(ip.dalpha, ip.gamma)] - param_covs["dalpha_c"] = p_cov[np.ix_(ip.dalpha, ip.c)] - param_covs["tafw_dalpha"] = ip.get_taf_values( - pval=p_cov[ip.dalpha], - x=params["x"].values, - trans_att=params["trans_att"].values, - axis="", + param_covs["gamma_dalpha"] = (tuple(), p_cov[np.ix_(ip.dalpha, ip.gamma)][0, 0]) + param_covs["dalpha_c"] = (("time",), p_cov[np.ix_(ip.dalpha, ip.c)][0, :]) + param_covs["tafw_dalpha"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_cov[ip.dalpha], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ), ) params["talpha_fw_full"] = ( ("x", "time"), ip.get_taf_values( - pval=p_val, x=params["x"].values, trans_att=params["trans_att"].values, axis="" + pval=p_val, + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", ), ) param_covs["talpha_fw_full"] = ( ("x", "time"), ip.get_taf_values( - pval=p_var, x=params["x"].values, trans_att=params["trans_att"].values, axis="" + pval=p_var, + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", ), ) - - param_covs["gamma_c"] = p_cov[np.ix_(ip.gamma, ip.c)] - param_covs["tafw_gamma"] = ip.get_taf_values( - pval=p_cov[ip.gamma], - x=params["x"].values, - trans_att=params["trans_att"].values, - axis="", + param_covs["gamma_c"] = (("time",), p_cov[np.ix_(ip.gamma, ip.c)][0, :]) + param_covs["tafw_gamma"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_cov[ip.gamma], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="", + ), ) - param_covs["tafw_c"] = ip.get_taf_values( - pval=p_cov[ip.c], - x=params["x"].values, - trans_att=params["trans_att"].values, - axis="time", + param_covs["tafw_c"] = ( + ("x", "time"), + ip.get_taf_values( + pval=p_cov[ip.c], + x=params["x"].values, + trans_att=params["trans_att"].values, + axis="time", + ), ) return params, param_covs diff --git a/src/dtscalibration/io.py b/src/dtscalibration/io.py index edbeac63..413851f5 100644 --- a/src/dtscalibration/io.py +++ b/src/dtscalibration/io.py @@ -423,10 +423,9 @@ def read_apsensing_files( else: warnings.warn( - "AP sensing device " - '"{0}"'.format(device) - + " has not been tested.\nPlease open an issue on github" - + " and provide an example file" + "AP sensing device {device}" + " has not been tested.\nPlease open an issue on github" + " and provide an example file" ) data_vars, coords, attrs = read_apsensing_files_routine( @@ -544,10 +543,9 @@ def read_sensornet_files( else: flip_reverse_measurements = False warnings.warn( - "Sensornet .dff version " - '"{0}"'.format(ddf_version) - + " has not been tested.\nPlease open an issue on github" - + " and provide an example file" + f"Sensornet .dff version {ddf_version}" + " has not been tested.\nPlease open an issue on github" + " and provide an example file" ) data_vars, coords, attrs = read_sensornet_files_routine_v3( diff --git a/src/dtscalibration/io_utils.py b/src/dtscalibration/io_utils.py index 4b480742..ca5dbb3a 100644 --- a/src/dtscalibration/io_utils.py +++ b/src/dtscalibration/io_utils.py @@ -1589,9 +1589,7 @@ def grab_data_per_file(file_handle): elif name in dim_attrs_apsensing: data_vars[tld[name]] = (["x", "time"], data_arri, dim_attrs_apsensing[name]) else: - raise ValueError( - "Dont know what to do with the" + f" {name} data column" - ) + raise ValueError("Dont know what to do with the" + f" {name} data column") _time_dtype = [("filename_tstamp", np.int64), ("acquisitionTime", " Date: Tue, 22 Aug 2023 21:34:43 +0200 Subject: [PATCH 16/22] Please ruff --- src/dtscalibration/datastore.py | 2 +- src/dtscalibration/plot.py | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 9e1449c2..1882aa64 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -155,7 +155,7 @@ def __repr__(self): # print sections vl = [ - "{0:.2f}{2} - {1:.2f}{2}".format(vi.start, vi.stop, unit) + f"{vi.start:.2f}{unit} - {vi.stop:.2f}{unit}" for vi in v ] preamble_new += " and ".join(vl) + "\n" diff --git a/src/dtscalibration/plot.py b/src/dtscalibration/plot.py index 7e487879..13c0c111 100755 --- a/src/dtscalibration/plot.py +++ b/src/dtscalibration/plot.py @@ -630,8 +630,7 @@ def plot_sigma_report( ) else: ax1.set_title( - "Projected uncertainty at t={} compared to standard error " - "in baths".format(itimes) + f"Projected uncertainty at t={itimes} compared to standard error in baths" ) ax1.legend() ax1.set_ylabel(r"Temperature [$^\circ$C]") From 49bbae3f820a567e546a24da2dbd92188d17eb75 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 3 Sep 2023 15:49:12 +0200 Subject: [PATCH 17/22] Refactored Datastore class --- src/dtscalibration/datastore.py | 580 +++++++++++++++++--------------- tests/test_dtscalibration.py | 6 +- 2 files changed, 306 insertions(+), 280 deletions(-) diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 1882aa64..db555798 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -154,10 +154,7 @@ def __repr__(self): preamble_new += sec_stat # print sections - vl = [ - f"{vi.start:.2f}{unit} - {vi.stop:.2f}{unit}" - for vi in v - ] + vl = [f"{vi.start:.2f}{unit} - {vi.stop:.2f}{unit}" for vi in v] preamble_new += " and ".join(vl) + "\n" else: @@ -2376,7 +2373,9 @@ def average_single_ended( if var_only_sections is not None: raise NotImplementedError() - self.conf_int_single_ended( + out = xr.Dataset() + + mcparams = self.conf_int_single_ended( p_val=p_val, p_cov=p_cov, st_var=st_var, @@ -2388,85 +2387,86 @@ def average_single_ended( reduce_memory_usage=reduce_memory_usage, **kwargs, ) + mcparams["tmpf"] = self["tmpf"] if ci_avg_time_sel is not None: time_dim2 = "time" + "_avg" x_dim2 = "x" - self.coords[time_dim2] = ( + mcparams.coords[time_dim2] = ( (time_dim2,), - self["time"].sel(**{"time": ci_avg_time_sel}).data, + mcparams["time"].sel(**{"time": ci_avg_time_sel}).data, ) - self["tmpf_avgsec"] = ( + mcparams["tmpf_avgsec"] = ( ("x", time_dim2), - self["tmpf"].sel(**{"time": ci_avg_time_sel}).data, + mcparams["tmpf"].sel(**{"time": ci_avg_time_sel}).data, ) - self["tmpf_mc_set"] = ( + mcparams["tmpf_mc_set"] = ( ("mc", "x", time_dim2), - self["tmpf" + "_mc_set"].sel(**{"time": ci_avg_time_sel}).data, + mcparams["tmpf" + "_mc_set"].sel(**{"time": ci_avg_time_sel}).data, ) elif ci_avg_time_isel is not None: time_dim2 = "time" + "_avg" x_dim2 = "x" - self.coords[time_dim2] = ( + mcparams.coords[time_dim2] = ( (time_dim2,), - self["time"].isel(**{"time": ci_avg_time_isel}).data, + mcparams["time"].isel(**{"time": ci_avg_time_isel}).data, ) - self["tmpf_avgsec"] = ( + mcparams["tmpf_avgsec"] = ( ("x", time_dim2), - self["tmpf"].isel(**{"time": ci_avg_time_isel}).data, + mcparams["tmpf"].isel(**{"time": ci_avg_time_isel}).data, ) - self["tmpf_mc_set"] = ( + mcparams["tmpf_mc_set"] = ( ("mc", "x", time_dim2), - self["tmpf" + "_mc_set"].isel(**{"time": ci_avg_time_isel}).data, + mcparams["tmpf" + "_mc_set"].isel(**{"time": ci_avg_time_isel}).data, ) elif ci_avg_x_sel is not None: time_dim2 = "time" x_dim2 = "x_avg" - self.coords[x_dim2] = ((x_dim2,), self.x.sel(x=ci_avg_x_sel).data) - self["tmpf_avgsec"] = ( + mcparams.coords[x_dim2] = ((x_dim2,), mcparams.x.sel(x=ci_avg_x_sel).data) + mcparams["tmpf_avgsec"] = ( (x_dim2, "time"), - self["tmpf"].sel(x=ci_avg_x_sel).data, + mcparams["tmpf"].sel(x=ci_avg_x_sel).data, ) - self["tmpf_mc_set"] = ( + mcparams["tmpf_mc_set"] = ( ("mc", x_dim2, "time"), - self["tmpf_mc_set"].sel(x=ci_avg_x_sel).data, + mcparams["tmpf_mc_set"].sel(x=ci_avg_x_sel).data, ) elif ci_avg_x_isel is not None: time_dim2 = "time" x_dim2 = "x_avg" - self.coords[x_dim2] = ((x_dim2,), self.x.isel(x=ci_avg_x_isel).data) - self["tmpf_avgsec"] = ( + mcparams.coords[x_dim2] = ((x_dim2,), mcparams.x.isel(x=ci_avg_x_isel).data) + mcparams["tmpf_avgsec"] = ( (x_dim2, time_dim2), - self["tmpf"].isel(x=ci_avg_x_isel).data, + mcparams["tmpf"].isel(x=ci_avg_x_isel).data, ) - self["tmpf_mc_set"] = ( + mcparams["tmpf_mc_set"] = ( ("mc", x_dim2, time_dim2), - self["tmpf_mc_set"].isel(x=ci_avg_x_isel).data, + mcparams["tmpf_mc_set"].isel(x=ci_avg_x_isel).data, ) else: - self["tmpf_avgsec"] = self["tmpf"] + mcparams["tmpf_avgsec"] = mcparams["tmpf"] x_dim2 = "x" time_dim2 = "time" # subtract the mean temperature - q = self["tmpf_mc_set"] - self["tmpf_avgsec"] - self["tmpf_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) + q = mcparams["tmpf_mc_set"] - mcparams["tmpf_avgsec"] + out["tmpf_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) if ci_avg_x_flag1: # unweighted mean - self["tmpf_avgx1"] = self["tmpf" + "_avgsec"].mean(dim=x_dim2) + out["tmpf_avgx1"] = mcparams["tmpf" + "_avgsec"].mean(dim=x_dim2) - q = self["tmpf_mc_set"] - self["tmpf_avgsec"] + q = mcparams["tmpf_mc_set"] - mcparams["tmpf_avgsec"] qvar = q.var(dim=["mc", x_dim2], ddof=1) - self["tmpf_mc_avgx1_var"] = qvar + out["tmpf_mc_avgx1_var"] = qvar if conf_ints: - new_chunks = (len(conf_ints), self["tmpf_mc_set"].chunks[2]) - avg_axis = self["tmpf_mc_set"].get_axis_num(["mc", x_dim2]) - q = self["tmpf_mc_set"].data.map_blocks( + new_chunks = (len(conf_ints), mcparams["tmpf_mc_set"].chunks[2]) + avg_axis = mcparams["tmpf_mc_set"].get_axis_num(["mc", x_dim2]) + q = mcparams["tmpf_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks, # drop_axis=avg_axis, @@ -2474,28 +2474,28 @@ def average_single_ended( new_axis=0, ) # The new CI dim is added as firsaxis - self["tmpf_mc_avgx1"] = (("CI", time_dim2), q) + out["tmpf_mc_avgx1"] = (("CI", time_dim2), q) if ci_avg_x_flag2: - q = self["tmpf_mc_set"] - self["tmpf_avgsec"] + q = mcparams["tmpf_mc_set"] - mcparams["tmpf_avgsec"] qvar = q.var(dim=["mc"], ddof=1) # Inverse-variance weighting avg_x_var = 1 / (1 / qvar).sum(dim=x_dim2) - self["tmpf_mc_avgx2_var"] = avg_x_var + out["tmpf_mc_avgx2_var"] = avg_x_var - self["tmpf" + "_mc_avgx2_set"] = (self["tmpf_mc_set"] / qvar).sum( + mcparams["tmpf" + "_mc_avgx2_set"] = (mcparams["tmpf_mc_set"] / qvar).sum( dim=x_dim2 ) * avg_x_var - self["tmpf" + "_avgx2"] = self["tmpf" + "_mc_avgx2_set"].mean(dim="mc") + out["tmpf" + "_avgx2"] = mcparams["tmpf" + "_mc_avgx2_set"].mean(dim="mc") if conf_ints: - new_chunks = (len(conf_ints), self["tmpf_mc_set"].chunks[2]) - avg_axis_avgx = self["tmpf_mc_set"].get_axis_num("mc") + new_chunks = (len(conf_ints), mcparams["tmpf_mc_set"].chunks[2]) + avg_axis_avgx = mcparams["tmpf_mc_set"].get_axis_num("mc") - qq = self["tmpf_mc_avgx2_set"].data.map_blocks( + qq = mcparams["tmpf_mc_avgx2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avgx), chunks=new_chunks, # drop_axis=avg_axis_avgx, @@ -2504,20 +2504,20 @@ def average_single_ended( dtype=float, ) # The new CI dimension is added as # firsaxis - self["tmpf_mc_avgx2"] = (("CI", time_dim2), qq) + out["tmpf_mc_avgx2"] = (("CI", time_dim2), qq) if ci_avg_time_flag1 is not None: # unweighted mean - self["tmpf_avg1"] = self["tmpf" + "_avgsec"].mean(dim=time_dim2) + out["tmpf_avg1"] = mcparams["tmpf_avgsec"].mean(dim=time_dim2) - q = self["tmpf_mc_set"] - self["tmpf_avgsec"] + q = mcparams["tmpf_mc_set"] - mcparams["tmpf_avgsec"] qvar = q.var(dim=["mc", time_dim2], ddof=1) - self["tmpf_mc_avg1_var"] = qvar + out["tmpf_mc_avg1_var"] = qvar if conf_ints: - new_chunks = (len(conf_ints), self["tmpf_mc_set"].chunks[1]) - avg_axis = self["tmpf_mc_set"].get_axis_num(["mc", time_dim2]) - q = self["tmpf_mc_set"].data.map_blocks( + new_chunks = (len(conf_ints), mcparams["tmpf_mc_set"].chunks[1]) + avg_axis = mcparams["tmpf_mc_set"].get_axis_num(["mc", time_dim2]) + q = mcparams["tmpf_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks, # drop_axis=avg_axis, @@ -2525,28 +2525,28 @@ def average_single_ended( new_axis=0, ) # The new CI dim is added as firsaxis - self["tmpf_mc_avg1"] = (("CI", x_dim2), q) + out["tmpf_mc_avg1"] = (("CI", x_dim2), q) if ci_avg_time_flag2: - q = self["tmpf_mc_set"] - self["tmpf_avgsec"] + q = mcparams["tmpf_mc_set"] - mcparams["tmpf_avgsec"] qvar = q.var(dim=["mc"], ddof=1) # Inverse-variance weighting avg_time_var = 1 / (1 / qvar).sum(dim=time_dim2) - self["tmpf_mc_avg2_var"] = avg_time_var + out["tmpf_mc_avg2_var"] = avg_time_var - self["tmpf" + "_mc_avg2_set"] = (self["tmpf_mc_set"] / qvar).sum( + mcparams["tmpf" + "_mc_avg2_set"] = (mcparams["tmpf_mc_set"] / qvar).sum( dim=time_dim2 ) * avg_time_var - self["tmpf_avg2"] = self["tmpf" + "_mc_avg2_set"].mean(dim="mc") + out["tmpf_avg2"] = mcparams["tmpf" + "_mc_avg2_set"].mean(dim="mc") if conf_ints: - new_chunks = (len(conf_ints), self["tmpf_mc_set"].chunks[1]) - avg_axis_avg2 = self["tmpf_mc_set"].get_axis_num("mc") + new_chunks = (len(conf_ints), mcparams["tmpf_mc_set"].chunks[1]) + avg_axis_avg2 = mcparams["tmpf_mc_set"].get_axis_num("mc") - qq = self["tmpf_mc_avg2_set"].data.map_blocks( + qq = mcparams["tmpf_mc_avg2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avg2), chunks=new_chunks, # drop_axis=avg_axis_avg2, @@ -2555,7 +2555,7 @@ def average_single_ended( dtype=float, ) # The new CI dimension is added as # firsaxis - self["tmpf_mc_avg2"] = (("CI", x_dim2), qq) + out["tmpf_mc_avg2"] = (("CI", x_dim2), qq) # Clean up the garbage. All arrays with a Monte Carlo dimension. if mc_remove_set_flag: @@ -2577,8 +2577,10 @@ def average_single_ended( remove_mc_set.append("tmpf_mc_avgsec_var") for k in remove_mc_set: - if k in self: - del self[k] + if k in out: + del out[k] + + self.update(out) pass def average_double_ended( @@ -2767,7 +2769,9 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): else: pass - self.conf_int_double_ended( + out = xr.Dataset() + + mcparams = self.conf_int_double_ended( sections=sections, p_val=p_val, p_cov=p_cov, @@ -2787,83 +2791,89 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): if ci_avg_time_sel is not None: time_dim2 = "time" + "_avg" x_dim2 = "x" - self.coords[time_dim2] = ( + mcparams.coords[time_dim2] = ( (time_dim2,), - self["time"].sel(**{"time": ci_avg_time_sel}).data, + mcparams["time"].sel(**{"time": ci_avg_time_sel}).data, ) - self[label + "_avgsec"] = ( + mcparams[label + "_avgsec"] = ( ("x", time_dim2), self[label].sel(**{"time": ci_avg_time_sel}).data, ) - self[label + "_mc_set"] = ( + mcparams[label + "_mc_set"] = ( ("mc", "x", time_dim2), - self[label + "_mc_set"].sel(**{"time": ci_avg_time_sel}).data, + mcparams[label + "_mc_set"].sel(**{"time": ci_avg_time_sel}).data, ) elif ci_avg_time_isel is not None: time_dim2 = "time" + "_avg" x_dim2 = "x" - self.coords[time_dim2] = ( + mcparams.coords[time_dim2] = ( (time_dim2,), - self["time"].isel(**{"time": ci_avg_time_isel}).data, + mcparams["time"].isel(**{"time": ci_avg_time_isel}).data, ) - self[label + "_avgsec"] = ( + mcparams[label + "_avgsec"] = ( ("x", time_dim2), self[label].isel(**{"time": ci_avg_time_isel}).data, ) - self[label + "_mc_set"] = ( + mcparams[label + "_mc_set"] = ( ("mc", "x", time_dim2), - self[label + "_mc_set"].isel(**{"time": ci_avg_time_isel}).data, + mcparams[label + "_mc_set"].isel(**{"time": ci_avg_time_isel}).data, ) elif ci_avg_x_sel is not None: time_dim2 = "time" x_dim2 = "x_avg" - self.coords[x_dim2] = ((x_dim2,), self.x.sel(x=ci_avg_x_sel).data) - self[label + "_avgsec"] = ( + mcparams.coords[x_dim2] = ( + (x_dim2,), + mcparams.x.sel(x=ci_avg_x_sel).data, + ) + mcparams[label + "_avgsec"] = ( (x_dim2, "time"), self[label].sel(x=ci_avg_x_sel).data, ) - self[label + "_mc_set"] = ( + mcparams[label + "_mc_set"] = ( ("mc", x_dim2, "time"), - self[label + "_mc_set"].sel(x=ci_avg_x_sel).data, + mcparams[label + "_mc_set"].sel(x=ci_avg_x_sel).data, ) elif ci_avg_x_isel is not None: time_dim2 = "time" x_dim2 = "x_avg" - self.coords[x_dim2] = ((x_dim2,), self.x.isel(x=ci_avg_x_isel).data) - self[label + "_avgsec"] = ( + mcparams.coords[x_dim2] = ( + (x_dim2,), + mcparams.x.isel(x=ci_avg_x_isel).data, + ) + mcparams[label + "_avgsec"] = ( (x_dim2, time_dim2), self[label].isel(x=ci_avg_x_isel).data, ) - self[label + "_mc_set"] = ( + mcparams[label + "_mc_set"] = ( ("mc", x_dim2, time_dim2), - self[label + "_mc_set"].isel(x=ci_avg_x_isel).data, + mcparams[label + "_mc_set"].isel(x=ci_avg_x_isel).data, ) else: - self[label + "_avgsec"] = self[label] + mcparams[label + "_avgsec"] = self[label] x_dim2 = "x" time_dim2 = "time" - memchunk = self[label + "_mc_set"].chunks + memchunk = mcparams[label + "_mc_set"].chunks # subtract the mean temperature - q = self[label + "_mc_set"] - self[label + "_avgsec"] - self[label + "_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) + q = mcparams[label + "_mc_set"] - mcparams[label + "_avgsec"] + out[label + "_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) if ci_avg_x_flag1: # unweighted mean - self[label + "_avgx1"] = self[label + "_avgsec"].mean(dim=x_dim2) + out[label + "_avgx1"] = mcparams[label + "_avgsec"].mean(dim=x_dim2) - q = self[label + "_mc_set"] - self[label + "_avgsec"] + q = mcparams[label + "_mc_set"] - mcparams[label + "_avgsec"] qvar = q.var(dim=["mc", x_dim2], ddof=1) - self[label + "_mc_avgx1_var"] = qvar + out[label + "_mc_avgx1_var"] = qvar if conf_ints: - new_chunks = (len(conf_ints), self[label + "_mc_set"].chunks[2]) - avg_axis = self[label + "_mc_set"].get_axis_num(["mc", x_dim2]) - q = self[label + "_mc_set"].data.map_blocks( + new_chunks = (len(conf_ints), mcparams[label + "_mc_set"].chunks[2]) + avg_axis = mcparams[label + "_mc_set"].get_axis_num(["mc", x_dim2]) + q = mcparams[label + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks, # drop_axis=avg_axis, @@ -2871,28 +2881,28 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): new_axis=0, ) # The new CI dim is added as firsaxis - self[label + "_mc_avgx1"] = (("CI", time_dim2), q) + out[label + "_mc_avgx1"] = (("CI", time_dim2), q) if ci_avg_x_flag2: - q = self[label + "_mc_set"] - self[label + "_avgsec"] + q = mcparams[label + "_mc_set"] - mcparams[label + "_avgsec"] qvar = q.var(dim=["mc"], ddof=1) # Inverse-variance weighting avg_x_var = 1 / (1 / qvar).sum(dim=x_dim2) - self[label + "_mc_avgx2_var"] = avg_x_var + out[label + "_mc_avgx2_var"] = avg_x_var - self[label + "_mc_avgx2_set"] = (self[label + "_mc_set"] / qvar).sum( - dim=x_dim2 - ) * avg_x_var - self[label + "_avgx2"] = self[label + "_mc_avgx2_set"].mean(dim="mc") + mcparams[label + "_mc_avgx2_set"] = ( + mcparams[label + "_mc_set"] / qvar + ).sum(dim=x_dim2) * avg_x_var + out[label + "_avgx2"] = mcparams[label + "_mc_avgx2_set"].mean(dim="mc") if conf_ints: - new_chunks = (len(conf_ints), self[label + "_mc_set"].chunks[2]) - avg_axis_avgx = self[label + "_mc_set"].get_axis_num("mc") + new_chunks = (len(conf_ints), mcparams[label + "_mc_set"].chunks[2]) + avg_axis_avgx = mcparams[label + "_mc_set"].get_axis_num("mc") - qq = self[label + "_mc_avgx2_set"].data.map_blocks( + qq = mcparams[label + "_mc_avgx2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avgx), chunks=new_chunks, # drop_axis=avg_axis_avgx, @@ -2901,20 +2911,22 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): dtype=float, ) # The new CI dimension is added as # firsaxis - self[label + "_mc_avgx2"] = (("CI", time_dim2), qq) + out[label + "_mc_avgx2"] = (("CI", time_dim2), qq) if ci_avg_time_flag1 is not None: # unweighted mean - self[label + "_avg1"] = self[label + "_avgsec"].mean(dim=time_dim2) + out[label + "_avg1"] = mcparams[label + "_avgsec"].mean(dim=time_dim2) - q = self[label + "_mc_set"] - self[label + "_avgsec"] + q = mcparams[label + "_mc_set"] - mcparams[label + "_avgsec"] qvar = q.var(dim=["mc", time_dim2], ddof=1) - self[label + "_mc_avg1_var"] = qvar + out[label + "_mc_avg1_var"] = qvar if conf_ints: - new_chunks = (len(conf_ints), self[label + "_mc_set"].chunks[1]) - avg_axis = self[label + "_mc_set"].get_axis_num(["mc", time_dim2]) - q = self[label + "_mc_set"].data.map_blocks( + new_chunks = (len(conf_ints), mcparams[label + "_mc_set"].chunks[1]) + avg_axis = mcparams[label + "_mc_set"].get_axis_num( + ["mc", time_dim2] + ) + q = mcparams[label + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks, # drop_axis=avg_axis, @@ -2922,28 +2934,28 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): new_axis=0, ) # The new CI dim is added as firsaxis - self[label + "_mc_avg1"] = (("CI", x_dim2), q) + out[label + "_mc_avg1"] = (("CI", x_dim2), q) if ci_avg_time_flag2: - q = self[label + "_mc_set"] - self[label + "_avgsec"] + q = mcparams[label + "_mc_set"] - mcparams[label + "_avgsec"] qvar = q.var(dim=["mc"], ddof=1) # Inverse-variance weighting avg_time_var = 1 / (1 / qvar).sum(dim=time_dim2) - self[label + "_mc_avg2_var"] = avg_time_var + out[label + "_mc_avg2_var"] = avg_time_var - self[label + "_mc_avg2_set"] = (self[label + "_mc_set"] / qvar).sum( - dim=time_dim2 - ) * avg_time_var - self[label + "_avg2"] = self[label + "_mc_avg2_set"].mean(dim="mc") + mcparams[label + "_mc_avg2_set"] = ( + mcparams[label + "_mc_set"] / qvar + ).sum(dim=time_dim2) * avg_time_var + out[label + "_avg2"] = mcparams[label + "_mc_avg2_set"].mean(dim="mc") if conf_ints: - new_chunks = (len(conf_ints), self[label + "_mc_set"].chunks[1]) - avg_axis_avg2 = self[label + "_mc_set"].get_axis_num("mc") + new_chunks = (len(conf_ints), mcparams[label + "_mc_set"].chunks[1]) + avg_axis_avg2 = mcparams[label + "_mc_set"].get_axis_num("mc") - qq = self[label + "_mc_avg2_set"].data.map_blocks( + qq = mcparams[label + "_mc_avg2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avg2), chunks=new_chunks, # drop_axis=avg_axis_avg2, @@ -2952,40 +2964,40 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): dtype=float, ) # The new CI dimension is added as # firsaxis - self[label + "_mc_avg2"] = (("CI", x_dim2), qq) + out[label + "_mc_avg2"] = (("CI", x_dim2), qq) # Weighted mean of the forward and backward tmpw_var = 1 / ( - 1 / self["tmpf_mc" + "_avgsec_var"] + 1 / self["tmpb_mc" + "_avgsec_var"] + 1 / out["tmpf_mc" + "_avgsec_var"] + 1 / out["tmpb_mc" + "_avgsec_var"] ) q = ( - self["tmpf_mc_set"] / self["tmpf_mc" + "_avgsec_var"] - + self["tmpb_mc_set"] / self["tmpb_mc" + "_avgsec_var"] + mcparams["tmpf_mc_set"] / out["tmpf_mc" + "_avgsec_var"] + + mcparams["tmpb_mc_set"] / out["tmpb_mc" + "_avgsec_var"] ) * tmpw_var - self["tmpw" + "_mc_set"] = q # + mcparams["tmpw" + "_mc_set"] = q # - # self["tmpw"] = self["tmpw" + '_mc_set'].mean(dim='mc') - self["tmpw" + "_avgsec"] = ( - self["tmpf_avgsec"] / self["tmpf_mc" + "_avgsec_var"] - + self["tmpb_avgsec"] / self["tmpb_mc" + "_avgsec_var"] + # out["tmpw"] = out["tmpw" + '_mc_set'].mean(dim='mc') + out["tmpw" + "_avgsec"] = ( + mcparams["tmpf_avgsec"] / out["tmpf_mc" + "_avgsec_var"] + + mcparams["tmpb_avgsec"] / out["tmpb_mc" + "_avgsec_var"] ) * tmpw_var - q = self["tmpw" + "_mc_set"] - self["tmpw" + "_avgsec"] - self["tmpw" + "_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) + q = mcparams["tmpw" + "_mc_set"] - out["tmpw_avgsec"] + out["tmpw" + "_mc" + "_avgsec_var"] = q.var(dim="mc", ddof=1) if ci_avg_time_flag1: - self["tmpw" + "_avg1"] = self["tmpw" + "_avgsec"].mean(dim=time_dim2) + out["tmpw" + "_avg1"] = out["tmpw" + "_avgsec"].mean(dim=time_dim2) - self["tmpw" + "_mc_avg1_var"] = self["tmpw" + "_mc_set"].var( + out["tmpw" + "_mc_avg1_var"] = mcparams["tmpw" + "_mc_set"].var( dim=["mc", time_dim2] ) if conf_ints: new_chunks_weighted = ((len(conf_ints),),) + (memchunk[1],) - avg_axis = self["tmpw" + "_mc_set"].get_axis_num(["mc", time_dim2]) - q2 = self["tmpw" + "_mc_set"].data.map_blocks( + avg_axis = mcparams["tmpw" + "_mc_set"].get_axis_num(["mc", time_dim2]) + q2 = mcparams["tmpw" + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks_weighted, # Explicitly define output chunks @@ -2994,32 +3006,32 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): dtype=float, ) # The new CI dimension is added as # first axis - self["tmpw" + "_mc_avg1"] = (("CI", x_dim2), q2) + out["tmpw" + "_mc_avg1"] = (("CI", x_dim2), q2) if ci_avg_time_flag2: tmpw_var_avg2 = 1 / ( - 1 / self["tmpf_mc_avg2_var"] + 1 / self["tmpb_mc_avg2_var"] + 1 / out["tmpf_mc_avg2_var"] + 1 / out["tmpb_mc_avg2_var"] ) q = ( - self["tmpf_mc_avg2_set"] / self["tmpf_mc_avg2_var"] - + self["tmpb_mc_avg2_set"] / self["tmpb_mc_avg2_var"] + mcparams["tmpf_mc_avg2_set"] / out["tmpf_mc_avg2_var"] + + mcparams["tmpb_mc_avg2_set"] / out["tmpb_mc_avg2_var"] ) * tmpw_var_avg2 - self["tmpw" + "_mc_avg2_set"] = q # + mcparams["tmpw" + "_mc_avg2_set"] = q # - self["tmpw" + "_avg2"] = ( - self["tmpf_avg2"] / self["tmpf_mc_avg2_var"] - + self["tmpb_avg2"] / self["tmpb_mc_avg2_var"] + out["tmpw" + "_avg2"] = ( + out["tmpf_avg2"] / out["tmpf_mc_avg2_var"] + + out["tmpb_avg2"] / out["tmpb_mc_avg2_var"] ) * tmpw_var_avg2 - self["tmpw" + "_mc_avg2_var"] = tmpw_var_avg2 + out["tmpw" + "_mc_avg2_var"] = tmpw_var_avg2 if conf_ints: # We first need to know the x-dim-chunk-size new_chunks_weighted = ((len(conf_ints),),) + (memchunk[1],) - avg_axis_avg2 = self["tmpw" + "_mc_avg2_set"].get_axis_num("mc") - q2 = self["tmpw" + "_mc_avg2_set"].data.map_blocks( + avg_axis_avg2 = mcparams["tmpw" + "_mc_avg2_set"].get_axis_num("mc") + q2 = mcparams["tmpw" + "_mc_avg2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avg2), chunks=new_chunks_weighted, # Explicitly define output chunks @@ -3027,17 +3039,17 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): new_axis=0, dtype=float, ) # The new CI dimension is added as firstax - self["tmpw" + "_mc_avg2"] = (("CI", x_dim2), q2) + out["tmpw" + "_mc_avg2"] = (("CI", x_dim2), q2) if ci_avg_x_flag1: - self["tmpw" + "_avgx1"] = self["tmpw" + "_avgsec"].mean(dim=x_dim2) + out["tmpw" + "_avgx1"] = out["tmpw" + "_avgsec"].mean(dim=x_dim2) - self["tmpw" + "_mc_avgx1_var"] = self["tmpw" + "_mc_set"].var(dim=x_dim2) + out["tmpw" + "_mc_avgx1_var"] = mcparams["tmpw" + "_mc_set"].var(dim=x_dim2) if conf_ints: new_chunks_weighted = ((len(conf_ints),),) + (memchunk[2],) - avg_axis = self["tmpw" + "_mc_set"].get_axis_num(["mc", x_dim2]) - q2 = self["tmpw" + "_mc_set"].data.map_blocks( + avg_axis = mcparams["tmpw" + "_mc_set"].get_axis_num(["mc", x_dim2]) + q2 = mcparams["tmpw" + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks_weighted, # Explicitly define output chunks @@ -3046,32 +3058,32 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): dtype=float, ) # The new CI dimension is added as # first axis - self["tmpw" + "_mc_avgx1"] = (("CI", time_dim2), q2) + out["tmpw" + "_mc_avgx1"] = (("CI", time_dim2), q2) if ci_avg_x_flag2: tmpw_var_avgx2 = 1 / ( - 1 / self["tmpf_mc_avgx2_var"] + 1 / self["tmpb_mc_avgx2_var"] + 1 / out["tmpf_mc_avgx2_var"] + 1 / out["tmpb_mc_avgx2_var"] ) q = ( - self["tmpf_mc_avgx2_set"] / self["tmpf_mc_avgx2_var"] - + self["tmpb_mc_avgx2_set"] / self["tmpb_mc_avgx2_var"] + mcparams["tmpf_mc_avgx2_set"] / out["tmpf_mc_avgx2_var"] + + mcparams["tmpb_mc_avgx2_set"] / out["tmpb_mc_avgx2_var"] ) * tmpw_var_avgx2 - self["tmpw" + "_mc_avgx2_set"] = q # + mcparams["tmpw" + "_mc_avgx2_set"] = q # - self["tmpw" + "_avgx2"] = ( - self["tmpf_avgx2"] / self["tmpf_mc_avgx2_var"] - + self["tmpb_avgx2"] / self["tmpb_mc_avgx2_var"] + out["tmpw" + "_avgx2"] = ( + out["tmpf_avgx2"] / out["tmpf_mc_avgx2_var"] + + out["tmpb_avgx2"] / out["tmpb_mc_avgx2_var"] ) * tmpw_var_avgx2 - self["tmpw" + "_mc_avgx2_var"] = tmpw_var_avgx2 + out["tmpw" + "_mc_avgx2_var"] = tmpw_var_avgx2 if conf_ints: # We first need to know the x-dim-chunk-size new_chunks_weighted = ((len(conf_ints),),) + (memchunk[2],) - avg_axis_avgx2 = self["tmpw" + "_mc_avgx2_set"].get_axis_num("mc") - q2 = self["tmpw" + "_mc_avgx2_set"].data.map_blocks( + avg_axis_avgx2 = mcparams["tmpw" + "_mc_avgx2_set"].get_axis_num("mc") + q2 = mcparams["tmpw" + "_mc_avgx2_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis_avgx2), chunks=new_chunks_weighted, # Explicitly define output chunks @@ -3079,7 +3091,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): new_axis=0, dtype=float, ) # The new CI dimension is added as firstax - self["tmpw" + "_mc_avgx2"] = (("CI", time_dim2), q2) + out["tmpw" + "_mc_avgx2"] = (("CI", time_dim2), q2) # Clean up the garbage. All arrays with a Monte Carlo dimension. if mc_remove_set_flag: @@ -3104,13 +3116,16 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): remove_mc_set.append(i + "_mc_avgx2_set") remove_mc_set.append(i + "_mc_avgsec_var") - if self.trans_att.size: + if "trans_att" in mcparams and mcparams.trans_att.size: remove_mc_set.append('talpha"_fw_mc') remove_mc_set.append('talpha"_bw_mc') for k in remove_mc_set: - if k in self: - del self[k] + if k in out: + print(f"Removed from results: {k}") + del out[k] + + self.update(out) pass def conf_int_single_ended( @@ -3194,6 +3209,9 @@ def conf_int_single_ended( """ check_deprecated_kwargs(kwargs) + out = xr.Dataset() + params = xr.Dataset() + if da_random_state: state = da_random_state else: @@ -3202,6 +3220,7 @@ def conf_int_single_ended( no, nt = self.st.data.shape if "trans_att" in self.keys(): nta = self.trans_att.size + else: nta = 0 @@ -3219,10 +3238,13 @@ def conf_int_single_ended( else: raise Exception("The size of `p_val` is not what I expected") - self.coords["mc"] = range(mc_sample_size) + params.coords["mc"] = range(mc_sample_size) + params.coords["x"] = self.x + params.coords["time"] = self.time if conf_ints: - self.coords["CI"] = conf_ints + out.coords["CI"] = conf_ints + params.coords["CI"] = conf_ints # WLS if isinstance(p_cov, str): @@ -3232,20 +3254,20 @@ def conf_int_single_ended( p_mc = sst.multivariate_normal.rvs(mean=p_val, cov=p_cov, size=mc_sample_size) if fixed_alpha: - self["alpha_mc"] = (("mc", "x"), p_mc[:, 1 : no + 1]) - self["c_mc"] = (("mc", "time"), p_mc[:, 1 + no : 1 + no + nt]) + params["alpha_mc"] = (("mc", "x"), p_mc[:, 1 : no + 1]) + params["c_mc"] = (("mc", "time"), p_mc[:, 1 + no : 1 + no + nt]) else: - self["dalpha_mc"] = (("mc",), p_mc[:, 1]) - self["c_mc"] = (("mc", "time"), p_mc[:, 2 : nt + 2]) + params["dalpha_mc"] = (("mc",), p_mc[:, 1]) + params["c_mc"] = (("mc", "time"), p_mc[:, 2 : nt + 2]) - self["gamma_mc"] = (("mc",), p_mc[:, 0]) + params["gamma_mc"] = (("mc",), p_mc[:, 0]) if nta: - self["ta_mc"] = ( + params["ta_mc"] = ( ("mc", "trans_att", "time"), np.reshape(p_mc[:, -nt * nta :], (mc_sample_size, nta, nt)), ) - rsize = (self.mc.size, self.x.size, self.time.size) + rsize = (params.mc.size, params.x.size, params.time.size) if reduce_memory_usage: memchunk = da.ones( @@ -3292,7 +3314,7 @@ def conf_int_single_ended( else: st_vari_da = da.from_array(st_vari, chunks=memchunk[1:]) - self[k] = ( + params[k] = ( ("mc", "x", "time"), state.normal( loc=loc, # has chunks=memchunk[1:] @@ -3305,52 +3327,50 @@ def conf_int_single_ended( ta_arr = np.zeros((mc_sample_size, no, nt)) if nta: - for ii, ta in enumerate(self["ta_mc"]): + for ii, ta in enumerate(params["ta_mc"]): for tai, taxi in zip(ta.values, self.trans_att.values): ta_arr[ii, self.x.values >= taxi] = ( ta_arr[ii, self.x.values >= taxi] + tai ) - self["ta_mc_arr"] = (("mc", "x", "time"), ta_arr) + params["ta_mc_arr"] = (("mc", "x", "time"), ta_arr) if fixed_alpha: - self["tmpf_mc_set"] = ( - self["gamma_mc"] + params["tmpf_mc_set"] = ( + params["gamma_mc"] / ( ( - np.log(self["r_st"]) - - np.log(self["r_ast"]) - + (self["c_mc"] + self["ta_mc_arr"]) + np.log(params["r_st"]) + - np.log(params["r_ast"]) + + (params["c_mc"] + params["ta_mc_arr"]) ) - + self["alpha_mc"] + + params["alpha_mc"] ) - 273.15 ) else: - self["tmpf_mc_set"] = ( - self["gamma_mc"] + params["tmpf_mc_set"] = ( + params["gamma_mc"] / ( ( - np.log(self["r_st"]) - - np.log(self["r_ast"]) - + (self["c_mc"] + self["ta_mc_arr"]) + np.log(params["r_st"]) + - np.log(params["r_ast"]) + + (params["c_mc"] + params["ta_mc_arr"]) ) - + (self["dalpha_mc"] * self.x) + + (params["dalpha_mc"] * params.x) ) - 273.15 ) avg_dims = ["mc"] - - avg_axis = self["tmpf_mc_set"].get_axis_num(avg_dims) - - self["tmpf_mc_var"] = (self["tmpf_mc_set"] - self["tmpf"]).var( + avg_axis = params["tmpf_mc_set"].get_axis_num(avg_dims) + out["tmpf_mc_var"] = (params["tmpf_mc_set"] - self["tmpf"]).var( dim=avg_dims, ddof=1 ) if conf_ints: - new_chunks = ((len(conf_ints),),) + self["tmpf_mc_set"].chunks[1:] + new_chunks = ((len(conf_ints),),) + params["tmpf_mc_set"].chunks[1:] - qq = self["tmpf_mc_set"] + qq = params["tmpf_mc_set"] q = qq.data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), @@ -3359,25 +3379,13 @@ def conf_int_single_ended( new_axis=0, ) # The new CI dimension is added as first axis - self["tmpf_mc"] = (("CI", "x", "time"), q) + out["tmpf_mc"] = (("CI", "x", "time"), q) - if mc_remove_set_flag: - drop_var = [ - "gamma_mc", - "dalpha_mc", - "alpha_mc", - "c_mc", - "mc", - "r_st", - "r_ast", - "tmpf_mc_set", - "ta_mc_arr", - ] - for k in drop_var: - if k in self: - del self[k] + if not mc_remove_set_flag: + out.update(params) - pass + self.update(out) + return out def conf_int_double_ended( self, @@ -3549,6 +3557,9 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): check_deprecated_kwargs(kwargs) + out = xr.Dataset() + params = xr.Dataset() + if da_random_state: # In testing environments assert isinstance(da_random_state, da.random.RandomState) @@ -3578,9 +3589,13 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): (mc_sample_size, no, nt), chunks={0: -1, 1: "auto", 2: "auto"} ).chunks - self.coords["mc"] = range(mc_sample_size) + params.coords["mc"] = range(mc_sample_size) + params.coords["x"] = self.x + params.coords["time"] = self.time + if conf_ints: self.coords["CI"] = conf_ints + params.coords["CI"] = conf_ints assert isinstance(p_val, (str, np.ndarray, np.generic)) if isinstance(p_val, str): @@ -3600,10 +3615,10 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): d_bw = p_val[1 + nt : 2 * nt + 1] alpha = p_val[2 * nt + 1 : 2 * nt + 1 + no] - self["gamma_mc"] = (tuple(), gamma) - self["alpha_mc"] = (("x",), alpha) - self["df_mc"] = (("time",), d_fw) - self["db_mc"] = (("time",), d_bw) + params["gamma_mc"] = (tuple(), gamma) + params["alpha_mc"] = (("x",), alpha) + params["df_mc"] = (("time",), d_fw) + params["db_mc"] = (("time",), d_bw) if nta: ta = p_val[2 * nt + 1 + no :].reshape((nt, 2, nta), order="F") @@ -3611,19 +3626,19 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): ta_bw = ta[:, 1, :] ta_fw_arr = np.zeros((no, nt)) - for tai, taxi in zip(ta_fw.T, self.coords["trans_att"].values): - ta_fw_arr[self.x.values >= taxi] = ( - ta_fw_arr[self.x.values >= taxi] + tai + for tai, taxi in zip(ta_fw.T, params.coords["trans_att"].values): + ta_fw_arr[params.x.values >= taxi] = ( + ta_fw_arr[params.x.values >= taxi] + tai ) ta_bw_arr = np.zeros((no, nt)) - for tai, taxi in zip(ta_bw.T, self.coords["trans_att"].values): - ta_bw_arr[self.x.values < taxi] = ( - ta_bw_arr[self.x.values < taxi] + tai + for tai, taxi in zip(ta_bw.T, params.coords["trans_att"].values): + ta_bw_arr[params.x.values < taxi] = ( + ta_bw_arr[params.x.values < taxi] + tai ) - self["talpha_fw_mc"] = (("x", "time"), ta_fw_arr) - self["talpha_bw_mc"] = (("x", "time"), ta_bw_arr) + params["talpha_fw_mc"] = (("x", "time"), ta_fw_arr) + params["talpha_bw_mc"] = (("x", "time"), ta_bw_arr) elif isinstance(p_cov, bool) and p_cov: raise NotImplementedError("Not an implemented option. Check p_cov argument") @@ -3657,9 +3672,9 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): d_fw = po_mc[:, 1 : nt + 1] d_bw = po_mc[:, 1 + nt : 2 * nt + 1] - self["gamma_mc"] = (("mc",), gamma) - self["df_mc"] = (("mc", "time"), d_fw) - self["db_mc"] = (("mc", "time"), d_bw) + params["gamma_mc"] = (("mc",), gamma) + params["df_mc"] = (("mc", "time"), d_fw) + params["db_mc"] = (("mc", "time"), d_bw) # calculate alpha seperately alpha = np.zeros((mc_sample_size, no), dtype=float) @@ -3679,7 +3694,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): alpha[:, not_ix_sec] = not_alpha_mc - self["alpha_mc"] = (("mc", "x"), alpha) + params["alpha_mc"] = (("mc", "x"), alpha) if nta: ta = po_mc[:, 2 * nt + 1 + nx_sec :].reshape( @@ -3692,10 +3707,10 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): (mc_sample_size, no, nt), chunks=memchunk, dtype=float ) for tai, taxi in zip( - ta_fw.swapaxes(0, 2), self.coords["trans_att"].values + ta_fw.swapaxes(0, 2), params.coords["trans_att"].values ): # iterate over the splices - i_splice = sum(self.x.values < taxi) + i_splice = sum(params.x.values < taxi) mask = create_da_ta2(no, i_splice, direction="fw", chunks=memchunk) ta_fw_arr += mask * tai.T[:, None, :] @@ -3704,15 +3719,15 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): (mc_sample_size, no, nt), chunks=memchunk, dtype=float ) for tai, taxi in zip( - ta_bw.swapaxes(0, 2), self.coords["trans_att"].values + ta_bw.swapaxes(0, 2), params.coords["trans_att"].values ): - i_splice = sum(self.x.values < taxi) + i_splice = sum(params.x.values < taxi) mask = create_da_ta2(no, i_splice, direction="bw", chunks=memchunk) ta_bw_arr += mask * tai.T[:, None, :] - self["talpha_fw_mc"] = (("mc", "x", "time"), ta_fw_arr) - self["talpha_bw_mc"] = (("mc", "x", "time"), ta_bw_arr) + params["talpha_fw_mc"] = (("mc", "x", "time"), ta_fw_arr) + params["talpha_bw_mc"] = (("mc", "x", "time"), ta_bw_arr) # Draw from the normal distributions for the Stokes intensities for k, st_labeli, st_vari in zip( @@ -3752,7 +3767,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): else: st_vari_da = da.from_array(st_vari, chunks=memchunk[1:]) - self[k] = ( + params[k] = ( ("mc", "x", "time"), state.normal( loc=loc, # has chunks=memchunk[1:] @@ -3766,45 +3781,45 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): if "tmpw" or label: if label == "tmpf": if nta: - self["tmpf_mc_set"] = ( - self["gamma_mc"] + params["tmpf_mc_set"] = ( + params["gamma_mc"] / ( - np.log(self["r_st"] / self["r_ast"]) - + self["df_mc"] - + self["alpha_mc"] - + self["talpha_fw_mc"] + np.log(params["r_st"] / params["r_ast"]) + + params["df_mc"] + + params["alpha_mc"] + + params["talpha_fw_mc"] ) - 273.15 ) else: - self["tmpf_mc_set"] = ( - self["gamma_mc"] + params["tmpf_mc_set"] = ( + params["gamma_mc"] / ( - np.log(self["r_st"] / self["r_ast"]) - + self["df_mc"] - + self["alpha_mc"] + np.log(params["r_st"] / params["r_ast"]) + + params["df_mc"] + + params["alpha_mc"] ) - 273.15 ) else: if nta: - self["tmpb_mc_set"] = ( - self["gamma_mc"] + params["tmpb_mc_set"] = ( + params["gamma_mc"] / ( - np.log(self["r_rst"] / self["r_rast"]) - + self["db_mc"] - - self["alpha_mc"] - + self["talpha_bw_mc"] + np.log(params["r_rst"] / params["r_rast"]) + + params["db_mc"] + - params["alpha_mc"] + + params["talpha_bw_mc"] ) - 273.15 ) else: - self["tmpb_mc_set"] = ( - self["gamma_mc"] + params["tmpb_mc_set"] = ( + params["gamma_mc"] / ( - np.log(self["r_rst"] / self["r_rast"]) - + self["db_mc"] - - self["alpha_mc"] + np.log(params["r_rst"] / params["r_rast"]) + + params["db_mc"] + - params["alpha_mc"] ) - 273.15 ) @@ -3814,19 +3829,21 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): xi = self.ufunc_per_section( sections=sections, x_indices=True, calc_per="all" ) - x_mask_ = [True if ix in xi else False for ix in range(self.x.size)] + x_mask_ = [ + True if ix in xi else False for ix in range(params.x.size) + ] x_mask = np.reshape(x_mask_, (1, -1, 1)) - self[label + "_mc_set"] = self[label + "_mc_set"].where(x_mask) + params[label + "_mc_set"] = params[label + "_mc_set"].where(x_mask) # subtract the mean temperature - q = self[label + "_mc_set"] - self[label] - self[label + "_mc_var"] = q.var(dim="mc", ddof=1) + q = params[label + "_mc_set"] - self[label] + out[label + "_mc_var"] = q.var(dim="mc", ddof=1) if conf_ints: - new_chunks = list(self[label + "_mc_set"].chunks) + new_chunks = list(params[label + "_mc_set"].chunks) new_chunks[0] = (len(conf_ints),) - avg_axis = self[label + "_mc_set"].get_axis_num("mc") - q = self[label + "_mc_set"].data.map_blocks( + avg_axis = params[label + "_mc_set"].get_axis_num("mc") + q = params[label + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks, # drop_axis=avg_axis, @@ -3834,37 +3851,37 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): new_axis=0, ) # The new CI dimension is added as firsaxis - self[label + "_mc"] = (("CI", "x", "time"), q) + out[label + "_mc"] = (("CI", "x", "time"), q) # Weighted mean of the forward and backward - tmpw_var = 1 / (1 / self["tmpf_mc_var"] + 1 / self["tmpb_mc_var"]) + tmpw_var = 1 / (1 / out["tmpf_mc_var"] + 1 / out["tmpb_mc_var"]) q = ( - self["tmpf_mc_set"] / self["tmpf_mc_var"] - + self["tmpb_mc_set"] / self["tmpb_mc_var"] + params["tmpf_mc_set"] / out["tmpf_mc_var"] + + params["tmpb_mc_set"] / out["tmpb_mc_var"] ) * tmpw_var - self["tmpw" + "_mc_set"] = q # + params["tmpw" + "_mc_set"] = q # - self["tmpw"] = ( - self["tmpf"] / self["tmpf_mc_var"] + self["tmpb"] / self["tmpb_mc_var"] + out["tmpw"] = ( + self["tmpf"] / out["tmpf_mc_var"] + self["tmpb"] / out["tmpb_mc_var"] ) * tmpw_var - q = self["tmpw" + "_mc_set"] - self["tmpw"] - self["tmpw" + "_mc_var"] = q.var(dim="mc", ddof=1) + q = params["tmpw" + "_mc_set"] - self["tmpw"] + out["tmpw" + "_mc_var"] = q.var(dim="mc", ddof=1) # Calculate the CI of the weighted MC_set if conf_ints: new_chunks_weighted = ((len(conf_ints),),) + memchunk[1:] - avg_axis = self["tmpw" + "_mc_set"].get_axis_num("mc") - q2 = self["tmpw" + "_mc_set"].data.map_blocks( + avg_axis = params["tmpw" + "_mc_set"].get_axis_num("mc") + q2 = params["tmpw" + "_mc_set"].data.map_blocks( lambda x: np.percentile(x, q=conf_ints, axis=avg_axis), chunks=new_chunks_weighted, # Explicitly define output chunks drop_axis=avg_axis, # avg dimensions are dropped new_axis=0, dtype=float, ) # The new CI dimension is added as first axis - self["tmpw" + "_mc"] = (("CI", "x", "time"), q2) + out["tmpw" + "_mc"] = (("CI", "x", "time"), q2) # Clean up the garbage. All arrays with a Monte Carlo dimension. if mc_remove_set_flag: @@ -3887,9 +3904,14 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): remove_mc_set.append('talpha"_bw_mc') for k in remove_mc_set: - if k in self: - del self[k] - pass + if k in out: + del out[k] + + if not mc_remove_set_flag: + out.update(params) + + self.update(out) + return out def in_confidence_interval(self, ci_label, conf_ints=None, sections=None): """ diff --git a/tests/test_dtscalibration.py b/tests/test_dtscalibration.py index 23516b09..b354c584 100644 --- a/tests/test_dtscalibration.py +++ b/tests/test_dtscalibration.py @@ -130,7 +130,11 @@ def test_variance_input_types_single(): ) ds.conf_int_single_ended( - st_var=st_var, ast_var=st_var, mc_sample_size=100, da_random_state=state + st_var=st_var, + ast_var=st_var, + mc_sample_size=100, + da_random_state=state, + mc_remove_set_flag=False, ) assert_almost_equal_verbose( From a515f8b3e1e73752a1181278f4f015b49b67feab Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 3 Sep 2023 15:57:19 +0200 Subject: [PATCH 18/22] Updated sections example notebook --- docs/notebooks/03Define_sections.ipynb | 35 +++++--------------------- 1 file changed, 6 insertions(+), 29 deletions(-) diff --git a/docs/notebooks/03Define_sections.ipynb b/docs/notebooks/03Define_sections.ipynb index ef9d97b2..8b1000ea 100644 --- a/docs/notebooks/03Define_sections.ipynb +++ b/docs/notebooks/03Define_sections.ipynb @@ -23,7 +23,7 @@ "source": [ "import os\n", "\n", - "from dtscalibration import read_silixa_files" + "from dtscalibration import read_silixa_files\n" ] }, { @@ -40,7 +40,7 @@ "outputs": [], "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"double_ended2\")\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" ] }, { @@ -66,7 +66,8 @@ "outputs": [], "source": [ "print(ds.timeseries_keys) # list the available timeseeries\n", - "ds.probe1Temperature.plot(figsize=(12, 8)); # plot one of the timeseries" + "ds.probe1Temperature.plot(figsize=(12, 8))\n", + "# plot one of the timeseries\n" ] }, { @@ -115,24 +116,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n", - "ds.sections = sections" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "execution": { - "iopub.execute_input": "2022-04-06T08:09:10.309314Z", - "iopub.status.busy": "2022-04-06T08:09:10.308985Z", - "iopub.status.idle": "2022-04-06T08:09:10.323444Z", - "shell.execute_reply": "2022-04-06T08:09:10.322874Z" - } - }, - "outputs": [], - "source": [ - "ds.sections" + "}\n" ] }, { @@ -141,13 +125,6 @@ "source": [ "NetCDF files do not support reading/writing python dictionaries. Internally the sections dictionary is stored in `ds._sections` as a string encoded with yaml, which can be saved to a netCDF file. Each time the sections dictionary is requested, yaml decodes the string and evaluates it to the Python dictionary. " ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -166,7 +143,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.11" + "version": "3.10.10" } }, "nbformat": 4, From 39f71c9c5384b182767f1f12630683a49c1cdc71 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 3 Sep 2023 16:26:53 +0200 Subject: [PATCH 19/22] Updated example notebooks --- .../04Calculate_variance_Stokes.ipynb | 13 ++++--- docs/notebooks/07Calibrate_single_ended.ipynb | 20 +++++------ docs/notebooks/08Calibrate_double_ended.ipynb | 36 ++++++++++--------- .../12Datastore_from_numpy_arrays.ipynb | 23 ++++++------ .../13Fixed_parameter_calibration.ipynb | 19 +++++----- docs/notebooks/14Lossy_splices.ipynb | 29 +++++++-------- docs/notebooks/15Matching_sections.ipynb | 32 +++++++++-------- docs/notebooks/16Averaging_temperatures.ipynb | 26 ++++++++------ ...Temperature_uncertainty_single_ended.ipynb | 18 +++++----- src/dtscalibration/datastore.py | 1 + 10 files changed, 115 insertions(+), 102 deletions(-) diff --git a/docs/notebooks/04Calculate_variance_Stokes.ipynb b/docs/notebooks/04Calculate_variance_Stokes.ipynb index 2e0eba4b..0de14167 100644 --- a/docs/notebooks/04Calculate_variance_Stokes.ipynb +++ b/docs/notebooks/04Calculate_variance_Stokes.ipynb @@ -55,7 +55,7 @@ "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"double_ended2\")\n", "\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" ] }, { @@ -81,8 +81,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n", - "ds.sections = sections" + "}\n" ] }, { @@ -120,13 +119,13 @@ }, "outputs": [], "source": [ - "I_var, residuals = ds.variance_stokes_constant(st_label=\"st\")\n", + "I_var, residuals = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", "print(\n", " \"The variance of the Stokes signal along the reference sections \"\n", " \"is approximately {:.2f} on a {:.1f} sec acquisition time\".format(\n", " I_var, ds.userAcquisitionTimeFW.data[0]\n", " )\n", - ")" + ")\n" ] }, { @@ -153,7 +152,7 @@ " robust=True,\n", " units=\"\",\n", " method=\"single\",\n", - ")" + ")\n" ] }, { @@ -187,7 +186,7 @@ "x = np.linspace(mean - 3 * sigma, mean + 3 * sigma, 100)\n", "approximated_normal_fit = scipy.stats.norm.pdf(x, mean, sigma)\n", "residuals.plot.hist(bins=50, figsize=(12, 8), density=True)\n", - "plt.plot(x, approximated_normal_fit);" + "plt.plot(x, approximated_normal_fit)\n" ] }, { diff --git a/docs/notebooks/07Calibrate_single_ended.ipynb b/docs/notebooks/07Calibrate_single_ended.ipynb index 28616d9e..4c591ede 100644 --- a/docs/notebooks/07Calibrate_single_ended.ipynb +++ b/docs/notebooks/07Calibrate_single_ended.ipynb @@ -72,7 +72,7 @@ "outputs": [], "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"single_ended\")\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" ] }, { @@ -90,10 +90,10 @@ "outputs": [], "source": [ "ds = ds.sel(x=slice(-30, 101)) # dismiss parts of the fiber that are not interesting\n", - "ds.sections = {\n", + "sections = {\n", " \"probe1Temperature\": [slice(20, 25.5)], # warm bath\n", " \"probe2Temperature\": [slice(5.5, 15.5)], # cold bath\n", - "}" + "}\n" ] }, { @@ -119,8 +119,8 @@ }, "outputs": [], "source": [ - "st_var, resid = ds.variance_stokes_constant(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes_constant(st_label=\"ast\")" + "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n" ] }, { @@ -136,7 +136,7 @@ "metadata": {}, "outputs": [], "source": [ - "resid.plot(figsize=(12, 4));" + "resid.plot(figsize=(12, 4))\n" ] }, { @@ -160,7 +160,7 @@ }, "outputs": [], "source": [ - "ds.calibration_single_ended(st_var=st_var, ast_var=ast_var)" + "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)\n" ] }, { @@ -177,7 +177,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds.tmpf.plot(figsize=(12, 4));" + "ds.tmpf.plot(figsize=(12, 4))\n" ] }, { @@ -189,7 +189,7 @@ "ds1 = ds.isel(time=0)\n", "ds1.tmpf.plot(figsize=(12, 4))\n", "(ds1.tmpf_var**0.5).plot(figsize=(12, 4))\n", - "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\");" + "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n" ] }, { @@ -217,7 +217,7 @@ "outputs": [], "source": [ "ds1.st.plot(figsize=(12, 8))\n", - "ds1.ast.plot();" + "ds1.ast.plot()\n" ] }, { diff --git a/docs/notebooks/08Calibrate_double_ended.ipynb b/docs/notebooks/08Calibrate_double_ended.ipynb index 2ed0f35b..21ba4b45 100644 --- a/docs/notebooks/08Calibrate_double_ended.ipynb +++ b/docs/notebooks/08Calibrate_double_ended.ipynb @@ -70,7 +70,7 @@ "\n", "ds_notaligned = read_silixa_files(\n", " directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\"\n", - ")" + ")\n" ] }, { @@ -102,7 +102,7 @@ }, "outputs": [], "source": [ - "ds_notaligned = ds_notaligned.sel(x=slice(0, 100)) # only calibrate parts of the fiber" + "ds_notaligned = ds_notaligned.sel(x=slice(0, 100)) # only calibrate parts of the fiber\n" ] }, { @@ -124,7 +124,7 @@ " plot_result=True,\n", " figsize=(12, 6),\n", ")\n", - "print(suggested_shift)" + "print(suggested_shift)\n" ] }, { @@ -140,7 +140,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds = shift_double_ended(ds_notaligned, suggested_shift[0])" + "ds = shift_double_ended(ds_notaligned, suggested_shift[0])\n" ] }, { @@ -157,10 +157,10 @@ "metadata": {}, "outputs": [], "source": [ - "ds.sections = {\n", + "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}" + "}\n" ] }, { @@ -186,10 +186,10 @@ }, "outputs": [], "source": [ - "st_var, resid = ds.variance_stokes_constant(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes_constant(st_label=\"ast\")\n", - "rst_var, _ = ds.variance_stokes_constant(st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes_constant(st_label=\"rast\")" + "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"rast\")\n" ] }, { @@ -212,7 +212,7 @@ }, "outputs": [], "source": [ - "resid.plot(figsize=(12, 4));" + "resid.plot(figsize=(12, 4))\n" ] }, { @@ -237,11 +237,12 @@ "outputs": [], "source": [ "ds.calibration_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", " rast_var=rast_var,\n", - ")" + ")\n" ] }, { @@ -257,7 +258,7 @@ }, "outputs": [], "source": [ - "ds.tmpw.plot(figsize=(12, 4));" + "ds.tmpw.plot(figsize=(12, 4))\n" ] }, { @@ -295,7 +296,7 @@ "ds.tmpw_var.plot(figsize=(12, 4))\n", "ds1 = ds.isel(time=-1) # take only the first timestep\n", "(ds1.tmpw_var**0.5).plot(figsize=(12, 4))\n", - "plt.gca().set_ylabel(\"Standard error ($^\\circ$C)\");" + "plt.gca().set_ylabel(\"Standard error ($^\\circ$C)\")\n" ] }, { @@ -323,6 +324,7 @@ "outputs": [], "source": [ "ds.conf_int_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -331,7 +333,7 @@ " mc_sample_size=500,\n", ") # < increase sample size for better approximation\n", "\n", - "ds.tmpw_mc_var.plot(figsize=(12, 4));" + "ds.tmpw_mc_var.plot(figsize=(12, 4))\n" ] }, { @@ -351,7 +353,7 @@ "ds1.tmpw.plot(linewidth=0.7, figsize=(12, 4))\n", "ds1.tmpw_mc.isel(CI=0).plot(linewidth=0.7, label=\"CI: 2.5%\")\n", "ds1.tmpw_mc.isel(CI=1).plot(linewidth=0.7, label=\"CI: 97.5%\")\n", - "plt.legend(fontsize=\"small\")" + "plt.legend(fontsize=\"small\")\n" ] }, { @@ -380,7 +382,7 @@ "(ds1.tmpb_var**0.5).plot()\n", "(ds1.tmpw_var**0.5).plot()\n", "(ds1.tmpw_mc_var**0.5).plot()\n", - "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\");" + "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n" ] }, { diff --git a/docs/notebooks/12Datastore_from_numpy_arrays.ipynb b/docs/notebooks/12Datastore_from_numpy_arrays.ipynb index 2ddb9a20..1f3d3ce7 100644 --- a/docs/notebooks/12Datastore_from_numpy_arrays.ipynb +++ b/docs/notebooks/12Datastore_from_numpy_arrays.ipynb @@ -26,7 +26,7 @@ "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "\n", - "from dtscalibration import DataStore, read_silixa_files" + "from dtscalibration import DataStore, read_silixa_files\n" ] }, { @@ -61,7 +61,7 @@ "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"single_ended\")\n", "\n", - "ds_silixa = read_silixa_files(directory=filepath, silent=True)" + "ds_silixa = read_silixa_files(directory=filepath, silent=True)\n" ] }, { @@ -89,7 +89,7 @@ "x = ds_silixa.x.values\n", "time = ds_silixa.time.values\n", "ST = ds_silixa.st.values\n", - "AST = ds_silixa.ast.values" + "AST = ds_silixa.ast.values\n" ] }, { @@ -116,7 +116,7 @@ "ds[\"x\"] = (\"x\", x)\n", "ds[\"time\"] = (\"time\", time)\n", "ds[\"st\"] = ([\"x\", \"time\"], ST)\n", - "ds[\"ast\"] = ([\"x\", \"time\"], AST)" + "ds[\"ast\"] = ([\"x\", \"time\"], AST)\n" ] }, { @@ -133,7 +133,7 @@ "outputs": [], "source": [ "ds = DataStore(ds)\n", - "print(ds)" + "print(ds)\n" ] }, { @@ -169,7 +169,7 @@ "ds[\"temp1\"] = ds_silixa[\"probe1Temperature\"]\n", "ds[\"temp2\"] = ds_silixa[\"probe2Temperature\"]\n", "\n", - "ds.attrs[\"isDoubleEnded\"] = \"0\"" + "ds.attrs[\"isDoubleEnded\"] = \"0\"\n" ] }, { @@ -197,13 +197,12 @@ " \"temp1\": [slice(20, 25.5)], # warm bath\n", " \"temp2\": [slice(5.5, 15.5)], # cold bath\n", "}\n", - "ds.sections = sections\n", "\n", - "st_var, resid = ds.variance_stokes_constant(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes_constant(st_label=\"ast\")\n", - "ds.calibration_single_ended(st_var=st_var, ast_var=ast_var)\n", + "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", + "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)\n", "\n", - "ds.isel(time=0).tmpf.plot();" + "ds.isel(time=0).tmpf.plot()\n" ] }, { @@ -212,7 +211,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds" + "ds\n" ] }, { diff --git a/docs/notebooks/13Fixed_parameter_calibration.ipynb b/docs/notebooks/13Fixed_parameter_calibration.ipynb index c5b9fb47..cf5a78c1 100644 --- a/docs/notebooks/13Fixed_parameter_calibration.ipynb +++ b/docs/notebooks/13Fixed_parameter_calibration.ipynb @@ -39,8 +39,7 @@ " \"probe1Temperature\": [\n", " slice(20, 25.5)\n", " ], # we only use the warm bath in this notebook\n", - "}\n", - "ds100.sections = sections" + "}" ] }, { @@ -69,11 +68,15 @@ "fix_gamma = (481.9, 0) # (gamma value, gamma variance)\n", "fix_dalpha = (-2.014e-5, 0) # (alpha value, alpha variance)\n", "\n", - "st_var, resid = ds100.variance_stokes_constant(st_label=\"st\")\n", - "ast_var, _ = ds100.variance_stokes_constant(st_label=\"ast\")\n", + "st_var, resid = ds100.variance_stokes_constant(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds100.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", "ds100.calibration_single_ended(\n", - " st_var=st_var, ast_var=ast_var, fix_gamma=fix_gamma, fix_dalpha=fix_dalpha\n", - ")" + " sections=sections,\n", + " st_var=st_var,\n", + " ast_var=ast_var,\n", + " fix_gamma=fix_gamma,\n", + " fix_dalpha=fix_dalpha,\n", + ")\n" ] }, { @@ -97,7 +100,7 @@ "outputs": [], "source": [ "print(\"gamma used in calibration:\", ds100.gamma.values)\n", - "print(\"dalpha used in calibration:\", ds100.dalpha.values)" + "print(\"dalpha used in calibration:\", ds100.dalpha.values)\n" ] }, { @@ -129,7 +132,7 @@ " linewidth=1, label=\"Device calibrated\"\n", ") # plot the temperature calibrated by the device\n", "plt.title(\"Temperature at the first time step\")\n", - "plt.legend();" + "plt.legend()\n" ] }, { diff --git a/docs/notebooks/14Lossy_splices.ipynb b/docs/notebooks/14Lossy_splices.ipynb index c2c98e6e..5e41b774 100644 --- a/docs/notebooks/14Lossy_splices.ipynb +++ b/docs/notebooks/14Lossy_splices.ipynb @@ -72,8 +72,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n", - "ds.sections = sections" + "}\n" ] }, { @@ -102,7 +101,7 @@ "ds[\"ast\"] = ds.ast.where(ds.x < 50, ds.ast * 0.82)\n", "\n", "ds[\"rst\"] = ds.rst.where(ds.x > 50, ds.rst * 0.85)\n", - "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)" + "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)\n" ] }, { @@ -122,7 +121,7 @@ "ds.isel(time=0).ast.plot(label=\"ast\")\n", "ds.isel(time=0).rst.plot(label=\"rst\")\n", "ds.isel(time=0).rast.plot(label=\"rast\")\n", - "plt.legend();" + "plt.legend()\n" ] }, { @@ -147,19 +146,20 @@ "source": [ "ds_a = ds.copy(deep=True)\n", "\n", - "st_var, resid = ds_a.variance_stokes(st_label=\"st\")\n", - "ast_var, _ = ds_a.variance_stokes(st_label=\"ast\")\n", - "rst_var, _ = ds_a.variance_stokes(st_label=\"rst\")\n", - "rast_var, _ = ds_a.variance_stokes(st_label=\"rast\")\n", + "st_var, resid = ds_a.variance_stokes(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"rast\")\n", "\n", "ds_a.calibration_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", " rast_var=rast_var,\n", ")\n", "\n", - "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\");" + "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")\n" ] }, { @@ -184,12 +184,13 @@ }, "outputs": [], "source": [ - "st_var, resid = ds.variance_stokes(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes(st_label=\"ast\")\n", - "rst_var, _ = ds.variance_stokes(st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes(st_label=\"rast\")\n", + "st_var, resid = ds.variance_stokes(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds.variance_stokes(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds.variance_stokes(sections=sections, st_label=\"rast\")\n", "\n", "ds.calibration_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -199,7 +200,7 @@ "\n", "ds_a.isel(time=0).tmpw.plot(label=\"no trans. att.\")\n", "ds.isel(time=0).tmpw.plot(label=\"with trans. att.\")\n", - "plt.legend();" + "plt.legend()\n" ] }, { diff --git a/docs/notebooks/15Matching_sections.ipynb b/docs/notebooks/15Matching_sections.ipynb index a83c3730..cb709b50 100644 --- a/docs/notebooks/15Matching_sections.ipynb +++ b/docs/notebooks/15Matching_sections.ipynb @@ -68,8 +68,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0)], # warm bath\n", - "}\n", - "ds.sections = sections" + "}\n" ] }, { @@ -96,7 +95,7 @@ "ds[\"ast\"] = ds.ast.where(ds.x < 50, ds.ast * 0.82)\n", "\n", "ds[\"rst\"] = ds.rst.where(ds.x > 50, ds.rst * 0.85)\n", - "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)" + "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)\n" ] }, { @@ -123,16 +122,20 @@ "source": [ "ds_a = ds.copy(deep=True)\n", "\n", - "st_var, resid = ds_a.variance_stokes(st_label=\"st\")\n", - "ast_var, _ = ds_a.variance_stokes(st_label=\"ast\")\n", - "rst_var, _ = ds_a.variance_stokes(st_label=\"rst\")\n", - "rast_var, _ = ds_a.variance_stokes(st_label=\"rast\")\n", + "st_var, resid = ds_a.variance_stokes(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds_a.variance_stokes(sections=sections, st_label=\"rast\")\n", "\n", "ds_a.calibration_double_ended(\n", - " st_var=st_var, ast_var=ast_var, rst_var=rst_var, rast_var=rast_var\n", + " sections=sections,\n", + " st_var=st_var,\n", + " ast_var=ast_var,\n", + " rst_var=rst_var,\n", + " rast_var=rast_var,\n", ")\n", "\n", - "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")" + "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")\n" ] }, { @@ -165,12 +168,13 @@ "source": [ "matching_sections = [(slice(7.5, 17.6), slice(69, 79.1), False)]\n", "\n", - "st_var, resid = ds.variance_stokes(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes(st_label=\"ast\")\n", - "rst_var, _ = ds.variance_stokes(st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes(st_label=\"rast\")\n", + "st_var, resid = ds.variance_stokes(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds.variance_stokes(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds.variance_stokes(sections=sections, st_label=\"rast\")\n", "\n", "ds.calibration_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -181,7 +185,7 @@ "\n", "ds_a.isel(time=0).tmpw.plot(label=\"normal calibration\")\n", "ds.isel(time=0).tmpw.plot(label=\"matching sections\")\n", - "plt.legend()" + "plt.legend()\n" ] }, { diff --git a/docs/notebooks/16Averaging_temperatures.ipynb b/docs/notebooks/16Averaging_temperatures.ipynb index 6079bd9f..76925061 100644 --- a/docs/notebooks/16Averaging_temperatures.ipynb +++ b/docs/notebooks/16Averaging_temperatures.ipynb @@ -63,8 +63,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n", - "ds.sections = sections" + "}\n" ] }, { @@ -80,10 +79,10 @@ }, "outputs": [], "source": [ - "st_var, resid = ds.variance_stokes(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes(st_label=\"ast\")\n", - "rst_var, _ = ds.variance_stokes(st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes(st_label=\"rast\")" + "st_var, resid = ds.variance_stokes(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes(sections=sections, st_label=\"ast\")\n", + "rst_var, _ = ds.variance_stokes(sections=sections, st_label=\"rst\")\n", + "rast_var, _ = ds.variance_stokes(sections=sections, st_label=\"rast\")\n" ] }, { @@ -100,11 +99,12 @@ "outputs": [], "source": [ "ds.calibration_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", " rast_var=rast_var,\n", - ")" + ")\n" ] }, { @@ -179,6 +179,7 @@ "outputs": [], "source": [ "ds.average_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -190,7 +191,7 @@ " ci_avg_time_isel=[0, 1, 2, 3, 4, 5],\n", " ci_avg_time_sel=None,\n", ")\n", - "ds.tmpw_mc_avg1.plot(hue=\"CI\", linewidth=0.8);" + "ds.tmpw_mc_avg1.plot(hue=\"CI\", linewidth=0.8)\n" ] }, { @@ -239,6 +240,7 @@ "outputs": [], "source": [ "ds.average_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -250,7 +252,7 @@ " ci_avg_time_isel=[0, 1, 2, 3, 4, 5],\n", " ci_avg_time_sel=None,\n", ")\n", - "ds.tmpw_mc_avg2.plot(hue=\"CI\", linewidth=0.8);" + "ds.tmpw_mc_avg2.plot(hue=\"CI\", linewidth=0.8)\n" ] }, { @@ -299,6 +301,7 @@ "outputs": [], "source": [ "ds.average_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -310,7 +313,7 @@ " ci_avg_x_sel=slice(7.5, 17.0),\n", " ci_avg_x_isel=None,\n", ")\n", - "ds.tmpw_mc_avgx1.plot(hue=\"CI\", linewidth=0.8);" + "ds.tmpw_mc_avgx1.plot(hue=\"CI\", linewidth=0.8)\n" ] }, { @@ -359,6 +362,7 @@ "outputs": [], "source": [ "ds.average_double_ended(\n", + " sections=sections,\n", " st_var=st_var,\n", " ast_var=ast_var,\n", " rst_var=rst_var,\n", @@ -370,7 +374,7 @@ " ci_avg_x_sel=slice(7.5, 17.0),\n", " ci_avg_x_isel=None,\n", ")\n", - "ds.tmpw_mc_avgx2.plot(hue=\"CI\", linewidth=0.8);" + "ds.tmpw_mc_avgx2.plot(hue=\"CI\", linewidth=0.8)\n" ] }, { diff --git a/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb b/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb index f270404f..478a490e 100644 --- a/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb +++ b/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb @@ -32,15 +32,15 @@ "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n", "\n", "ds = ds.sel(x=slice(-30, 101)) # dismiss parts of the fiber that are not interesting\n", - "ds.sections = {\n", + "sections = {\n", " \"probe1Temperature\": [slice(20, 25.5)], # warm bath\n", " \"probe2Temperature\": [slice(5.5, 15.5)], # cold bath\n", "}\n", "\n", - "st_var, resid = ds.variance_stokes_constant(st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes_constant(st_label=\"ast\")\n", + "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", + "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", "\n", - "ds.calibration_single_ended(st_var=st_var, ast_var=ast_var)" + "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)" ] }, { @@ -73,7 +73,7 @@ "plt.fill_between(ds1.x, y1=ds1.tmpf_var, y2=stast_var, label=\"Parameter estimation\")\n", "plt.suptitle(\"Variance of the estimated temperature\")\n", "plt.ylabel(\"$\\sigma^2$ ($^\\circ$C$^2$)\")\n", - "plt.legend(fontsize=\"small\");" + "plt.legend(fontsize=\"small\")\n" ] }, { @@ -84,7 +84,7 @@ "outputs": [], "source": [ "# The effects of the parameter uncertainty can be further inspected\n", - "ds1.var_fw_da.plot(hue=\"comp_fw\", figsize=(12, 4));" + "ds1.var_fw_da.plot(hue=\"comp_fw\", figsize=(12, 4))\n" ] }, { @@ -121,7 +121,7 @@ "source": [ "ds.conf_int_single_ended(\n", " st_var=st_var, ast_var=ast_var, conf_ints=[2.5, 97.5], mc_sample_size=500\n", - ")" + ")\n" ] }, { @@ -145,7 +145,7 @@ "(ds1.tmpf_mc_var**0.5).plot(figsize=(12, 4), label=\"Monte Carlo approx.\")\n", "(ds1.tmpf_var**0.5).plot(label=\"Linear error approx.\")\n", "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n", - "plt.legend(fontsize=\"small\");" + "plt.legend(fontsize=\"small\")\n" ] }, { @@ -166,7 +166,7 @@ "ds1.tmpf.plot(linewidth=0.7, figsize=(12, 4))\n", "ds1.tmpf_mc.sel(CI=2.5).plot(linewidth=0.7, label=\"CI: 2.5%\")\n", "ds1.tmpf_mc.sel(CI=97.5).plot(linewidth=0.7, label=\"CI: 97.5%\")\n", - "plt.legend(fontsize=\"small\");" + "plt.legend(fontsize=\"small\")\n" ] }, { diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index db555798..459a45dd 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -3649,6 +3649,7 @@ def create_da_ta2(no, i_splice, direction="fw", chunks=None): p_cov = self[p_cov].values assert p_cov.shape == (npar, npar) + assert sections is not None, "Define sections" ix_sec = self.ufunc_per_section( sections=sections, x_indices=True, calc_per="all" ) From a40c810b68bb4971cbc1c7e8f738cff94bfdf0f7 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 3 Sep 2023 16:31:26 +0200 Subject: [PATCH 20/22] Add api to docs --- docs/api/dtscalibration.DataStore.rst | 7 +++++++ docs/api/dtscalibration.check_dims.rst | 6 ++++++ docs/api/dtscalibration.check_timestep_allclose.rst | 6 ++++++ docs/api/dtscalibration.get_netcdf_encoding.rst | 6 ++++++ docs/api/dtscalibration.merge_double_ended.rst | 6 ++++++ docs/api/dtscalibration.open_datastore.rst | 6 ++++++ docs/api/dtscalibration.open_mf_datastore.rst | 6 ++++++ docs/api/dtscalibration.plot_accuracy.rst | 6 ++++++ ...dtscalibration.plot_location_residuals_double_ended.rst | 6 ++++++ .../dtscalibration.plot_residuals_reference_sections.rst | 6 ++++++ ...alibration.plot_residuals_reference_sections_single.rst | 6 ++++++ docs/api/dtscalibration.plot_sigma_report.rst | 6 ++++++ docs/api/dtscalibration.read_apsensing_files.rst | 6 ++++++ docs/api/dtscalibration.read_sensornet_files.rst | 6 ++++++ docs/api/dtscalibration.read_sensortran_files.rst | 6 ++++++ docs/api/dtscalibration.read_silixa_files.rst | 6 ++++++ docs/api/dtscalibration.shift_double_ended.rst | 6 ++++++ .../dtscalibration.suggest_cable_shift_double_ended.rst | 6 ++++++ 18 files changed, 109 insertions(+) create mode 100644 docs/api/dtscalibration.DataStore.rst create mode 100644 docs/api/dtscalibration.check_dims.rst create mode 100644 docs/api/dtscalibration.check_timestep_allclose.rst create mode 100644 docs/api/dtscalibration.get_netcdf_encoding.rst create mode 100644 docs/api/dtscalibration.merge_double_ended.rst create mode 100644 docs/api/dtscalibration.open_datastore.rst create mode 100644 docs/api/dtscalibration.open_mf_datastore.rst create mode 100644 docs/api/dtscalibration.plot_accuracy.rst create mode 100644 docs/api/dtscalibration.plot_location_residuals_double_ended.rst create mode 100644 docs/api/dtscalibration.plot_residuals_reference_sections.rst create mode 100644 docs/api/dtscalibration.plot_residuals_reference_sections_single.rst create mode 100644 docs/api/dtscalibration.plot_sigma_report.rst create mode 100644 docs/api/dtscalibration.read_apsensing_files.rst create mode 100644 docs/api/dtscalibration.read_sensornet_files.rst create mode 100644 docs/api/dtscalibration.read_sensortran_files.rst create mode 100644 docs/api/dtscalibration.read_silixa_files.rst create mode 100644 docs/api/dtscalibration.shift_double_ended.rst create mode 100644 docs/api/dtscalibration.suggest_cable_shift_double_ended.rst diff --git a/docs/api/dtscalibration.DataStore.rst b/docs/api/dtscalibration.DataStore.rst new file mode 100644 index 00000000..8ed756d9 --- /dev/null +++ b/docs/api/dtscalibration.DataStore.rst @@ -0,0 +1,7 @@ +DataStore +========= + +.. currentmodule:: dtscalibration + +.. autoclass:: DataStore + :show-inheritance: diff --git a/docs/api/dtscalibration.check_dims.rst b/docs/api/dtscalibration.check_dims.rst new file mode 100644 index 00000000..1dc68e7e --- /dev/null +++ b/docs/api/dtscalibration.check_dims.rst @@ -0,0 +1,6 @@ +check_dims +========== + +.. currentmodule:: dtscalibration + +.. autofunction:: check_dims diff --git a/docs/api/dtscalibration.check_timestep_allclose.rst b/docs/api/dtscalibration.check_timestep_allclose.rst new file mode 100644 index 00000000..d655bfe4 --- /dev/null +++ b/docs/api/dtscalibration.check_timestep_allclose.rst @@ -0,0 +1,6 @@ +check_timestep_allclose +======================= + +.. currentmodule:: dtscalibration + +.. autofunction:: check_timestep_allclose diff --git a/docs/api/dtscalibration.get_netcdf_encoding.rst b/docs/api/dtscalibration.get_netcdf_encoding.rst new file mode 100644 index 00000000..820b0a0f --- /dev/null +++ b/docs/api/dtscalibration.get_netcdf_encoding.rst @@ -0,0 +1,6 @@ +get_netcdf_encoding +=================== + +.. currentmodule:: dtscalibration + +.. autofunction:: get_netcdf_encoding diff --git a/docs/api/dtscalibration.merge_double_ended.rst b/docs/api/dtscalibration.merge_double_ended.rst new file mode 100644 index 00000000..c8e512b5 --- /dev/null +++ b/docs/api/dtscalibration.merge_double_ended.rst @@ -0,0 +1,6 @@ +merge_double_ended +================== + +.. currentmodule:: dtscalibration + +.. autofunction:: merge_double_ended diff --git a/docs/api/dtscalibration.open_datastore.rst b/docs/api/dtscalibration.open_datastore.rst new file mode 100644 index 00000000..95a987bd --- /dev/null +++ b/docs/api/dtscalibration.open_datastore.rst @@ -0,0 +1,6 @@ +open_datastore +============== + +.. currentmodule:: dtscalibration + +.. autofunction:: open_datastore diff --git a/docs/api/dtscalibration.open_mf_datastore.rst b/docs/api/dtscalibration.open_mf_datastore.rst new file mode 100644 index 00000000..edab4ebe --- /dev/null +++ b/docs/api/dtscalibration.open_mf_datastore.rst @@ -0,0 +1,6 @@ +open_mf_datastore +================= + +.. currentmodule:: dtscalibration + +.. autofunction:: open_mf_datastore diff --git a/docs/api/dtscalibration.plot_accuracy.rst b/docs/api/dtscalibration.plot_accuracy.rst new file mode 100644 index 00000000..a2e41fc2 --- /dev/null +++ b/docs/api/dtscalibration.plot_accuracy.rst @@ -0,0 +1,6 @@ +plot_accuracy +============= + +.. currentmodule:: dtscalibration + +.. autofunction:: plot_accuracy diff --git a/docs/api/dtscalibration.plot_location_residuals_double_ended.rst b/docs/api/dtscalibration.plot_location_residuals_double_ended.rst new file mode 100644 index 00000000..ad0d16db --- /dev/null +++ b/docs/api/dtscalibration.plot_location_residuals_double_ended.rst @@ -0,0 +1,6 @@ +plot_location_residuals_double_ended +==================================== + +.. currentmodule:: dtscalibration + +.. autofunction:: plot_location_residuals_double_ended diff --git a/docs/api/dtscalibration.plot_residuals_reference_sections.rst b/docs/api/dtscalibration.plot_residuals_reference_sections.rst new file mode 100644 index 00000000..45f2529b --- /dev/null +++ b/docs/api/dtscalibration.plot_residuals_reference_sections.rst @@ -0,0 +1,6 @@ +plot_residuals_reference_sections +================================= + +.. currentmodule:: dtscalibration + +.. autofunction:: plot_residuals_reference_sections diff --git a/docs/api/dtscalibration.plot_residuals_reference_sections_single.rst b/docs/api/dtscalibration.plot_residuals_reference_sections_single.rst new file mode 100644 index 00000000..725d77a0 --- /dev/null +++ b/docs/api/dtscalibration.plot_residuals_reference_sections_single.rst @@ -0,0 +1,6 @@ +plot_residuals_reference_sections_single +======================================== + +.. currentmodule:: dtscalibration + +.. autofunction:: plot_residuals_reference_sections_single diff --git a/docs/api/dtscalibration.plot_sigma_report.rst b/docs/api/dtscalibration.plot_sigma_report.rst new file mode 100644 index 00000000..b047cdeb --- /dev/null +++ b/docs/api/dtscalibration.plot_sigma_report.rst @@ -0,0 +1,6 @@ +plot_sigma_report +================= + +.. currentmodule:: dtscalibration + +.. autofunction:: plot_sigma_report diff --git a/docs/api/dtscalibration.read_apsensing_files.rst b/docs/api/dtscalibration.read_apsensing_files.rst new file mode 100644 index 00000000..04bd67d6 --- /dev/null +++ b/docs/api/dtscalibration.read_apsensing_files.rst @@ -0,0 +1,6 @@ +read_apsensing_files +==================== + +.. currentmodule:: dtscalibration + +.. autofunction:: read_apsensing_files diff --git a/docs/api/dtscalibration.read_sensornet_files.rst b/docs/api/dtscalibration.read_sensornet_files.rst new file mode 100644 index 00000000..541f00cb --- /dev/null +++ b/docs/api/dtscalibration.read_sensornet_files.rst @@ -0,0 +1,6 @@ +read_sensornet_files +==================== + +.. currentmodule:: dtscalibration + +.. autofunction:: read_sensornet_files diff --git a/docs/api/dtscalibration.read_sensortran_files.rst b/docs/api/dtscalibration.read_sensortran_files.rst new file mode 100644 index 00000000..35d92c94 --- /dev/null +++ b/docs/api/dtscalibration.read_sensortran_files.rst @@ -0,0 +1,6 @@ +read_sensortran_files +===================== + +.. currentmodule:: dtscalibration + +.. autofunction:: read_sensortran_files diff --git a/docs/api/dtscalibration.read_silixa_files.rst b/docs/api/dtscalibration.read_silixa_files.rst new file mode 100644 index 00000000..d5adee72 --- /dev/null +++ b/docs/api/dtscalibration.read_silixa_files.rst @@ -0,0 +1,6 @@ +read_silixa_files +================= + +.. currentmodule:: dtscalibration + +.. autofunction:: read_silixa_files diff --git a/docs/api/dtscalibration.shift_double_ended.rst b/docs/api/dtscalibration.shift_double_ended.rst new file mode 100644 index 00000000..0a879ea1 --- /dev/null +++ b/docs/api/dtscalibration.shift_double_ended.rst @@ -0,0 +1,6 @@ +shift_double_ended +================== + +.. currentmodule:: dtscalibration + +.. autofunction:: shift_double_ended diff --git a/docs/api/dtscalibration.suggest_cable_shift_double_ended.rst b/docs/api/dtscalibration.suggest_cable_shift_double_ended.rst new file mode 100644 index 00000000..14a135c0 --- /dev/null +++ b/docs/api/dtscalibration.suggest_cable_shift_double_ended.rst @@ -0,0 +1,6 @@ +suggest_cable_shift_double_ended +================================ + +.. currentmodule:: dtscalibration + +.. autofunction:: suggest_cable_shift_double_ended From cfc4c1b0df46fb596d15c484c3ab00fdc2f361da Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Sun, 3 Sep 2023 16:35:19 +0200 Subject: [PATCH 21/22] Unfortenate formatting of notebooks --- docs/notebooks/03Define_sections.ipynb | 8 +++--- .../04Calculate_variance_Stokes.ipynb | 10 +++---- docs/notebooks/07Calibrate_single_ended.ipynb | 16 ++++++------ docs/notebooks/08Calibrate_double_ended.ipynb | 26 +++++++++---------- .../12Datastore_from_numpy_arrays.ipynb | 16 ++++++------ .../13Fixed_parameter_calibration.ipynb | 6 ++--- docs/notebooks/14Lossy_splices.ipynb | 10 +++---- docs/notebooks/15Matching_sections.ipynb | 8 +++--- docs/notebooks/16Averaging_temperatures.ipynb | 14 +++++----- ...Temperature_uncertainty_single_ended.ipynb | 10 +++---- 10 files changed, 62 insertions(+), 62 deletions(-) diff --git a/docs/notebooks/03Define_sections.ipynb b/docs/notebooks/03Define_sections.ipynb index 8b1000ea..c342cd53 100644 --- a/docs/notebooks/03Define_sections.ipynb +++ b/docs/notebooks/03Define_sections.ipynb @@ -23,7 +23,7 @@ "source": [ "import os\n", "\n", - "from dtscalibration import read_silixa_files\n" + "from dtscalibration import read_silixa_files" ] }, { @@ -40,7 +40,7 @@ "outputs": [], "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"double_ended2\")\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" ] }, { @@ -67,7 +67,7 @@ "source": [ "print(ds.timeseries_keys) # list the available timeseeries\n", "ds.probe1Temperature.plot(figsize=(12, 8))\n", - "# plot one of the timeseries\n" + "# plot one of the timeseries" ] }, { @@ -116,7 +116,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n" + "}" ] }, { diff --git a/docs/notebooks/04Calculate_variance_Stokes.ipynb b/docs/notebooks/04Calculate_variance_Stokes.ipynb index 0de14167..a167e999 100644 --- a/docs/notebooks/04Calculate_variance_Stokes.ipynb +++ b/docs/notebooks/04Calculate_variance_Stokes.ipynb @@ -55,7 +55,7 @@ "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"double_ended2\")\n", "\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" ] }, { @@ -81,7 +81,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n" + "}" ] }, { @@ -125,7 +125,7 @@ " \"is approximately {:.2f} on a {:.1f} sec acquisition time\".format(\n", " I_var, ds.userAcquisitionTimeFW.data[0]\n", " )\n", - ")\n" + ")" ] }, { @@ -152,7 +152,7 @@ " robust=True,\n", " units=\"\",\n", " method=\"single\",\n", - ")\n" + ")" ] }, { @@ -186,7 +186,7 @@ "x = np.linspace(mean - 3 * sigma, mean + 3 * sigma, 100)\n", "approximated_normal_fit = scipy.stats.norm.pdf(x, mean, sigma)\n", "residuals.plot.hist(bins=50, figsize=(12, 8), density=True)\n", - "plt.plot(x, approximated_normal_fit)\n" + "plt.plot(x, approximated_normal_fit)" ] }, { diff --git a/docs/notebooks/07Calibrate_single_ended.ipynb b/docs/notebooks/07Calibrate_single_ended.ipynb index 4c591ede..7edf0a8a 100644 --- a/docs/notebooks/07Calibrate_single_ended.ipynb +++ b/docs/notebooks/07Calibrate_single_ended.ipynb @@ -72,7 +72,7 @@ "outputs": [], "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"single_ended\")\n", - "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")\n" + "ds = read_silixa_files(directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\")" ] }, { @@ -93,7 +93,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(20, 25.5)], # warm bath\n", " \"probe2Temperature\": [slice(5.5, 15.5)], # cold bath\n", - "}\n" + "}" ] }, { @@ -120,7 +120,7 @@ "outputs": [], "source": [ "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", - "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n" + "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")" ] }, { @@ -136,7 +136,7 @@ "metadata": {}, "outputs": [], "source": [ - "resid.plot(figsize=(12, 4))\n" + "resid.plot(figsize=(12, 4))" ] }, { @@ -160,7 +160,7 @@ }, "outputs": [], "source": [ - "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)\n" + "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)" ] }, { @@ -177,7 +177,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds.tmpf.plot(figsize=(12, 4))\n" + "ds.tmpf.plot(figsize=(12, 4))" ] }, { @@ -189,7 +189,7 @@ "ds1 = ds.isel(time=0)\n", "ds1.tmpf.plot(figsize=(12, 4))\n", "(ds1.tmpf_var**0.5).plot(figsize=(12, 4))\n", - "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n" + "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")" ] }, { @@ -217,7 +217,7 @@ "outputs": [], "source": [ "ds1.st.plot(figsize=(12, 8))\n", - "ds1.ast.plot()\n" + "ds1.ast.plot()" ] }, { diff --git a/docs/notebooks/08Calibrate_double_ended.ipynb b/docs/notebooks/08Calibrate_double_ended.ipynb index 21ba4b45..5f36c954 100644 --- a/docs/notebooks/08Calibrate_double_ended.ipynb +++ b/docs/notebooks/08Calibrate_double_ended.ipynb @@ -70,7 +70,7 @@ "\n", "ds_notaligned = read_silixa_files(\n", " directory=filepath, timezone_netcdf=\"UTC\", file_ext=\"*.xml\"\n", - ")\n" + ")" ] }, { @@ -102,7 +102,7 @@ }, "outputs": [], "source": [ - "ds_notaligned = ds_notaligned.sel(x=slice(0, 100)) # only calibrate parts of the fiber\n" + "ds_notaligned = ds_notaligned.sel(x=slice(0, 100)) # only calibrate parts of the fiber" ] }, { @@ -124,7 +124,7 @@ " plot_result=True,\n", " figsize=(12, 6),\n", ")\n", - "print(suggested_shift)\n" + "print(suggested_shift)" ] }, { @@ -140,7 +140,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds = shift_double_ended(ds_notaligned, suggested_shift[0])\n" + "ds = shift_double_ended(ds_notaligned, suggested_shift[0])" ] }, { @@ -160,7 +160,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n" + "}" ] }, { @@ -189,7 +189,7 @@ "st_var, resid = ds.variance_stokes_constant(sections=sections, st_label=\"st\")\n", "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", "rst_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"rast\")\n" + "rast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"rast\")" ] }, { @@ -212,7 +212,7 @@ }, "outputs": [], "source": [ - "resid.plot(figsize=(12, 4))\n" + "resid.plot(figsize=(12, 4))" ] }, { @@ -242,7 +242,7 @@ " ast_var=ast_var,\n", " rst_var=rst_var,\n", " rast_var=rast_var,\n", - ")\n" + ")" ] }, { @@ -258,7 +258,7 @@ }, "outputs": [], "source": [ - "ds.tmpw.plot(figsize=(12, 4))\n" + "ds.tmpw.plot(figsize=(12, 4))" ] }, { @@ -296,7 +296,7 @@ "ds.tmpw_var.plot(figsize=(12, 4))\n", "ds1 = ds.isel(time=-1) # take only the first timestep\n", "(ds1.tmpw_var**0.5).plot(figsize=(12, 4))\n", - "plt.gca().set_ylabel(\"Standard error ($^\\circ$C)\")\n" + "plt.gca().set_ylabel(\"Standard error ($^\\circ$C)\")" ] }, { @@ -333,7 +333,7 @@ " mc_sample_size=500,\n", ") # < increase sample size for better approximation\n", "\n", - "ds.tmpw_mc_var.plot(figsize=(12, 4))\n" + "ds.tmpw_mc_var.plot(figsize=(12, 4))" ] }, { @@ -353,7 +353,7 @@ "ds1.tmpw.plot(linewidth=0.7, figsize=(12, 4))\n", "ds1.tmpw_mc.isel(CI=0).plot(linewidth=0.7, label=\"CI: 2.5%\")\n", "ds1.tmpw_mc.isel(CI=1).plot(linewidth=0.7, label=\"CI: 97.5%\")\n", - "plt.legend(fontsize=\"small\")\n" + "plt.legend(fontsize=\"small\")" ] }, { @@ -382,7 +382,7 @@ "(ds1.tmpb_var**0.5).plot()\n", "(ds1.tmpw_var**0.5).plot()\n", "(ds1.tmpw_mc_var**0.5).plot()\n", - "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n" + "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")" ] }, { diff --git a/docs/notebooks/12Datastore_from_numpy_arrays.ipynb b/docs/notebooks/12Datastore_from_numpy_arrays.ipynb index 1f3d3ce7..9cd233a5 100644 --- a/docs/notebooks/12Datastore_from_numpy_arrays.ipynb +++ b/docs/notebooks/12Datastore_from_numpy_arrays.ipynb @@ -26,7 +26,7 @@ "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "\n", - "from dtscalibration import DataStore, read_silixa_files\n" + "from dtscalibration import DataStore, read_silixa_files" ] }, { @@ -61,7 +61,7 @@ "source": [ "filepath = os.path.join(\"..\", \"..\", \"tests\", \"data\", \"single_ended\")\n", "\n", - "ds_silixa = read_silixa_files(directory=filepath, silent=True)\n" + "ds_silixa = read_silixa_files(directory=filepath, silent=True)" ] }, { @@ -89,7 +89,7 @@ "x = ds_silixa.x.values\n", "time = ds_silixa.time.values\n", "ST = ds_silixa.st.values\n", - "AST = ds_silixa.ast.values\n" + "AST = ds_silixa.ast.values" ] }, { @@ -116,7 +116,7 @@ "ds[\"x\"] = (\"x\", x)\n", "ds[\"time\"] = (\"time\", time)\n", "ds[\"st\"] = ([\"x\", \"time\"], ST)\n", - "ds[\"ast\"] = ([\"x\", \"time\"], AST)\n" + "ds[\"ast\"] = ([\"x\", \"time\"], AST)" ] }, { @@ -133,7 +133,7 @@ "outputs": [], "source": [ "ds = DataStore(ds)\n", - "print(ds)\n" + "print(ds)" ] }, { @@ -169,7 +169,7 @@ "ds[\"temp1\"] = ds_silixa[\"probe1Temperature\"]\n", "ds[\"temp2\"] = ds_silixa[\"probe2Temperature\"]\n", "\n", - "ds.attrs[\"isDoubleEnded\"] = \"0\"\n" + "ds.attrs[\"isDoubleEnded\"] = \"0\"" ] }, { @@ -202,7 +202,7 @@ "ast_var, _ = ds.variance_stokes_constant(sections=sections, st_label=\"ast\")\n", "ds.calibration_single_ended(sections=sections, st_var=st_var, ast_var=ast_var)\n", "\n", - "ds.isel(time=0).tmpf.plot()\n" + "ds.isel(time=0).tmpf.plot()" ] }, { @@ -211,7 +211,7 @@ "metadata": {}, "outputs": [], "source": [ - "ds\n" + "ds" ] }, { diff --git a/docs/notebooks/13Fixed_parameter_calibration.ipynb b/docs/notebooks/13Fixed_parameter_calibration.ipynb index cf5a78c1..fc2642ae 100644 --- a/docs/notebooks/13Fixed_parameter_calibration.ipynb +++ b/docs/notebooks/13Fixed_parameter_calibration.ipynb @@ -76,7 +76,7 @@ " ast_var=ast_var,\n", " fix_gamma=fix_gamma,\n", " fix_dalpha=fix_dalpha,\n", - ")\n" + ")" ] }, { @@ -100,7 +100,7 @@ "outputs": [], "source": [ "print(\"gamma used in calibration:\", ds100.gamma.values)\n", - "print(\"dalpha used in calibration:\", ds100.dalpha.values)\n" + "print(\"dalpha used in calibration:\", ds100.dalpha.values)" ] }, { @@ -132,7 +132,7 @@ " linewidth=1, label=\"Device calibrated\"\n", ") # plot the temperature calibrated by the device\n", "plt.title(\"Temperature at the first time step\")\n", - "plt.legend()\n" + "plt.legend()" ] }, { diff --git a/docs/notebooks/14Lossy_splices.ipynb b/docs/notebooks/14Lossy_splices.ipynb index 5e41b774..977969f3 100644 --- a/docs/notebooks/14Lossy_splices.ipynb +++ b/docs/notebooks/14Lossy_splices.ipynb @@ -72,7 +72,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n" + "}" ] }, { @@ -101,7 +101,7 @@ "ds[\"ast\"] = ds.ast.where(ds.x < 50, ds.ast * 0.82)\n", "\n", "ds[\"rst\"] = ds.rst.where(ds.x > 50, ds.rst * 0.85)\n", - "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)\n" + "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)" ] }, { @@ -121,7 +121,7 @@ "ds.isel(time=0).ast.plot(label=\"ast\")\n", "ds.isel(time=0).rst.plot(label=\"rst\")\n", "ds.isel(time=0).rast.plot(label=\"rast\")\n", - "plt.legend()\n" + "plt.legend()" ] }, { @@ -159,7 +159,7 @@ " rast_var=rast_var,\n", ")\n", "\n", - "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")\n" + "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")" ] }, { @@ -200,7 +200,7 @@ "\n", "ds_a.isel(time=0).tmpw.plot(label=\"no trans. att.\")\n", "ds.isel(time=0).tmpw.plot(label=\"with trans. att.\")\n", - "plt.legend()\n" + "plt.legend()" ] }, { diff --git a/docs/notebooks/15Matching_sections.ipynb b/docs/notebooks/15Matching_sections.ipynb index cb709b50..3017c04b 100644 --- a/docs/notebooks/15Matching_sections.ipynb +++ b/docs/notebooks/15Matching_sections.ipynb @@ -68,7 +68,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0)], # warm bath\n", - "}\n" + "}" ] }, { @@ -95,7 +95,7 @@ "ds[\"ast\"] = ds.ast.where(ds.x < 50, ds.ast * 0.82)\n", "\n", "ds[\"rst\"] = ds.rst.where(ds.x > 50, ds.rst * 0.85)\n", - "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)\n" + "ds[\"rast\"] = ds.rast.where(ds.x > 50, ds.rast * 0.81)" ] }, { @@ -135,7 +135,7 @@ " rast_var=rast_var,\n", ")\n", "\n", - "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")\n" + "ds_a.isel(time=0).tmpw.plot(label=\"calibrated\")" ] }, { @@ -185,7 +185,7 @@ "\n", "ds_a.isel(time=0).tmpw.plot(label=\"normal calibration\")\n", "ds.isel(time=0).tmpw.plot(label=\"matching sections\")\n", - "plt.legend()\n" + "plt.legend()" ] }, { diff --git a/docs/notebooks/16Averaging_temperatures.ipynb b/docs/notebooks/16Averaging_temperatures.ipynb index 76925061..2d97bb85 100644 --- a/docs/notebooks/16Averaging_temperatures.ipynb +++ b/docs/notebooks/16Averaging_temperatures.ipynb @@ -63,7 +63,7 @@ "sections = {\n", " \"probe1Temperature\": [slice(7.5, 17.0), slice(70.0, 80.0)], # cold bath\n", " \"probe2Temperature\": [slice(24.0, 34.0), slice(85.0, 95.0)], # warm bath\n", - "}\n" + "}" ] }, { @@ -82,7 +82,7 @@ "st_var, resid = ds.variance_stokes(sections=sections, st_label=\"st\")\n", "ast_var, _ = ds.variance_stokes(sections=sections, st_label=\"ast\")\n", "rst_var, _ = ds.variance_stokes(sections=sections, st_label=\"rst\")\n", - "rast_var, _ = ds.variance_stokes(sections=sections, st_label=\"rast\")\n" + "rast_var, _ = ds.variance_stokes(sections=sections, st_label=\"rast\")" ] }, { @@ -104,7 +104,7 @@ " ast_var=ast_var,\n", " rst_var=rst_var,\n", " rast_var=rast_var,\n", - ")\n" + ")" ] }, { @@ -191,7 +191,7 @@ " ci_avg_time_isel=[0, 1, 2, 3, 4, 5],\n", " ci_avg_time_sel=None,\n", ")\n", - "ds.tmpw_mc_avg1.plot(hue=\"CI\", linewidth=0.8)\n" + "ds.tmpw_mc_avg1.plot(hue=\"CI\", linewidth=0.8)" ] }, { @@ -252,7 +252,7 @@ " ci_avg_time_isel=[0, 1, 2, 3, 4, 5],\n", " ci_avg_time_sel=None,\n", ")\n", - "ds.tmpw_mc_avg2.plot(hue=\"CI\", linewidth=0.8)\n" + "ds.tmpw_mc_avg2.plot(hue=\"CI\", linewidth=0.8)" ] }, { @@ -313,7 +313,7 @@ " ci_avg_x_sel=slice(7.5, 17.0),\n", " ci_avg_x_isel=None,\n", ")\n", - "ds.tmpw_mc_avgx1.plot(hue=\"CI\", linewidth=0.8)\n" + "ds.tmpw_mc_avgx1.plot(hue=\"CI\", linewidth=0.8)" ] }, { @@ -374,7 +374,7 @@ " ci_avg_x_sel=slice(7.5, 17.0),\n", " ci_avg_x_isel=None,\n", ")\n", - "ds.tmpw_mc_avgx2.plot(hue=\"CI\", linewidth=0.8)\n" + "ds.tmpw_mc_avgx2.plot(hue=\"CI\", linewidth=0.8)" ] }, { diff --git a/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb b/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb index 478a490e..301e4963 100644 --- a/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb +++ b/docs/notebooks/17Temperature_uncertainty_single_ended.ipynb @@ -73,7 +73,7 @@ "plt.fill_between(ds1.x, y1=ds1.tmpf_var, y2=stast_var, label=\"Parameter estimation\")\n", "plt.suptitle(\"Variance of the estimated temperature\")\n", "plt.ylabel(\"$\\sigma^2$ ($^\\circ$C$^2$)\")\n", - "plt.legend(fontsize=\"small\")\n" + "plt.legend(fontsize=\"small\")" ] }, { @@ -84,7 +84,7 @@ "outputs": [], "source": [ "# The effects of the parameter uncertainty can be further inspected\n", - "ds1.var_fw_da.plot(hue=\"comp_fw\", figsize=(12, 4))\n" + "ds1.var_fw_da.plot(hue=\"comp_fw\", figsize=(12, 4))" ] }, { @@ -121,7 +121,7 @@ "source": [ "ds.conf_int_single_ended(\n", " st_var=st_var, ast_var=ast_var, conf_ints=[2.5, 97.5], mc_sample_size=500\n", - ")\n" + ")" ] }, { @@ -145,7 +145,7 @@ "(ds1.tmpf_mc_var**0.5).plot(figsize=(12, 4), label=\"Monte Carlo approx.\")\n", "(ds1.tmpf_var**0.5).plot(label=\"Linear error approx.\")\n", "plt.ylabel(\"$\\sigma$ ($^\\circ$C)\")\n", - "plt.legend(fontsize=\"small\")\n" + "plt.legend(fontsize=\"small\")" ] }, { @@ -166,7 +166,7 @@ "ds1.tmpf.plot(linewidth=0.7, figsize=(12, 4))\n", "ds1.tmpf_mc.sel(CI=2.5).plot(linewidth=0.7, label=\"CI: 2.5%\")\n", "ds1.tmpf_mc.sel(CI=97.5).plot(linewidth=0.7, label=\"CI: 97.5%\")\n", - "plt.legend(fontsize=\"small\")\n" + "plt.legend(fontsize=\"small\")" ] }, { From 044eedf949a9f1f5357d26f49a2ccf732594b55e Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 16 Oct 2023 15:01:25 +0200 Subject: [PATCH 22/22] Fixes to minor comments from Bart --- pyproject.toml | 12 +++++------- src/dtscalibration/averaging_utils.py | 7 +------ src/dtscalibration/datastore.py | 2 +- 3 files changed, 7 insertions(+), 14 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index af0bd1fd..682d0c3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,12 +24,10 @@ readme = "README.rst" license = "BSD-3-Clause" requires-python = ">=3.9, <3.12" authors = [ - {name = "Bas des Tombe", email = "bdestombe@gmail.com"}, - {name = "Bart Schilperoort", email = "b.schilperoort@gmail.com"}, + {name = "Bas des Tombe, Bart Schilperoort"}, ] maintainers = [ - {name = "Bas des Tombe", email = "bdestombe@gmail.com"}, - {name = "Bart Schilperoort", email = "b.schilperoort@gmail.com"}, + {name = "Bas des Tombe, Bart Schilperoort"}, ] keywords = [ "DTS", @@ -56,15 +54,15 @@ dependencies = [ "dask", "pandas", "xarray[parallel]", # numbagg (llvmlite) is a pain to install with pip - "bottleneck", # optional, speed up Xarray - "flox", # optional, speed up Xarray + "bottleneck", # speeds up Xarray + "flox", # speeds up Xarray "pyyaml>=6.0.1", "xmltodict", "scipy", "statsmodels", "matplotlib", "netCDF4>=1.6.4", - "nc-time-axis>=1.4.1" # optional plot dependency of xarray + "nc-time-axis>=1.4.1" # plot dependency of xarray ] dynamic = ["version"] diff --git a/src/dtscalibration/averaging_utils.py b/src/dtscalibration/averaging_utils.py index c94194ba..d17b1f56 100644 --- a/src/dtscalibration/averaging_utils.py +++ b/src/dtscalibration/averaging_utils.py @@ -7,12 +7,7 @@ def inverse_variance_weighted_mean( tmpw_store="tmpw", tmpw_var_store="tmpw_var", ): - """ - Average two temperature datasets with the inverse of the variance as - weights. The two - temperature datasets `tmp1` and `tmp2` with their variances - `tmp1_var` and `tmp2_var`, - respectively. Are averaged and stored in the DataStore. + """Compute inverse variance weighted average, and add result in-place. Parameters ---------- diff --git a/src/dtscalibration/datastore.py b/src/dtscalibration/datastore.py index 459a45dd..d25b4ff6 100644 --- a/src/dtscalibration/datastore.py +++ b/src/dtscalibration/datastore.py @@ -218,7 +218,7 @@ def sections(self, value): "Not possible anymore. Instead, pass the sections as an argument to \n" "ds.dts.calibrate_single_ended() or ds.dts.calibrate_double_ended()." ) - raise NotImplementedError(msg) + raise DeprecationWarning(msg) def check_reference_section_values(self): """