Merge pull request #139 from NREL/postfix

Fix compute fitted kla

Author: malihass

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)
 

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,

bird/version.py

@@ -1,3 +1,3 @@
 """Bio reactor design version"""
 
-__version__ = "0.0.43"
+__version__ = "0.0.45"

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.
 
 
 

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()