From d5b287f97519fdab67dc475de2347ecdaacbc25d Mon Sep 17 00:00:00 2001 From: Malik Date: Fri, 3 Oct 2025 09:17:31 -0600 Subject: [PATCH 1/3] remove rogue breakpoint --- bird/postprocess/kla_utils.py | 4 ++-- bird/postprocess/post_quantities.py | 16 +++++++++++++--- bird/version.py | 2 +- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/bird/postprocess/kla_utils.py b/bird/postprocess/kla_utils.py index 7720ee27..b08bfd65 100644 --- a/bird/postprocess/kla_utils.py +++ b/bird/postprocess/kla_utils.py @@ -147,14 +147,14 @@ def check_data_size(time_obs: np.ndarray, concentration_obs: np.ndarray): def compute_kla( data_t: np.ndarray, data_c: np.ndarray, + num_warmup: int = 4000, + num_samples: int = 1000, max_chop: int | None = None, bootstrap: bool = True, ) -> dict: """ Main loop to compute kla """ - num_warmup = 4000 - num_samples = 1000 rng_key = random.PRNGKey(0) rng_key, rng_key_ = random.split(rng_key) diff --git a/bird/postprocess/post_quantities.py b/bird/postprocess/post_quantities.py index a7fbf0a6..13456859 100644 --- a/bird/postprocess/post_quantities.py +++ b/bird/postprocess/post_quantities.py @@ -1117,6 +1117,8 @@ def compute_fitted_kla( species_names: str | list[str], n_cells: int | None = None, volume_time: str | None = None, + num_warmup: int = 4000, + num_samples: int = 1000, field_dict: dict | None = None, ) -> tuple[dict, dict, dict]: r""" @@ -1149,6 +1151,12 @@ def compute_fitted_kla( volume_time : str | None Time folder to read to get the cell volumes. If None, finds volume time automatically + num_warmup: int + Number of MCMC samples in the warmup phase + Defaults to 4000 + num_samples: int + Number of posterior MCMC samples generated + Defaults to 1000 field_dict : dict Dictionary of fields used to avoid rereading the same fields to calculate different quantities @@ -1228,16 +1236,18 @@ def compute_fitted_kla( ) c_history[species_name][itime] = c_liq - breakpoint() - logger.info("Doing kla fit") # Compute kla kla_spec = {} cstar_spec = {} for species_name in species_names: kla_res = compute_kla( - np.array(time_float_sorted), c_history[species_name] + np.array(time_float_sorted), + c_history[species_name], + num_warmup=num_warmup, + num_samples=num_samples, ) + # Convert to h-1 kla_spec[species_name] = { "mean": kla_res["kla"] * 3600, "std": kla_res["kla_err"] * 3600, diff --git a/bird/version.py b/bird/version.py index 9154fcd5..68c45d06 100644 --- a/bird/version.py +++ b/bird/version.py @@ -1,3 +1,3 @@ """Bio reactor design version""" -__version__ = "0.0.43" +__version__ = "0.0.45" From 9ba625427024c47fbab994c1ba45ec831d538ca3 Mon Sep 17 00:00:00 2001 From: Malik Date: Fri, 3 Oct 2025 09:34:19 -0600 Subject: [PATCH 2/3] add test for kla compute --- tests/postprocess/test_post_quantities.py | 58 ++++++++++++++++++++++- 1 file changed, 56 insertions(+), 2 deletions(-) diff --git a/tests/postprocess/test_post_quantities.py b/tests/postprocess/test_post_quantities.py index 05d0f6c1..9ecde554 100644 --- a/tests/postprocess/test_post_quantities.py +++ b/tests/postprocess/test_post_quantities.py @@ -315,8 +315,62 @@ def test_instantaneous_kla(): ) +def test_fitted_kla(): + """ + Test for fitted kla calculation + """ + case_folder = os.path.join( + Path(__file__).parent, + "..", + "..", + "bird", + "postprocess", + "data_conditional_mean", + ) + field_dict = {} + # Check that assertion error is sent if too few snapshots + try: + kla_spec1, cstar_spec1, _ = compute_fitted_kla( + species_names=["CO2"], + case_folder=case_folder, + ) + except AssertionError: + pass + + # Create dummy time folders + for time_folder in [str(entry) for entry in range(81, 89)]: + shutil.copytree( + os.path.join(case_folder, "80"), + os.path.join(case_folder, time_folder), + ) + kla_spec1, cstar_spec1, _ = compute_fitted_kla( + species_names=["CO2"], + case_folder=case_folder, + num_warmup=100, + num_samples=100, + ) + for time_folder in [str(entry) for entry in range(81, 89)]: + shutil.rmtree(os.path.join(case_folder, time_folder)) + + # Create dummy time folders to trigger bootstrapping + for time_folder in [str(entry) for entry in range(81, 91)]: + shutil.copytree( + os.path.join(case_folder, "80"), + os.path.join(case_folder, time_folder), + ) + kla_spec1, cstar_spec1, _ = compute_fitted_kla( + species_names=["CO2"], + case_folder=case_folder, + num_warmup=100, + num_samples=100, + ) + for time_folder in [str(entry) for entry in range(81, 91)]: + shutil.rmtree(os.path.join(case_folder, time_folder)) + + if __name__ == "__main__": - test_compute_superficial_gas_velocity() - test_compute_gh() + # test_compute_superficial_gas_velocity() + # test_compute_gh() # test_ave_y_liq() # test_ave_conc_liq() + test_fitted_kla() From 77558502f486e5b6ec4b81046da17ffe2b8cd09a Mon Sep 17 00:00:00 2001 From: Malik Date: Fri, 3 Oct 2025 09:36:40 -0600 Subject: [PATCH 3/3] explicit doc --- docs/source/python_interface_tut.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/python_interface_tut.rst b/docs/source/python_interface_tut.rst index a1547d41..5a6a0c0a 100644 --- a/docs/source/python_interface_tut.rst +++ b/docs/source/python_interface_tut.rst @@ -50,7 +50,7 @@ If we want to read the gas volume fraction ``alpha.gas``, the mass fraction of C print("cell alpha gas shape = ", alpha_gas.shape) print("cell u liq shape = ", u_liq.shape) -The function ``read_field`` automatically detects whether the field is a vector or a scalar field. +The function ``read_field`` automatically detects whether the field is a vector or a scalar field. The first returned results in (N,) for a scalar field and (N,3) for a vector field.