43 improvement of sensitivity analysis

.coveragerc

@@ -15,4 +15,5 @@ exclude_lines =
 	if __name__ == '__main__':
 	raise ImportError
 	except Exception as error:
-	except ImportError
+	except ImportError
+	plt.show()

.gitignore

@@ -5,7 +5,8 @@
 *.bak
 .project
 .pydevproject
-aixcalibuha/examples/testzone
+examples/testzone
+examples/results
 
 # pip install #
 #####################

.gitlab-ci.yml

@@ -9,7 +9,7 @@ stages:
 variables:
   COVERAGE_TYPE: "FMU"
   PAGES_BRANCH: master
-  EXLUDE_PYTHON: 37
+  EXCLUDE_PYTHON: 37
   GIT_REPO: RWTH-EBC/AixCaliBuHA
   PYTHON_PACKAGE_NAME: "aixcalibuha"
   PYTHON_VERSION: "registry.git.rwth-aachen.de/ebc/ebc_all/gitlab_ci/templates:python_3.9"

CHANGELOG.md

@@ -51,3 +51,12 @@
 
 - **v0.3.1**
    - Issue 41: Fix logging and add kwarg
+
+-  **v1.0.0**
+    - Issue 43: Improvement of sensitivity analysis
+    - Enables verbose sensitivity analysis and the reuse of simulations
+    - It is now possible to use verbose sensitivity analysis for an automatic selection of tuner parameters
+    - Enables multiprocessing for the entire sensitivity process
+    - Sensitivity analysis is now usable for large models and data
+    - Add time dependent sensitivity analysis
+    - Ends support vor python 3.7

aixcalibuha/__init__.py

@@ -4,5 +4,5 @@
 """
 from .data_types import CalibrationClass, TunerParas, Goals
 from .calibration import Calibrator, MultipleClassCalibrator
-from .sensitivity_analysis import SobolAnalyzer, MorrisAnalyzer
-__version__ = "0.3.1"
+from .sensitivity_analysis import SobolAnalyzer, MorrisAnalyzer, FASTAnalyzer, PAWNAnalyzer
+__version__ = "1.0.0"

aixcalibuha/calibration/calibrator.py

@@ -10,6 +10,7 @@
 from typing import Dict
 from copy import copy
 import numpy as np
+import pandas as pd
 from ebcpy import data_types, Optimizer
 from ebcpy.simulationapi import SimulationAPI
 from aixcalibuha.utils import visualizer, MaxIterationsReached
@@ -243,13 +244,13 @@ def obj(self, xk, *args):
             )
             return self.ret_val_on_error
 
-        total_res = self._kpi_and_logging_calculation(
+        total_res, unweighted_objective = self._kpi_and_logging_calculation(
             xk_descaled=xk_descaled,
             counter=self._counter,
             results=sim_target_data
         )
 
-        return total_res
+        return total_res, unweighted_objective
 
     def mp_obj(self, x, *args):
         # Initialize list for results
@@ -304,7 +305,7 @@ def mp_obj(self, x, *args):
             if result is None:
                 total_res_list[idx] = self.ret_val_on_error
                 continue
-            total_res = self._kpi_and_logging_calculation(
+            total_res, unweighted_objective = self._kpi_and_logging_calculation(
                 xk_descaled=xk_descaled_list[idx],
                 counter=self._counter,
                 results=result
@@ -370,7 +371,7 @@ def _kpi_and_logging_calculation(self, *, xk_descaled, counter, results):
                 f"of iterations {self.max_itercount} has been reached."
             )
 
-        return total_res
+        return total_res, unweighted_objective
 
     def calibrate(self, framework, method=None, **kwargs) -> dict:
         """
@@ -486,7 +487,7 @@ def save_results(self, parameter_values: dict, filename: str):
             with open(s_path, 'w') as json_file:
                 json.dump(parameter_values, json_file, indent=4)
 
