Demonstrator calibration

environment.yml

@@ -40,6 +40,7 @@ dependencies:
   - codecov
   - pydantic<2.0
   - seaborn
+  - pytorch==2.0.0
   - pip:
     - probeye==2.3.0
     - -e .

lebedigital/demonstrator_optimization_scripts/design_variable_to_kpi.py

@@ -0,0 +1,69 @@
+import json
+import os, sys
+from lebedigital.demonstrator_optimization_scripts.utils import load_json, update_json
+
+def design_var_to_kpi(workflow_path:str,X: dict, seed: int) -> dict:
+    """
+    Runs the snakemake workflow and the returns the KPIs for objective and constraints for a given value of the design
+    variables. The Random variables (b) x->b->KPIs are also sampled for a given value of seed.
+    Args:
+     X: (dict) with keys 'agg_ratio' (volume ratio of the aggregates) and 'slag_ratio'
+     seed: the seed parameter. This ensures that the sampled Random variable here is the same as the one passed in the
+     forward call
+
+    Returns:
+        y : dict with all the KPIs
+
+    """
+    # Pass the parameter to X to the input to forward. Meaning overwrrite the input.
+    # The design variables, aggregate ratio and the slag ratio needs to be updated.
+
+    #TODO: the below is hardcoded. fixit
+    #design_var_paths = {'aggregates_volume_fraction': workflow_path +'/Inputs/aggregates_volume_fraction.json',
+    #                    'sc_volume_fraction': workflow_path + '/Inputs/sc_fraction.json'}
+    design_var_paths = {'height': workflow_path + '/Inputs/geometry.json',
+                        'sc_volume_fraction': workflow_path + '/Inputs/sc_fraction.json'}
+
+    for key, value in X.items():
+        update_json(design_var_paths[key],key,value)
+
+    # pass the seed to the scripts for the RVs (see eqn 29 SVO paper)
+    # Updating the phi's which are input to the script.
+    phi_hydration_path = workflow_path + '/Inputs/phi_hydration.json'
+    phi_paste_path = workflow_path + '/Inputs/phi_paste.json'
+    update_json(phi_hydration_path, 'seed', seed)
+    update_json(phi_paste_path, 'seed', seed)
+
+    # Run the workflow using snakemake
+    # add the path to the workflow file and the path to the directory
+    workflow_file_path = workflow_path + '/Snakefile'
+    os.system(f'snakemake --cores 7 --snakefile {workflow_file_path} '
+              f'--directory {workflow_path}  workflow_targets --use-conda')
+
+    # Read in the KPIs in a dict
+    Results_path = workflow_path + '/Results/'
+    FEM_KPI = Results_path + 'kpi_from_fem.json'
+    gwp_KPI = Results_path + 'gwp_beam.json'
+    beam_design_KPI = Results_path + 'beam_design.json'
+    y = {}
+    for i, path in enumerate([FEM_KPI, gwp_KPI, beam_design_KPI]):
+        tmp = load_json(path)
+        y.update(tmp)
+
+    # return the KPIs
+    #TODO: this is specific to the constraints and objective choosen. careful
+    kpi = {
+        "gwp_beam": y["gwp_beam"]["value"],
+       # "check_steel_area": y["check_steel_area"]["value"],
+        "constraint_beam_design": y["constraint_beam_design"]["value"],
+        "max_reached_temperature": y["max_reached_temperature"]["value"],
+        "time_of_demoulding": y["time_of_demolding"]["value"]
+    }
+    return kpi
+
+if __name__ == '__main__':
+    path = '../../usecases/optimization_paper/1'
+    design_var = {'height': 260,
+                        'sc_volume_fraction': 0.35}
+    seed = 66
+    design_var_to_kpi(workflow_path=path,X=design_var,seed=seed)

lebedigital/demonstrator_optimization_scripts/farm_workflow.py

