Merge pull request #198 from nasa/hotfix/v1.3.1

hotfix v1.3.1

README.md

@@ -38,7 +38,7 @@ Use the following to cite this repository:
     month     = May,
     year      = 2022,
     version   = {1.3.0},
-    url       = {https://github.com/nasa/prog_algs}
+    url       = {https://github.com/nasa/prog\_algs}
     }
 ```
 

requirements.txt

@@ -1,4 +1,4 @@
 filterpy>=1.4.5
 matplotlib>=3.5.2
 scipy>=1.7.3
-prog_models>=1.3
+prog_models>=1.3.1

setup.py

@@ -10,7 +10,7 @@
 
 setup(
     name = 'prog_algs',
-    version = '1.3.0',
+    version = '1.3.1',
     description = "The NASA Prognostics Algorithm Package is a framework for model-based prognostics (computation of remaining useful life) of engineering systems. It includes algorithms for state estimation and prediction, including uncertainty propagation. The algorithms use prognostic models (see prog_models) to perform estimation and prediction. The package enables rapid development of prognostics solutions for given models of components and systems. Algorithms can be swapped for comparative studies and evaluations",
     long_description=long_description,
     long_description_content_type='text/markdown',

src/prog_algs/__init__.py

@@ -5,7 +5,7 @@
 
 import warnings
 
-__version__ = '1.3.0'
+__version__ = '1.3.1'
 
 def run_prog_playback(obs, pred, future_loading, output_measurements, **kwargs):
     warnings.warn("Depreciated in 1.2.0, will be removed in a future release.", DeprecationWarning)

src/prog_algs/predictors/monte_carlo.py

@@ -40,9 +40,12 @@ class MonteCarlo(Predictor):
 
     def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs) -> PredictionResults:
         if isinstance(state, dict) or isinstance(state, self.model.StateContainer):
-            # Convert to UnweightedSamples
             from prog_algs.uncertain_data import ScalarData
-            state = ScalarData(state)
+            state = ScalarData(state, _type = self.model.StateContainer)
+        elif isinstance(state, UncertainData):
+            state._type = self.model.StateContainer
+        else:
+            raise TypeError("state must be UncertainData, dict, or StateContainer")
 
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
@@ -65,7 +68,7 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
 
         # Perform prediction
         for x in state:
-            events_remaining = deepcopy(params['events'])
+            events_remaining = params['events'].copy()
             first_output = ouput_eqn(x)
 
             time_of_event = {}
@@ -121,7 +124,7 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
                 t0 = times.pop()
                 inputs.pop()
                 x = states.pop()
-                last_state[event] = deepcopy(x)
+                last_state[event] = x.copy()
                 outputs.pop()
                 event_states.pop()
 
@@ -147,7 +150,7 @@ def predict(self, state : UncertainData, future_loading_eqn : Callable, **kwargs
 
         # Transform final states:
         last_states = {
-            key: UnweightedSamples([sample[key] for sample in last_states]) for key in time_of_event.keys()
+            key: UnweightedSamples([sample[key] for sample in last_states], _type = self.model.StateContainer) for key in time_of_event.keys()
         }
         time_of_event.final_state = last_states
 

src/prog_algs/predictors/predictor.py

@@ -40,7 +40,7 @@ def __init__(self, model, **kwargs):
         self.model = model
 
         self.parameters = deepcopy(self.default_parameters)
-        self.parameters['events'] = deepcopy(self.model.events)  # Events to predict to
+        self.parameters['events'] = self.model.events.copy()  # Events to predict to
         self.parameters.update(kwargs)
 
     @abstractmethod

src/prog_algs/predictors/unscented_transform.py

@@ -3,11 +3,10 @@
 from typing import Callable
 from .prediction import Prediction, UnweightedSamplesPrediction, PredictionResults
 from .predictor import Predictor
-from numpy import diag, array, transpose
+from numpy import diag, array, transpose, isnan
 from copy import deepcopy
-from math import isnan
 from filterpy import kalman
-from prog_algs.uncertain_data import MultivariateNormalDist
+from prog_algs.uncertain_data import MultivariateNormalDist, UncertainData, ScalarData
 
 
 class LazyUTPrediction(Prediction):
@@ -64,6 +63,8 @@ class UnscentedTransformPredictor(Predictor):
     ------------------------------
     alpha, beta, kappa: float
         UKF Scaling parameters. See: https://en.wikipedia.org/wiki/Kalman_filter#Unscented_Kalman_filter
+    Q: np.array
+        Process noise covariance matrix [nStates x nStates]
     t0 : float
         Initial time at which prediction begins, e.g., 0
     dt : float
@@ -105,24 +106,21 @@ def __init__(self, model, **kwargs):
         if 'Q' not in self.parameters:
             # Default 
             self.parameters['Q'] = diag([1.0e-1 for i in range(num_states)])
-        if 'R' not in self.parameters:
-            self.parameters['R'] = diag([1.0e-1 for i in range(num_measurements)])
 
         def measure(x):
-            x = {key: value for (key, value) in zip(self.__state_keys, x)}
+            x = model.StateContainer({key: value for (key, value) in zip(self.__state_keys, x)})
             z = model.output(x)
-            return {array(list(z.values()))}
+            return model.OutputContainer({array(list(z.values()))})
 
         def state_transition(x, dt):
-            x = {key: value for (key, value) in zip(self.__state_keys, x)}
+            x = model.StateContainer({key: value for (key, value) in zip(self.__state_keys, x)})
             x = model.next_state(x, self.__input, dt)
             x = model.apply_limits(x)
             return array(list(x.values()))
 
         self.sigma_points = kalman.MerweScaledSigmaPoints(num_states, alpha=self.parameters['alpha'], beta=self.parameters['beta'], kappa=self.parameters['kappa'])
         self.filter = kalman.UnscentedKalmanFilter(num_states, num_measurements, self.parameters['dt'], measure, state_transition, self.sigma_points)
         self.filter.Q = self.parameters['Q']
-        self.filter.R = self.parameters['R']
 
     def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionResults:
         """
@@ -161,6 +159,13 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
             Estimated time where a predicted event will occur for each sample. Note: Mean and Covariance Matrix will both 
             be nan if every sigma point doesnt reach threshold within horizon
         """
+        if isinstance(state, dict) or isinstance(state, self.model.StateContainer) or isinstance(state, ScalarData):
+            raise TypeError("state must be a distribution (e.g., MultivariateNormalDist, UnweightedSamples), not scalar")
+        elif isinstance(state, UncertainData):
+            state._type = self.model.StateContainer
+        else:
+            raise TypeError("state must be UncertainData, dict, or StateContainer")
+
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
         events_to_predict = params['events']
@@ -172,6 +177,7 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
         sigma_points = self.sigma_points
         n_points = sigma_points.num_sigmas()
         threshold_met = model.threshold_met
+        StateContainer = model.StateContainer
 
         # Update State 
         self.__state_keys = state_keys = state.mean.keys()  # Used to maintain ordering as we strip keys and return
@@ -194,18 +200,16 @@ def predict(self, state, future_loading_eqn : Callable, **kwargs) -> PredictionR
         def update_all():
             times.append(t)
             inputs.append(deepcopy(self.__input))  # Avoid optimization where u is not copied
-            x_dict = MultivariateNormalDist(self.__state_keys, filt.x, filt.P)
+            x_dict = MultivariateNormalDist(self.__state_keys, filt.x, filt.P, _type = self.model.StateContainer)
             states.append(x_dict)  # Avoid optimization where x is not copied
 
-        # Optimization
-        state_keys = self.__state_keys
-
         # Simulation
+        self.__input = future_loading_eqn(t, state.mean)
         update_all()  # First State
         while t < params['horizon']:
             # Iterate through time
             t += dt
-            mean_state = {key: x for (key, x) in zip(state_keys, filt.x)}
+            mean_state = StateContainer({key: x for (key, x) in zip(state_keys, filt.x)})
             self.__input = future_loading_eqn(t, mean_state)
             filt.predict(dt=dt)
 
@@ -221,7 +225,7 @@ def update_all():
             points = sigma_points.sigma_points(filt.x, filt.P)
             all_failed = True
             for i, point in zip(range(n_points), points):
-                x = {key: x for (key, x) in zip(state_keys, point)}
+                x = StateContainer({key: x for (key, x) in zip(state_keys, point)})
                 t_met = threshold_met(x)
 
                 # Check Thresholds
@@ -230,7 +234,7 @@ def update_all():
                         if isnan(ToE[key][i]):
                             # First time event has been reached
                             ToE[key][i] = t
-                            last_state[key][i] = deepcopy(x)
+                            last_state[key][i] = x.copy()
                     else:
                         all_failed = False  # This event for this sigma point hasn't been met yet
             if all_failed:
@@ -248,7 +252,7 @@ def update_all():
             last_state_pts = array([[last_state_i[state_key] for state_key in state_keys] for last_state_i in last_state[event_key]])
             # last_state_pts = transpose(last_state_pts)
             last_state_mean, last_state_cov = kalman.unscented_transform(last_state_pts, sigma_points.Wm, sigma_points.Wc)
-            final_state[event_key] = MultivariateNormalDist(state_keys, last_state_mean, last_state_cov)
+            final_state[event_key] = MultivariateNormalDist(state_keys, last_state_mean, last_state_cov, _type = self.model.StateContainer)
 
         # At this point only time of event, inputs, and state are calculated 
         inputs_prediction = UnweightedSamplesPrediction(times, [inputs])

src/prog_algs/state_estimators/kalman_filter.py

@@ -160,4 +160,4 @@ def x(self) -> MultivariateNormalDist:
         -------
         state = observer.x
         """
-        return MultivariateNormalDist(self.model.states, self.filter.x.ravel(), self.filter.P)
+        return MultivariateNormalDist(self.model.states, self.filter.x.ravel(), self.filter.P, _type = self.model.StateContainer)

src/prog_algs/state_estimators/particle_filter.py

@@ -49,6 +49,7 @@ def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
         super().__init__(model, x0, measurement_eqn = measurement_eqn, **kwargs)
 
         if measurement_eqn:
+            warn("Warning: measurement_eqn depreciated as of v1.3.1, will be removed in v1.4. Use Model subclassing instead. See examples.measurement_eqn_example")
             # update output_container
             from prog_models.utils.containers import DictLikeMatrixWrapper
             z0 = measurement_eqn(x0)
@@ -72,7 +73,7 @@ def __measure(x):
             sample_gen = x0.sample(self.parameters['num_particles'])
             samples = [array(sample_gen.key(k)) for k in x0.keys()]
         elif paramters_x0_exist and (parameters_x0_dict or parameters_x0_num):
-            warn("Warning: Use UncertainData type if estimating filtering with uncertain data.")
+            warn("Warning: x0_uncertainty depreciated as of v1.3, will be removed in v1.4. Use UncertainData type if estimating filtering with uncertain data.")
             x = array(list(x0.values()))
             if parameters_x0_dict:
                 sd = array([self.parameters['x0_uncertainty'][key] for key in x0.keys()])
@@ -104,11 +105,11 @@ def estimate(self, t : float, u, z):
         next_state = self.model.next_state
         apply_process_noise = self.model.apply_process_noise
         output = self._measure
-        apply_measurement_noise = self.model.apply_measurement_noise
+        # apply_measurement_noise = self.model.apply_measurement_noise
         noise_params = self.model.parameters['measurement_noise']
         num_particles = self.parameters['num_particles']
         # Check which output keys are present (i.e., output of measurement function)
-        measurement_keys = output({key: particles[key][0] for key in particles.keys()}).keys()
+        measurement_keys = output(self.model.StateContainer({key: particles[key][0] for key in particles.keys()})).keys()
         zPredicted = {key: empty(num_particles) for key in measurement_keys}
 
         if self.model.is_vectorized:
@@ -120,7 +121,7 @@ def estimate(self, t : float, u, z):
         else:
             # Propogate and calculate weights
             for i in range(num_particles):
-                x = {key: particles[key][i] for key in particles.keys()}
+                x = self.model.StateContainer({key: particles[key][i] for key in particles.keys()})
                 x = next_state(x, u, dt) 
                 x = apply_process_noise(x, dt)
                 for key in particles.keys():
@@ -172,4 +173,4 @@ def x(self) -> UnweightedSamples:
         -------
         state = observer.x
         """
-        return UnweightedSamples(self.particles)
+        return UnweightedSamples(self.particles, _type = self.model.StateContainer)

src/prog_algs/state_estimators/state_estimator.py

@@ -40,7 +40,7 @@ def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
             raise ProgAlgTypeError("model must have `outputs` property")
         if not hasattr(model, 'states'):
             raise ProgAlgTypeError("model must have `states` property")
-        self.model = deepcopy(model)
+        self.model = model
 
         if measurement_eqn is not None:
             warn('Measurement_eqn was deprecated in v1.3 in favor of model subclassing. I will remove this in v1.4. See measurement_equation example for more information', DeprecationWarning)

src/prog_algs/state_estimators/unscented_kalman_filter.py

@@ -52,27 +52,26 @@ def __init__(self, model, x0, measurement_eqn : Callable = None, **kwargs):
 
         if measurement_eqn is None: 
             def measure(x):
-                x = {key: value for (key, value) in zip(x0.keys(), x)}
+                x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
                 R_err = model.parameters['measurement_noise'].copy()
                 model.parameters['measurement_noise'] = dict.fromkeys(R_err, 0)
                 z = model.output(x)
                 model.parameters['measurement_noise'] = R_err
                 return array(list(z.values())).ravel()
         else:
+            warn("Warning: measurement_eqn depreciated as of v1.3.1, will be removed in v1.4. Use Model subclassing instead. See examples.measurement_eqn_example")
             def measure(x):
-                x = {key: value for (key, value) in zip(x0.keys(), x)}
+                x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
                 z = measurement_eqn(x)
                 return array(list(z.values())).ravel()
 
         if 'Q' not in self.parameters:
             self.parameters['Q'] = diag([1.0e-3 for i in x0.keys()])
 
         def state_transition(x, dt):
-            x = {key: value for (key, value) in zip(x0.keys(), x)}
+            x = model.StateContainer({key: value for (key, value) in zip(x0.keys(), x)})
             Q_err = model.parameters['process_noise'].copy()
             model.parameters['process_noise'] = dict.fromkeys(Q_err, 0)
-            z = model.output(x)
-            model.parameters['process_noise'] = Q_err
             x = model.next_state(x, self.__input, dt)
             return array(list(x.values())).ravel()
 
@@ -82,7 +81,7 @@ def state_transition(x, dt):
         self.filter = kalman.UnscentedKalmanFilter(num_states, num_measurements, self.parameters['dt'], measure, state_transition, points)
 
         if isinstance(x0, dict) or isinstance(x0, model.StateContainer):
-            warn("Warning: Use UncertainData type if estimating filtering with uncertain data.")
+            warn("Warning: x0_uncertainty depreciated as of v1.3, will be removed in v1.4. Use UncertainData type if estimating filtering with uncertain data.")
             self.filter.x = array(list(x0.values()))
             self.filter.P = self.parameters['Q'] / 10
         elif isinstance(x0, UncertainData):
@@ -93,9 +92,9 @@ def state_transition(x, dt):
             raise TypeError("TypeError: x0 initial state must be of type {{dict, UncertainData}}")
 
         if 'R' not in self.parameters:
-                # Size of what's being measured (not output) 
-                # This is determined by running the measure function on the first state
-                self.parameters['R'] = diag([1.0e-3 for i in range(len(measure(self.filter.x)))])
+            # Size of what's being measured (not output) 
+            # This is determined by running the measure function on the first state
+            self.parameters['R'] = diag([1.0e-3 for i in range(len(measure(self.filter.x)))])
         self.filter.Q = self.parameters['Q']
         self.filter.R = self.parameters['R']
 
@@ -131,4 +130,4 @@ def x(self) -> MultivariateNormalDist:
         -------
         state = observer.x
         """
-        return MultivariateNormalDist(self.x0.keys(), self.filter.x, self.filter.P)
+        return MultivariateNormalDist(self.x0.keys(), self.filter.x, self.filter.P, _type = self.model.StateContainer)

src/prog_algs/uncertain_data/multivariate_normal_dist.py

@@ -14,13 +14,17 @@ class MultivariateNormalDist(UncertainData):
             mean (array[float]): Mean values for state in the same order as labels
             covar (array[array[float]]): Covariance matrix for state
     """
-    def __init__(self, labels, mean: array, covar : array):
+    def __init__(self, labels, mean: array, covar : array, _type = dict):
         self.__labels = list(labels)
         self.__mean = array(list(mean))
         self.__covar = array(list(covar))
+        super().__init__(_type)
+
+    def __reduce__(self):
+        return (MultivariateNormalDist, (self.__labels, self.__mean, self.__covar))
 
     def __eq__(self, other : "MultivariateNormalDist") -> bool:
-        return self.keys() == other.keys() and self.mean == other.mean and (self.cov == other.cov).all()
+        return isinstance(other, MultivariateNormalDist) and self.keys() == other.keys() and self.mean == other.mean and (self.cov == other.cov).all()
 
     def __add__(self, other : int) -> "UncertainData":
         if other == 0:
@@ -58,7 +62,7 @@ def sample(self, num_samples : int = 1) -> UnweightedSamples:
 
         samples = multivariate_normal(self.__mean, self.__covar, num_samples)
         samples = [{key: value for (key, value) in zip(self.__labels, x)} for x in samples]
-        return UnweightedSamples(samples)
+        return UnweightedSamples(samples, _type = self._type)
 
     def keys(self) -> list:
         return self.__labels
@@ -69,8 +73,8 @@ def median(self) -> float:
         return self.mean
 
     @property
-    def mean(self) -> array:
-        return {key: value for (key, value) in zip(self.__labels, self.__mean)}
+    def mean(self) -> dict:
+        return self._type({key: value for (key, value) in zip(self.__labels, self.__mean)})
 
     def __str__(self) -> str:
         return 'MultivariateNormalDist(mean: {}, covar: {})'.format(self.__mean, self.__covar)