-    def validate(self, validation_class: CalibrationClass, calibration_result: Dict):
+    def validate(self, validation_class: CalibrationClass, calibration_result: Dict, verbose=False):
         """
         Validate the given calibration class based on the given
         values for tuner_parameters.
@@ -501,10 +502,10 @@ def validate(self, validation_class: CalibrationClass, calibration_result: Dict)
         self.logger.log(f"Start validation of model: {self.sim_api.model_name} with "
                         f"framework-class {self.__class__.__name__}")
         # Use start-time of calibration class
+        self.calibration_class = validation_class
         start_time = self._apply_start_time_method(
             start_time=self.calibration_class.start_time
         )
-        self.calibration_class = validation_class
         old_tuner_paras = copy(self.calibration_class.tuner_paras)
         tuner_values = list(calibration_result.values())
         self.calibration_class.tuner_paras = TunerParas(
@@ -526,12 +527,29 @@ def validate(self, validation_class: CalibrationClass, calibration_result: Dict)
         # Scale the tuner parameters
         xk = self.tuner_paras.scale(tuner_values)
         # Evaluate objective
-        obj = self.obj(xk=xk)
+        obj, unweighted_objective = self.obj(xk=xk)
         self.logger.validation_callback_func(
             obj=obj
         )
         # Reset tuner_parameters to avoid unwanted behaviour
         self.calibration_class.tuner_paras = old_tuner_paras
+        if verbose:
+            weights = [1]
+            objectives = [obj]
+            goals = ['all']
+            for goal, val in unweighted_objective.items():
+                weights.append(val[0])
+                objectives.append(val[1])
+                goals.append(goal)
+            index = pd.MultiIndex.from_product(
+                [[validation_class.name], goals],
+                names=['Class', 'Goal']
+            )
+            obj_verbos = pd.DataFrame(
+                {'weight': weights, validation_class.goals.statistical_measure: objectives},
+                index=index
+            )
+            return obj_verbos
         return obj
 
     def _handle_error(self, error):

aixcalibuha/calibration/multi_class_calibrator.py

@@ -176,6 +176,8 @@ def calibrate(self, framework, method=None, **kwargs) -> dict:
         for cal_run in self._cal_history:
             for p_name in cal_run['res']['Parameters'].index:
                 parameter_values[p_name] = cal_run['res']['Parameters'][p_name]
+        for p_name, res_intersection in res_tuner.items():
+            parameter_values[p_name] = res_intersection
         self.save_results(parameter_values=parameter_values,
                           filename='MultiClassCalibrationResult')
 

aixcalibuha/data_types.py

@@ -213,10 +213,8 @@ def eval_difference(self, verbose=False, penaltyfactor=1):
             )
             # Apply penalty function
             _diff = _diff * penaltyfactor
-
-            _verbose_calculation[self._weightings[i]] = _diff
+            _verbose_calculation[goal_name] = (self._weightings[i], _diff)
             total_difference += self._weightings[i] * _diff
-
         if verbose:
             return total_difference, _verbose_calculation
         return total_difference

aixcalibuha/sensitivity_analysis/__init__.py

@@ -2,6 +2,8 @@
 This package contains classes to perform
 sensitivity analysis with.
 """
-from .sensitivity_analyzer import SenAnalyzer
+from .sensitivity_analyzer import SenAnalyzer, _del_duplicates, _rename_tuner_names
 from .sobol import SobolAnalyzer
 from .morris import MorrisAnalyzer
+from .fast import FASTAnalyzer
+from .pawn import PAWNAnalyzer

aixcalibuha/sensitivity_analysis/fast.py