@@ -0,0 +1,33 @@
+import os
+import shutil
+import numpy as np
+import sys
+
+def farm_workflow(path:str,seed:list):
+    """
+    Created multiple copies of snakemake workflow folder to be used later to parallelize
+    Parameters
+    ----------
+    path: should point to the optimization_paper path
+    seed
+
+    Returns
+    -------
+
+    """
+
+    for i,v in enumerate(seed):
+        new_dir_path = os.path.join(path,str(i+1))
+        src_path = path + '/optimization_workflow'
+        # copy to the newly created folder
+        if not os.path.exists(new_dir_path):
+            shutil.copytree(src=src_path,dst=new_dir_path)
+        else:
+            shutil.rmtree(new_dir_path)
+            shutil.copytree(src=src_path,dst=new_dir_path)
+
+if __name__ == '__main__':
+    path = '../../usecases/optimization_paper'
+    #seed = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
+    seed = np.load('../../usecases/demonstrator/Calibration/seed_tmp.npy').astype(int).tolist()
+    farm_workflow(path,seed)

lebedigital/demonstrator_optimization_scripts/parallel_compute_workflow.py

@@ -0,0 +1,45 @@
+import json
+import sys
+import numpy as np
+from lebedigital.demonstrator_optimization_scripts.design_variable_to_kpi import design_var_to_kpi
+
+def parallel_workflow(array_id,design_var:np.ndarray, seed:list):
+    """
+    For the given array_id of job, points the created folder, runs the workflow and
+    saves a kpi dictionary.
+    Parameters
+    ----------
+    array_id
+    design_var: rows are the samples and columns the design variable number.
+    seed
+
+    Returns
+    -------
+
+    """
+    path = '../../usecases/optimization_paper/' + str(array_id)
+    # TODO: pass the below aslo, not hardcode
+
+    idx = array_id -1
+    # design_var = {'aggregates_volume_fraction':design_var[idx,0], #0.4
+    #               'sc_volume_fraction': design_var[idx,1]} #0.35
+    design_var = {'height': design_var[idx, 0],  # 0.4
+                  'sc_volume_fraction': design_var[idx, 1]}  # 0.35
+    kpi = design_var_to_kpi(workflow_path=path, X=design_var, seed=seed[idx])
+    kpi_path = path + '/kpi.json'
+    with open(kpi_path, 'w') as f:
+        json.dump(kpi, f, indent=4, sort_keys=True)
+
+# to pass the job array number here.
+#TODO: read from the seed file
+#seeds = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]
+
+#seed = np.random.randint(666,size=5)
+#np.save('../../usecases/demonstrator/Calibration/seed_tmp.npy',seed)
+#des_var = np.random.uniform(size=(5,2))
+#np.save('../../usecases/demonstrator/Calibration/design_var_tmp.npy',des_var)
+
+seeds = np.load('../../usecases/demonstrator/Calibration/seed_tmp.npy').astype(int).tolist()
+design_var = np.load('../../usecases/demonstrator/Calibration/design_var_tmp.npy')
+parallel_workflow(int(sys.argv[1]),design_var=design_var,seed=seeds)
+#parallel_workflow(4,design_var=design_var,seed=seeds)

lebedigital/demonstrator_optimization_scripts/run_jobs.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+#SBATCH --job-name=LBD_optimization
+#SBATCH --nodes=1
+#SBATCH --ntasks=1
+#SBATCH --cpus-per-task=3
+#SBATCH --partition=batch_SKL,batch_SNB
+#SBATCH --array=1-100
+#SBATCH --output=slurm-%A_%a.out
+
+# 1-$1, the $1 is for the first argument of the sbatch run_jobs 100, where 100 samples total
+# Load Python Module
+source /home/atul/.bashrc
+
+. /home/atul/miniconda3/etc/profile.d/conda.sh
+
+conda activate lebedigital
+
+# Print to a file a message that includes the current $SLURM_ARRAY_TASK_ID, the same name, and the sex of the sample
+#echo "This is array task ${SLURM_ARRAY_TASK_ID}" >> output.txt
+
+# creates multiple folders for each sample
+#python farm_workflow.py
+
+# Run the workflow in folder parallely
+python parallel_compute_workflow.py $SLURM_ARRAY_TASK_ID
+
+echo "!!!! Job array completed.!!!!"

