Add autoscale enum and add alternative options

Author: oerc0122

MDMC/control/control.py

@@ -17,7 +17,7 @@
 from MDMC.MD.engine_facades.facade import MDEngineError
 from MDMC.MD.parameters import Parameters
 from MDMC.MD.simulation import Simulation
-from MDMC.refinement.FoM.FoM_abs import ObservablePair
+from MDMC.refinement.FoM.FoM_abs import AutoScale, ObservablePair
 from MDMC.refinement.FoM.FoM_factory import FoMFactory
 from MDMC.refinement.minimizers.minimizer_factory import MinimizerFactory
 from MDMC.resolution.resolution_factory import ResolutionFactory
@@ -64,10 +64,8 @@ class Control:
           - ``rescale_factor`` (`float`, optional, defaults to `1.`) applied to
             the experimental data when calculating the FoM to ensure it is on
             the same scale as the calculated observable
-          - ``auto_scale`` (`bool`, optional, defaults to `False`) set the
-            ``rescale_factor`` automatically to minimise the FoM, if both
-            ``rescale_factor`` and ``auto_scale`` are provided then a warning
-            is printed and ``auto_scale`` takes precedence
+          - ``auto_scale`` (str or :class:`AutoScale`, optional, defaults to `CONSTANT`)
+            See :class:`AutoScale` documentation for more information.
           - ``use_FFT`` (`bool`, optional, defaults to `True`) whether to use
             Fast Fourier Transforms in the calculation of dependent variables.
             FFT speeds up calculation but places restrictions on spacing in the
@@ -161,7 +159,7 @@ class Control:
           'reader':'GENERIC_READER',
           'weight':0.5,
           'resolution':{'gaussian':2.35}
-          'auto_scale':True}]
+          'auto_scale': 'minimise_fom'}]
 
     Attributes
     ----------