@@ -0,0 +1,90 @@
+"""
+Adds the FASTAnalyzer to the available
+classes of sensitivity analysis.
+"""
+from SALib.sample import fast_sampler as fast
+from SALib.analyze import fast as analyze_fast
+import numpy as np
+from aixcalibuha.sensitivity_analysis import SenAnalyzer
+from aixcalibuha import CalibrationClass
+import warnings
+
+
+class FASTAnalyzer(SenAnalyzer):
+    """
+    Sobol method from SALib https://salib.readthedocs.io/en/latest/api.html#fast-fourier-amplitude-sensitivity-test
+    A variance-based method which can compute the sensitivity measures
+    'S1' and 'ST'.
+
+    Additional arguments:
+
+    :keyword int M:
+        Default 4, used for the fast-method
+    :keyword seed:
+        Used for the fast-method
+    """
+
+    def __init__(self, sim_api, **kwargs):
+        super().__init__(
+            sim_api=sim_api,
+            **kwargs)
+        # Set additional kwargs
+        self.M = kwargs.pop("M", 4)
+        self.seed = kwargs.pop("seed", None)
+
+    @property
+    def analysis_variables(self):
+        """The analysis variables of the FAST method"""
+        return ['S1', 'ST']
+
+    def analysis_function(self, x, y):
+        """
+        Use the SALib.analyze.fast method to analyze the simulation results.
+
+        :param np.array x:
+            placeholder for the `X` parameter of the morris method not used for sobol
+        :param np.array y:
+            The NumPy array containing the model outputs
+        :return:
+            returns the result of the SALib.analyze.fast method (from the documentation:
+            Returns a dictionary with keys 'S1' and 'ST', where each entry is a list of
+            size D (the number of parameters) containing the indices in the same order
+            as the parameter file.)
+        """
+        return analyze_fast.analyze(self.problem, y,
+                                    M=self.M)
+
+    def create_sampler_demand(self):
+        """
+        Function to create the sampler parameters for the fast method
+        """
+        return {'M': self.M}
+
+    def generate_samples(self):
+        """
+        Run the sampler for fast and return the results.
+
+        :return:
+            The list of samples generated as a NumPy array with one row per sample
+            and each row containing one value for each variable name in `problem['names']`.
+        :rtype: np.ndarray
+        """
+        return fast.sample(self.problem,
+                           N=self.num_samples,
+                           **self.create_sampler_demand())
+
+    def _get_res_dict(self, result: dict, cal_class: CalibrationClass, analysis_variable: str):
+        """
+        Convert the result object to a dict with the key
+        being the variable name and the value being the result
+        associated to self.analysis_variable.
+        """
+        names = self.create_problem(cal_class.tuner_paras)['names']
+        if result is None:
+            return {var_name: np.abs(res_val)
+                    for var_name, res_val in zip(names,
+                                                 np.zeros(len(names)))}
+        else:
+            return {var_name: np.abs(res_val)
+                    for var_name, res_val in zip(names,
+                                                 result[analysis_variable])}

aixcalibuha/sensitivity_analysis/morris.py

@@ -2,6 +2,7 @@
 Adds the MorrisAnalyzer to the available
 classes of sensitivity analysis.
 """
+import numpy as np
 from SALib.sample import morris
 from SALib.analyze import morris as analyze_morris
 from aixcalibuha.sensitivity_analysis import SenAnalyzer
@@ -10,10 +11,14 @@
 
 class MorrisAnalyzer(SenAnalyzer):
     """
+    Moris method from SALib https://salib.readthedocs.io/en/latest/api.html#method-of-morris
+    An elementary effects (O-a-T) method which computes the sensitivity
+    measures 'mu', 'mu_star' and 'sigma' with a confidence interval for mu_star.
+
     Additional arguments:
 
     :keyword int num_levels:
-        Default 4, used for the morris-method
+        Default num_samples, used for the morris-method
     :keyword optimal_trajectories:
         Used for the morris-method
     :keyword bool local_optimization:
@@ -24,14 +29,14 @@ def __init__(self, sim_api, **kwargs):
             sim_api=sim_api,
             **kwargs)
         # Set additional kwargs
-        self.num_levels = kwargs.pop("num_levels", 4)
+        self.num_levels = kwargs.pop("num_levels", self.num_samples)
         self.optimal_trajectories = kwargs.pop("optimal_trajectories", None)
         self.local_optimization = kwargs.pop("local_optimization", True)
 
     @property
     def analysis_variables(self):
         """The analysis variables of the sobol method"""
-        return ['mu_star', 'sigma', 'mu_star_conf']
+        return ['mu_star', 'mu', 'sigma', 'mu_star_conf']
 
     def analysis_function(self, x, y):
         """
@@ -71,10 +76,14 @@ def generate_samples(self):
                              N=self.num_samples,
                              **self.create_sampler_demand())
 
-    def _get_res_dict(self, result: dict, cal_class: CalibrationClass):
+    def _get_res_dict(self, result: dict, cal_class: CalibrationClass, analysis_variable: str):
         """
         Convert the result object to a dict with the key
         being the variable name and the value being the result
         associated to self.analysis_variable.
         """
-        return dict(zip(result['names'], result[self.analysis_variable]))
+        if result is None:
+            names = cal_class.tuner_paras.get_names()
+            return dict(zip(names, np.zeros(len(names))))
+        else:
+            return dict(zip(result['names'], result[analysis_variable]))