lebedigital/demonstrator_optimization_scripts/utils.py

@@ -0,0 +1,95 @@
+import os, sys
+import json
+
+import numpy as np
+
+from lebedigital.demonstrator_optimization_scripts.farm_workflow import farm_workflow
+
+def python_fn_run_jobs(path_to_scripts_folder:str,no_samples:int):
+    """
+    Create copies of the workflow for each sample, and run slurm job array.
+    Parameters
+    ----------
+    path_to_scripts_folder
+    no_samples
+
+    Returns
+    -------
+
+    """
+    #current_dir = os.path.dirname(os.path.abspath(__file__))
+    #path_farm = path_to_scripts_folder + "/farm_workflow.py"
+    #path_jobs = path_to_scripts_folder + 'run_jobs.sh'
+
+    print("!!! Creating folders for copies of the workflow !!!")
+    #os.system(f'python {path_farm}')
+    seed = np.load(path_to_scripts_folder+'../../usecases/demonstrator/Calibration/seed_tmp.npy').astype(int).tolist()
+    farm_workflow(path=path_to_scripts_folder+'../../usecases/optimization_paper',
+                  seed=seed)
+
+    print("!!! folder creating DONE !!!")
+    print("!!! Run multiple jobs in cluster !!!")
+    # change directory to the file in which this fn is.
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    original_dir = os.getcwd()
+    os.chdir(script_dir)
+    os.system(f'sbatch --wait --array=1-{no_samples} run_jobs.sh')
+    if not os.path.exists(f'../../usecases/optimization_paper/{no_samples}/kpi.json'):
+        raise FileNotFoundError
+    print('All jobs finished')
+    # restore to the working directory
+    os.chdir(original_dir)
+
+# %%
+def load_json(path: str) -> dict:
+    if path[-5:] == '.json':
+        with open(path) as f:
+            data = json.load(f)
+    return data
+
+
+# %%
+def update_json(file_path: str, key: str, value):
+    # Read the JSON file
+    with open(file_path, 'r') as f:
+        data = json.load(f)
+    # TODO:will work only when 'value' key is present
+    # Update the value of the specified key
+    data[key]['value'] = value
+
+    # Write the updated data back to the JSON file
+    with open(file_path, 'w') as f:
+        json.dump(data, f, indent=4, sort_keys=True)
+
+def read_kpis(kpi_path:str):
+    """
+    Read in the kpis from the
+    Parameters
+    ----------
+    kpi_path
+
+    Returns
+    -------
+
+    """
+
+    if not os.path.exists(kpi_path):
+        print(f"Error: File {kpi_path} does not exist.")
+    data = load_json(kpi_path)
+    #TODO: the below is specific to the problem
+    # print("!!! Attention the KPIs are specific and can change. Careful.")
+    obj = data["gwp_beam"]
+    C_1 = data["constraint_beam_design"]
+    C_2 = data["max_reached_temperature"]
+    C_3 = data["time_of_demoulding"]
+
+    return obj, C_1, C_2, C_3
+
+if __name__=='__main__':
+    #o,c1,c2,c3 =read_kpis(kpi_path='../../usecases/optimization_paper/1/kpi.json')
+    #print(o,c1,c2,c3)
+
+    # test function
+    python_fn_run_jobs('./',5)
+
+

lebedigital/demonstrator_scripts/dummy_hydration_parameters.py

@@ -1,16 +1,22 @@
 import numpy as np
-
 from lebedigital.unit_registry import ureg
+import torch as th
+import numpy as np
+
+th.set_default_dtype(th.float64)
 
 
-@ureg.check("", "")
-def dummy_hydration_parameters(slag_ratio, phi_hydration):
+@ureg.check('', '', '', '')
+def dummy_hydration_parameters(slag_ratio, phi_mean: list, phi_cov: list, seed: int):
     """
     This is a dummy function to make the snakemake workflow work until the real function is ready
     It changes arbitrarily chosen values depending on slag content, not based on physics or anything.
 
     Parameters
     ----------
+    seed : int
+    phi_cov : list
+    phi_mean : list
     slag_ratio : float / pint unitless
         amount of slag compared to cement, value from 0 to 1
     phi_hydration: ??
@@ -32,15 +38,49 @@ def dummy_hydration_parameters(slag_ratio, phi_hydration):
         maximum potential hydration parameter
     """
 