@@ -253,39 +251,45 @@ def __init__(self, simulation: Simulation, exp_datasets: List[dict],
 
             self._validate_energy(MD_observable)
 
-            auto_scale = dset.get('auto_scale', False)
-            rescale_factor = dset.get('rescale_factor')
-            if auto_scale and rescale_factor and self.verbose != -1:
-                print('Both `rescale_factor` and `auto_scale` set for file {};'
-                      ' scaling will be automated to minimise FoM'
-                      ''.format(dset['file_name']))
-                rescale_factor = 1.
-            elif not rescale_factor:
-                rescale_factor = 1.
-
-            observable_pair = ObservablePair(exp_observable,
-                                             MD_observable,
-                                             dset['weight'],
-                                             rescale_factor=rescale_factor,
-                                             auto_scale=auto_scale)
+            auto_scale = dset.get("auto_scale", "CONSTANT")
+            rescale_factor = dset.get("rescale_factor", 1)
+
+            if auto_scale == "NONE" and rescale_factor != 1:
+                logging.warning("Requested no scaling, but provided `rescale_factor`. "
+                                "Disabling scaling.")
+                rescale_factor = 1
+
+            auto_scale = (AutoScale[auto_scale.upper()]
+                          if isinstance(auto_scale, str) else
+                          AutoScale(auto_scale))
+
+            if (
+                    rescale_factor != 1 and
+                    auto_scale is not AutoScale.CONSTANT and
+                    self.verbose != -1
+            ):
+                logging.warning("`rescale_factor` incompatible with %s set for file"
+                                " %s; scaling will be %s.",
+                                auto_scale.name, dset['file_name'], auto_scale.name)
+
+            observable_pair = ObservablePair(
+                exp_observable,
+                MD_observable,
+                dset['weight'],
+                rescale_factor=float(rescale_factor),
+                auto_scale=auto_scale,
+            )
             self.observable_pairs.append(observable_pair)
             self.recreated_independent_vars  = {}
             self.production_time_step = self.simulation.time_step
 
             # Take the largest minimum number of MD_steps needed by any dataset
-            min_MD_steps_dset = self._calculate_minimum_MD_steps(
-                observable_pair)
+            min_MD_steps_dset = self._calculate_minimum_MD_steps(observable_pair)
             minimum_MD_steps = max(minimum_MD_steps, min_MD_steps_dset)
 
-        if FoM_options is None or FoM_options.get('error') is None:
-            FoM_error = 'exp'
-        else:
-            FoM_error = FoM_options.get('error')
-
-        if FoM_options is None or FoM_options.get('norm') is None:
-            FoM_norm = 'data_points'
-        else:
-            FoM_norm = FoM_options.get('norm')
+        FoM_options = FoM_options if FoM_options is not None else {}
+        FoM_error = FoM_options.get("error", "exp")
+        FoM_norm = FoM_options.get("norm", "data_points")
 
         self.FoM_calculator = FoMFactory.create_FoM(FoM_error, self.observable_pairs,
                                                     norm=FoM_norm,

MDMC/refinement/FoM/ChiSquared_experror.py

@@ -14,6 +14,43 @@ class ChiSquaredExpError(FigureOfMerit):
     mathematical details.
     """
 
+    def _compute_unreduced(self, obs_pair: ObservablePair):
+        """
+        Compute the unreduced FoM value for the given observable pair.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated.
+
+        Returns
+        -------
+        float
+            Unreduced FoM value.
+        """
+        return np.sum((obs_pair.calculate_difference() /
+                       obs_pair.calculate_exp_errors()) ** 2)
+
+    def _minimise_factor(self, obs_pair: ObservablePair) -> float:
+        """
+        Minimise the FoM factor for the given FoM type.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated
+
+        Returns
+        -------
+        float
+            Computed auto_scale factor to minimise the FoM.
+        """
+        exp_errors = np.array(*obs_pair.exp_obs.errors.values())
+        exp_values = np.array(*obs_pair.exp_obs.dependent_variables.values())
+        MD_values = np.array(*obs_pair.MD_obs.dependent_variables.values())
+        return (np.sum((MD_values / exp_errors) ** 2) /
+                np.sum(MD_values * exp_values / exp_errors ** 2))
+
     def calculate_single_FoM(self, obs_pair: ObservablePair):
         """
         Calculates the chi-squared figure of merit for a single
@@ -30,16 +67,7 @@ def calculate_single_FoM(self, obs_pair: ObservablePair):
         float
             The FoM for the obs_pair
         """
-
-        if obs_pair.auto_scale:
-            exp_errors = np.array(*obs_pair.exp_obs.errors.values())
-            exp_values = np.array(
-                *obs_pair.exp_obs.dependent_variables.values())
-            MD_values = np.array(*obs_pair.MD_obs.dependent_variables.values())
-            obs_pair.rescale_factor = (np.sum((MD_values / exp_errors) ** 2) / np.sum(
-                MD_values * exp_values / exp_errors ** 2))
-
-        norm_factor = self.data_norm_factor(obs_pair=obs_pair)
-        value_unreduced = np.sum((obs_pair.calculate_difference()
-                                  / obs_pair.calculate_exp_errors()) ** 2)
+        obs_pair.rescale_factor = self.compute_rescale_factor(obs_pair)
+        norm_factor = self.data_norm_factor(obs_pair)
+        value_unreduced = self._compute_unreduced(obs_pair)
         return obs_pair.weight * value_unreduced / norm_factor

MDMC/refinement/FoM/FoM_abs.py

@@ -1,13 +1,59 @@
 """A module for Figure of Merits"""
 
 from abc import ABC, abstractmethod
+from enum import Enum, auto
 
 import numpy as np
 
 from MDMC.common.decorators import repr_decorator
 from MDMC.trajectory_analysis.observables.obs import Observable
 
 
+class AutoScale(Enum):
+    """
+    Autoscaling methods.
+
+    Notes
+    -----
+
+    - ``CONSTANT`` applies a constant scaling factor to the MD data
+      to align with the epxerimental data.
+
+    - ``MINIMISE_FOM`` uses the algorithm described in
+      :ref:`explanation/figure-of-merit:rescaling` to best minimise
+      the figure of merit.
+
+    - ``MATCH_MAXIMUM`` rescales the MD data such that the maximum is
+      the same as that of the experimental data.
+
+    - ``MATCH_ABS_MAXIMUM`` rescales the MD data such that the absolute maximum is
+      the same as that of the experimental data.
+
+    - ``MATCH_SUM`` is a naïve approximation to rescale the MD data such that
+      the integral area under the curve is the same under the assumumption
+      that the samples are evenly spaced.
+
+    - ``MATCH_ABS_SUM`` is like ``MATCH_SUM`` except that the data's
+      absolute values are used instead.
+    """
+    #: Constant valued scaling.
+    CONSTANT = auto()
+    #: Minimise the main FoM factor e.g. χ²
+    MINIMISE_FOM = auto()
+    #: Match the maxima of the values.
+    MATCH_MAXIMUM = auto()
+    #: Match the absolute maxima of the values.
+    MATCH_ABS_MAXIMUM = auto()
+    #: Match the integral area (approximated by a sum) of the values.
+    MATCH_SUM = auto()
+    #: Match the absolute integral area (approximated by an abs->sum) of the values.
+    MATCH_ABS_SUM = auto()
+
+    #: No scaling to be applied (CONSTANT w/ 1.).
+    NONE = CONSTANT
+    MINIMIZE_FOM = MINIMISE_FOM
+
+
 @repr_decorator('weight', 'exp_obs', 'MD_obs', 'rescale_factor', 'auto_scale')
 class ObservablePair:
 
@@ -35,7 +81,7 @@ class ObservablePair:
     """
 
     def __init__(self, exp_obs: Observable, MD_obs: Observable, weight: float,
-                 rescale_factor: float = 1., auto_scale: bool = False):
+                 rescale_factor: float = 1., auto_scale: AutoScale = AutoScale.NONE):
 
         self.exp_obs = exp_obs
         self.MD_obs = MD_obs
@@ -326,7 +372,6 @@ def calculate_exp_errors(self) -> np.ndarray:
 
         return np.array(*self.exp_obs.errors.values()) * self.rescale_factor
 
-
 @repr_decorator('value', 'obs_pairs')
 class FigureOfMerit(ABC):
 
@@ -421,8 +466,7 @@ def data_norm_factor(self, obs_pair: ObservablePair) -> int:
         """
 
         if self.norm:
-            norm_factor = np.size(
-                *obs_pair.MD_obs.dependent_variables.values())
+            norm_factor = np.size(*obs_pair.MD_obs.dependent_variables.values())
             norm_factor -= self.n_parameters
         else:
             norm_factor = 1
@@ -446,3 +490,69 @@ def calculate_single_FoM(self, obs_pair: ObservablePair) -> float:
         """
 
         raise NotImplementedError
+
+    @abstractmethod
+    def _compute_unreduced(self, obs_pair: ObservablePair) -> float:
+        """
+        Compute the unreduced FoM value for the given observable pair.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated.
+
+        Returns
+        -------
+        float
+            Unreduced FoM value.
+        """
+
+    @abstractmethod
+    def _minimise_factor(self, obs_pair: ObservablePair) -> float:
+        """
+        Minimise the FoM factor for the given FoM type.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated
+
+        Returns
+        -------
+        float
+            Computed auto_scale factor to minimise the FoM.
+        """
+
+    def compute_rescale_factor(self, obs_pair: ObservablePair) -> float:
+        """
+        Compute rescale factor for calculated observable to match experimental data.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated
+
+        Returns
+        -------
+        float
+            Computed rescale factor.
+        """
+        dep_vars = np.array(*obs_pair.exp_obs.dependent_variables.values())
+
+        match obs_pair.auto_scale:
+            case AutoScale.CONSTANT:
+                fac = obs_pair.rescale_factor
+            case AutoScale.MINIMISE_FOM:
+                fac = self._minimise_factor(obs_pair)
+            case AutoScale.MATCH_MAXIMUM:
+                fac = max(obs.max() for obs in dep_vars)
+            case AutoScale.MATCH_ABS_MAXIMUM:
+                fac = max(np.abs(obs).max() for obs in dep_vars)
+            case AutoScale.MATCH_SUM:
+                fac = sum(obs.sum() for obs in dep_vars)
+            case AutoScale.MATCH_ABS_SUM:
+                fac = sum(np.abs(obs.sum()) for obs in dep_vars)
+            case _:
+                fac = 1.
+
+        return fac

MDMC/refinement/FoM/RSquared_noneerror.py

@@ -33,6 +33,43 @@ class RSquared_noneerror(FigureOfMerit):
     simple linear scaling.
     """
 
+    def _compute_unreduced(self, obs_pair: ObservablePair):
+        """
+        Compute the unreduced FoM value for the given observable pair.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated.
+
+        Returns
+        -------
+        float
+            Unreduced FoM value.
+        """
+        return np.sum(obs_pair.calculate_difference() ** 2)
+
+    def _minimise_factor(self, obs_pair: ObservablePair) -> float:
+        """
+        Minimise the FoM factor for the given FoM type.
+
+        Parameters
+        ----------
+        obs_pair : ObservablePair
+            An ``ObservablePair`` for which the FoM is calculated
+
+        Returns
+        -------
+        float
+            Computed auto_scale factor to minimise the FoM.
+        """
+        exp_values = np.array(
+            *obs_pair.exp_obs.dependent_variables.values())
+        MD_values = np.array(*obs_pair.MD_obs.dependent_variables.values())
+        A = np.sum(MD_values * exp_values)
+        B = np.sum(exp_values ** 2)
+        return A / B
+
     def calculate_single_FoM(self, obs_pair: ObservablePair):
         # ignore line too long linting as it is necessary for LaTeX formatting
         # pylint: disable=line-too-long
@@ -68,14 +105,7 @@ def calculate_single_FoM(self, obs_pair: ObservablePair):
             The FoM for the obs_pair
         """
 
-        if obs_pair.auto_scale:
-            exp_values = np.array(
-                *obs_pair.exp_obs.dependent_variables.values())
-            MD_values = np.array(*obs_pair.MD_obs.dependent_variables.values())
-            A = np.sum(MD_values * exp_values)
-            B = np.sum(exp_values ** 2)
-            obs_pair.rescale_factor = A / B
-
-        norm_factor = self.data_norm_factor(obs_pair=obs_pair)
-        value_unreduced = np.sum(obs_pair.calculate_difference() ** 2)
+        obs_pair.rescale_factor = self.compute_rescale_factor(obs_pair)
+        norm_factor = self.data_norm_factor(obs_pair)
+        value_unreduced = self._compute_unreduced(obs_pair)
         return obs_pair.weight * value_unreduced / norm_factor

MDMC/refinement/FoM/__init__.py

@@ -1,9 +1,11 @@
 """A module for Figure of Merit calculation"""
 from . import ChiSquared_experror, FoM_abs, FoM_factory, RSquared_noneerror
+from .FoM_abs import AutoScale
 
 __all__ = [
     "ChiSquared_experror",
     "FoM_abs",
     "FoM_factory",
     "RSquared_noneerror",
+    "AutoScale",
 ]