Add rudimentary doc

.github/workflows/ci.yml

@@ -163,3 +163,15 @@ jobs:
             exit 2
           fi 
           cd ../../../
+      - name: NPE training
+        run: |
+          cd scripts/reg_tests/4.npe
+          bash run_change_recipe_github.sh
+          bash run.sh
+          if [ ! -f "models/post_test.txt" ]; then
+            exit 1
+          fi 
+          if [ ! -f "models/post_val.txt" ]; then
+            exit 2
+          fi 
+          cd ../../../

README.md

@@ -2,12 +2,64 @@
 
 ## Summary
 
+This package implements several data-based techniques for parameter fitting in Li-ion battery models. It uses [BATMODS-lite](https://github.com/NatLabRockies/batmods-lite) to generate the data. 
+
+The repository contains the code that is used for the paper "Neural posterior estimation is accurate, tractable and scalable for inverse parameter inference in Li-ion batteries", M. Hassanaly, C. R. Randall, P. J.  Weddle, P. J. Gasper, C. Kelly, T. Tanim, E. J. Dufek, K. Smith.
+
 ## Installation
 
+We recommend using a conda environment
+
+```
+conda create -n batfit python=3.12
+conda activate batfit
+```
 
+Once the files are available on your machine, use your terminal to navigate into the folder and execute one of the following depending on your installation preference.
 
 ```
 pip install .             (basic installation)
 pip install -e .[dev]     (editable installation with developer options)
 ```
 
+## Get Started
+
+The regression tests (run as part of the CI) show how to use the basic capabilities of the code.
+1. Data generation with BATMODS-lite (``scripts/reg_tests/1.gen_data``). This demonstrates how to generate data from a single particle model in parallel. This has been tested with up to 1200 workers but uses 4 workers here. 
+Once BATFIT is installed
+```
+bash run.sh
+``` 
+
+2. Train a surrogate of the physics-based model (``scripts/reg_tests/2.train_surrogate``). This demonstrate how to preprocess the data and train a surrogate model of the single particle model that can be used for simulation-based inference.
+Once BATFIT is installed, and ``run_change_recipe_local.sh`` has been changed to used your path.
+```
+bash run_change_recipe_local.sh 
+bash run.sh
+```
+
+3. Run MCMC to identify parameters with the trained surrogate (``scripts/reg_tests/3.surrogate_mcmc``). This demonstrates how to run MCMC with a data-based surrogate instead of a physics-based model.
+Once BATFIT is installed, and ``run_change_recipe_local.sh`` has been changed to used your path.
+```
+bash run_change_recipe_local.sh 
+bash run.sh
+```
+
+4. Use Neural Posterior Estimation to approximate the parameter posterior PDF (``scripts/reg_tests/4.npe``)
+Once BATFIT is installed, and ``run_change_recipe_local.sh`` has been changed to used your path.
+```
+bash run_change_recipe_local.sh
+bash run.sh
+```
+
+## Citing this Work
+
+```
+SWR BatFit
+
+arxiv citation
+```
+
+## Acknowledgements
+This work was authored by the National Laboratory of the Rockies (NLR) for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This work was supported by funding from DOE's Vehicle Technologies Office (VTO). The research was performed using computational resources sponsored by the Department of Energy's Office of Critical Minerals and Energy Innovation (CMEI) and located at the National Laboratory of the Rockies. The views expressed in the repository do not necessarily represent the views of the DOE or the U.S. Government.
+

batfit/calibration/cal_utils.py

@@ -1207,7 +1207,8 @@ def mcmc_iter(
         )
 
     if cal_sigma and read_sigma:
-        sys.exit("ERROR: cal sigma OR read sigma")
+        logger.error("Cal sigma OR read sigma")
+        sys.exit()
     if cal_sigma:
         theta.append(np.random.uniform(min_sigma, max_sigma))
     if read_sigma:

batfit/calibration/data_utils.py

@@ -183,7 +183,8 @@ def read_cycling_data(
         charging = True
         cccv = True
     else:
-        sys.exit(f"ERROR: mode {cyc_mode} not recognized")
+        logger.error(f"Mode {cyc_mode} not recognized")
+        sys.exit()
     if discharging:
         print("Discharging")
     if charging and cc:
@@ -292,14 +293,16 @@ def read_cycling_data(
         ind_max = np.argmax(c_arr)
         max_c_tmp = c_arr[ind_max]
         if min_c_tmp < 0.4:
-            print(f"ERROR: {cellFolder} for cycle {cyc}")
-            print(f"min c {min_c_tmp} index {ind_min}")
-            print(c_arr)
+            msg = f"ERROR: {cellFolder} for cycle {cyc}"
+            msg += f"\nmin c {min_c_tmp} index {ind_min}"
+            msg += f"\n{c_arr}"
+            logger.error(msg)
             sys.exit()
         if max_c_tmp > 15:
-            print(f"ERROR: {cellFolder} for cycle {cyc}")
-            print(f"max c {max_c_tmp} index {ind_max}")
-            print(c_arr)
+            msg = f"ERROR: {cellFolder} for cycle {cyc}"
+            msg += f"\nmin c {max_c_tmp} index {ind_max}"
+            msg += f"\n{c_arr}"
+            logger.error(msg)
             sys.exit()
 
         if max_c is None or max_c_tmp > max_c:
@@ -435,7 +438,8 @@ def read_cycling_data(
                         print("neighbors = ", neighbors)
                         print("cycle_extract = ", cycle_extract)
                         print("cycles keys ", list(cycles.keys()))
-                        sys.exit("ERROR: missing cycles")
+                        logger.error("Missing cycles")
+                        sys.exit()
                 else:
                     if target == "phis_c":
                         data_ext_x = np.hstack((data_ext_x, t_interp))
@@ -1604,7 +1608,8 @@ def mcmc_iter(
         )
 
     if cal_sigma and read_sigma:
-        sys.exit("ERROR: cal sigma OR read sigma")
+        logger.error("Cal sigma OR read sigma")
+        sys.exit()
     if cal_sigma:
         theta.append(np.random.uniform(min_sigma, max_sigma))
     if read_sigma:

batfit/model/paramNN.py

@@ -1902,9 +1902,10 @@ def compute_post(
                 if post_fn in [accuracy, rel_accuracy]:
                     post_val = post_fn(pred, batch[1].to(device))
                 elif post_fn in [identifiability]:
-                    sys.exit(
-                        "ERROR: identifiability can only be computed with probabilistic model"
+                    logger.error(
+                        "Identifiability can only be computed with probabilistic model"
                     )
+                    sys.exit()
             elif isinstance(model, ProbParamCNN) or isinstance(
                 model, ProbParamFCNN
             ):

batfit/preprocess/pickledb.py

@@ -3,6 +3,8 @@
 import sys
 import time
 
+from batfit import logger
+
 
 class PickleDB:
     def __init__(self, filename, read_from_existing=False):
@@ -45,9 +47,11 @@ def read(self, max_try=10):
                 try_count += 1
                 time.sleep(random.random() + 1)
             except FileNotFoundError:
-                sys.exit("ERROR: database must be created before it is read")
+                logger.error("Database must be created before it is read")
+                sys.exit()
         if try_count >= max_try:
-            sys.exit(f"ERROR: could not read DB after {max_try} tries")
+            logger.error(f"Could not read DB after {max_try} tries")
+            sys.exit()
 
         return records
 
@@ -67,8 +71,10 @@ def append(self, data, max_try=10):
                 try_count += 1
                 time.sleep(random.random() + 1)
             except FileNotFoundError:
-                sys.exit(
-                    "ERROR: database must be created before it is written to"
+                logger.error(
+                    "Database must be created before it is written to"
                 )
+                sys.exit()
         if try_count >= max_try:
-            sys.exit(f"ERROR: could not write to DB after {max_try} tries")
+            logger.error(f"Could not write to DB after {max_try} tries")
+            sys.exit()

batfit/preprocess/sim_setup.py

@@ -56,7 +56,7 @@ def parse_input(filename, parallel_env=None):
             for param_name in deg_param_names
         ]
     except KeyError:
-        print("ERROR: mismatch of parameters")
+        logger.error("Mismatch of parameters")
         raise KeyError
 
     if cyc_mode == "discharge-chargecc":
@@ -70,7 +70,7 @@ def parse_input(filename, parallel_env=None):
                 for param_name in deg_param_names_chcc
             ]
         except KeyError:
-            print("ERROR: mismatch of parameters")
+            logger.error("Mismatch of parameters")
             raise KeyError
         deg_param_names += deg_param_names_chcc_new
         deg_param_min += deg_param_min_chcc

batfit/preprocess/sol_gen.py

@@ -286,7 +286,8 @@ def robust_LHRH(sim, df, charge, protocol, sim_params, bat_model):
         elif bat_model.lower() == "p2d":
             sim_prot = bm.P2D.CycleSolution(*all_solns)
         else:
-            sys.exit("ERROR: battery model not recognized")
+            logger.error("Battery model not recognized")
+            sys.exit()
         rootsol = sim_prot
         assert all(rootsol.success)
     except:
@@ -335,7 +336,8 @@ def robust_LHRH(sim, df, charge, protocol, sim_params, bat_model):
                 elif bat_model.lower() == "p2d":
                     sim_prot = bm.P2D.CycleSolution(*all_solns)
                 else:
-                    sys.exit("ERROR: battery model not recognized")
+                    logger.error("Battery model not recognized")
+                    sys.exit()
                 rootsol = sim_prot
                 assert all(rootsol.success)
                 break
@@ -634,7 +636,8 @@ def single_run_save(
         run_p2d = False
         run_spm = True
     elif rootsol is not None:
-        sys.exit("ERROR: battery model not recognized")
+        logger.error("Battery model not recognized")
+        sys.exit()
     try:
         assert rootsol is not None
         if only_phi_CC:
@@ -713,7 +716,7 @@ def single_run_save(
 
         for key in diff_dict:
             if key in save_dict:
-                raise ValueError(f"ERROR: save_dict already contains {key}")
+                raise ValueError(f"Save_dict already contains {key}")
             save_dict[key] = diff_dict[key]
 
         # Reduce if needed

batfit/utils/data_utils.py

@@ -44,7 +44,7 @@ def get_max_time(data_root_folder):
 
 
 def from_sol_to_data(data_root_folder, filename, n_points):
-    print("ERROR: Use pkl file option instead")
+    logger.error("Use pkl file option instead")
     raise NotImplementedError
     sol = np.load(os.path.join(data_root_folder, filename))
     min_t = np.amin(sol["t"])

scripts/reg_tests/2.train_surrogate/run_change_recipe_local.sh

@@ -0,0 +1,6 @@
+in_yml=training_recipes/recipe_rel.yml
+out_yml=training_recipes/recipe.yml
+r_path=/Users/mhassana/Desktop/GitHub/BatFit/scripts/reg_tests
+s_path=/Users/mhassana/Desktop/GitHub/BatFit/batfit/default_exps
+
+python write_absolute_recipe.py $in_yml $out_yml $r_path $s_path

scripts/reg_tests/3.surrogate_mcmc/run_change_recipe_local.sh

@@ -0,0 +1,5 @@
+in_yml=cal_recipes/recipe_rel.yml
+out_yml=cal_recipes/recipe.yml
+r_path=/Users/mhassana/Desktop/GitHub/BatFit/scripts/reg_tests
+
+python write_absolute_recipe.py $in_yml $out_yml $r_path

scripts/reg_tests/3.surrogate_mcmc/test_mcmc.py

@@ -121,7 +121,7 @@ def forward(self, degradation_parameters: list):
         return output[:, 0]
 
 
-def load_synthetic_data(inp, noisy=True):
+def load_synthetic_data(inp):
     t = {}
     phi = {}
     truth = {}
@@ -136,23 +136,16 @@ def load_synthetic_data(inp, noisy=True):
         cyc_mode=inp.cyc_mode,
     )
     data_path = inp.data_path_discharge
-    tmp = np.load(os.path.join(data_path, "data_split.npz"))
-    with open(os.path.join(data_path, "scaler_X.pkl"), "rb") as f:
-        scaler = pickle.load(f)
-    if noisy:
-        batch_in_scaled = apply_noise(
-            torch.tensor(scaler.transform(tmp["X_test"])),
-            scaler_X=scaler,
-            noise_levels=noise_levels,
-            a_min=a_min,
-            a_max=a_max,
-        )
-    else:
-        batch_in_scaled = torch.tensor(scaler.transform(tmp["X_test"]))
-    batch_in_unscaled = scaler.inverse_transform(batch_in_scaled).numpy()
+    tmp = np.load(os.path.join(data_path, "assembled_data.npz"))
+    batch_in_unscaled = apply_noise_unscaled(
+        torch.tensor(tmp["X_data"]),
+        noise_levels=noise_levels,
+        a_min=a_min,
+        a_max=a_max,
+    )
     t["discharge"] = batch_in_unscaled[:, 0, :]
     phi["discharge"] = batch_in_unscaled[:, 1, :]
-    truth["discharge"] = tmp["Y_test"][:, :]
+    truth["discharge"] = tmp["Y_data"][:, :]
 
     return (
         t,
@@ -214,7 +207,7 @@ def compute_fit_error(inp, samples, models):
 
         for key in models:
             # fig=plt.figure()
-            true_v = data_phis_c[key]
+            true_v = np.array(data_phis_c[key])
             # plt.plot(true_v, color="k")
             for j_sample_test in range(samples.shape[1]):
                 pred_v = np.array(
@@ -260,7 +253,7 @@ def compute_fit_error(inp, samples, models):
         total_data_t,
         total_data_phi,
         total_truth,
-    ) = load_synthetic_data(inp, noisy=False)
+    ) = load_synthetic_data(inp)
 
     cycle_types = list(models.keys())
 

scripts/reg_tests/4.npe/run.sh

@@ -0,0 +1,4 @@
+# Train surrogate
+python train_nn.py training_recipes/recipe.yml
+# Test surrogate
+python test_nn.py training_recipes/recipe.yml

scripts/reg_tests/4.npe/run_change_recipe_github.sh

@@ -0,0 +1,6 @@
+in_yml=training_recipes/recipe_rel.yml
+out_yml=training_recipes/recipe.yml
+r_path=$GITHUB_WORKSPACE/scripts/reg_tests
+s_path=$GITHUB_WORKSPACE/batfit/default_exps
+
+python write_absolute_recipe.py $in_yml $out_yml $r_path $s_path

scripts/reg_tests/4.npe/run_change_recipe_local.sh

@@ -0,0 +1,6 @@
+in_yml=training_recipes/recipe_rel.yml
+out_yml=training_recipes/recipe.yml
+r_path=/Users/mhassana/Desktop/GitHub/BatFit/scripts/reg_tests
+s_path=/Users/mhassana/Desktop/GitHub/BatFit/batfit/default_exps
+
+python write_absolute_recipe.py $in_yml $out_yml $r_path $s_path