-    B1 = 2.916e-4 * ureg("1/s")  # in 1/s
-    B2 = 0.0024229 * ureg("")  # -
-    eta = 5.554 * ureg("")  # something about diffusion
-    E_act = 5653 * 8.3145 * ureg("J/mol")  # activation energy in Jmol^-1
+    # paste_youngs_modulus_min = 15 * ureg('GPa')
+    # paste_youngs_modulus_max = 25 * ureg('GPa')
+    # paste_youngs_modulus = paste_youngs_modulus_min + (paste_youngs_modulus_max-paste_youngs_modulus_min)*slag_ratio
+    #
+    # paste_strength_min = 8 * ureg('MPa')
+    # paste_strength_max = 13 * ureg('MPa')
+    # paste_strength = paste_strength_min + (paste_strength_max-paste_strength_min)*slag_ratio
+
+    B1 = 2.916E-4 * ureg('1/s')  # in 1/s
+    B2 = 0.0024229 * ureg('')  # -
+    eta = 5.554 * ureg('')  # something about diffusion
+    E_act = 3653 * 8.3145 * ureg('J/mol')  # activation energy in Jmol^-1
     Q_ = ureg.Quantity
     T_ref = Q_(25, ureg.degC)
 
-    Q_pot_min = 100000 * ureg("J/kg")
-    Q_pot_max = 300000 * ureg("J/kg")
+    Q_pot_min = 200000 * ureg('J/kg')
+    Q_pot_max = 300000 * ureg('J/kg')
     Q_pot = Q_pot_max - (Q_pot_max - Q_pot_min) * slag_ratio
 
-    return B1, B2, eta, E_act, Q_pot, T_ref
+    # ATUL : temporary q(b|x)~N(mu,cov), b=(B1,B2,eta,E_ect,Q_pot,T_ref), x = slag ratio
+    th.manual_seed(seed=seed)
+    no_of_parameter = 6
+    assert len(phi_mean) == no_of_parameter
+    assert np.array(phi_mean).ndim == 2
+
+    # TODO: move the below to methods. Custom here for the time being.
+    mean = np.matmul(np.array(phi_mean)[:, :-1], np.atleast_1d(slag_ratio)) + np.array(phi_mean)[:,
+-1]  # assuming linear
+    dist = th.distributions.MultivariateNormal(loc=th.as_tensor(mean), covariance_matrix=th.tensor(phi_cov))
+    B1, B2, eta, E_act, Q_pot, T_ref = dist.sample()
+
+    return B1.item() * ureg('1/s'), B2.item() * ureg(''), eta.item() * ureg(''), \
+           E_act.item() * ureg('J/mol'), Q_pot.item() * ureg('J/kg'), Q_(T_ref.item(),ureg.degC)
+
+
+if __name__ == "__main__":
+    # test while developing this
+    # slight slope for for all except Q. Q follows the same relation as Erik had proposed
+    # Q_pot = Q_pot_max - (Q_pot_max - Q_pot_min) * slag_ratio
+    # also guessing 5% noise. sigma^2 = (5/100)^2 = 0.0025
+
+    phi_mean = [[0.01 * 2.916E-4, 2.916E-4], [0.01 * 0.0024229, 0.0024229], [0.01 * 5.554, 5.554],
+                [0.01 * 5653 * 8.3145, 5653 * 8.3145], [-200000, 300000], [0.01 * 25, 25]]
+    phi_cov = np.diag(0.0025 * np.array(phi_mean)[:, 1]).tolist()
+    seed = 7
+    print(dummy_hydration_parameters(0.5, phi_mean=phi_mean, phi_cov=phi_cov, seed=seed))