Merge pull request #176 from nasa/dev

Fix ThrownObject

README.md

@@ -38,16 +38,16 @@ Use the following to cite this repository:
 @misc{2021_nasa_prog_models,
     author    = {Christopher Teubert and Matteo Corbetta and Chetan Kulkarni and Matthew Daigle},
     title     = {Prognostics Models Python Package},
-    month     = Nov,
+    month     = Dec,
     year      = 2021,
-    version   = {1.2},
+    version   = {1.2.2},
     url       = {https://github.com/nasa/prog_models}
     }
 ```
 
 The corresponding reference should look like this:
 
-C. Teubert, M. Corbetta, C. Kulkarni, and M. Daigle, Prognostics Models Python Package, v1.2, Nov. 2021. URL https://github.com/nasa/prog_models.
+C. Teubert, M. Corbetta, C. Kulkarni, and M. Daigle, Prognostics Models Python Package, v1.2.2, Dec. 2021. URL https://github.com/nasa/prog_models.
 
 ## Acknowledgements
 The structure and algorithms of this package are strongly inspired by the [MATLAB Prognostics Model Library](https://github.com/nasa/PrognosticsModelLibrary). We would like to recognize Matthew Daigle and the rest of the team that contributed to the Prognostics Model Library for the contributions their work on the MATLAB library made to the design of prog_models

examples/future_loading.py

@@ -7,7 +7,6 @@
 
 from prog_models.models import BatteryCircuit
 from statistics import mean
-from prog_models.visualize import plot_timeseries
 from numpy.random import normal
 
 def run_example(): 
@@ -35,8 +34,8 @@ def future_loading(t, x=None):
     (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, **options)
 
     # Now lets plot the inputs and event_states
-    plot_timeseries(times, inputs, options={'ylabel': 'Variable Load Current (amps)'})
-    plot_timeseries(times, event_states, options={'ylabel': 'Variable Load Event State'})
+    inputs.plot(ylabel = 'Variable Load Current (amps)')
+    event_states.plot(ylabel = 'Variable Load Event State')
 
     ## Example 2: Moving Average loading 
     # This is useful in cases where you are running reoccuring simulations, and are measuring the actual load on the system, 
@@ -71,8 +70,8 @@ def moving_avg(i):
     (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, **options)
 
     # Now lets plot the inputs and event_states
-    plot_timeseries(times, inputs, options={'ylabel': 'Moving Average Current (amps)'})
-    plot_timeseries(times, event_states, options={'ylabel': 'Moving Average Event State'})
+    inputs.plot(ylabel = 'Moving Average Current (amps)')
+    event_states.plot(ylabel = 'Moving Average Event State')
 
     # In this case, this estimate is wrong because loading will not be steady, but at least it would give you an approximation.
 
@@ -101,8 +100,8 @@ def future_loading(t, x=None):
     (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, **options)
 
     # Now lets plot the inputs and event_states
-    plot_timeseries(times, inputs, options={'ylabel': 'Variable Gaussian Current (amps)'})
-    plot_timeseries(times, event_states, options={'ylabel': 'Variable Gaussian Event State'})
+    inputs.plot(ylabel = 'Variable Gaussian Current (amps)')
+    event_states.plot(ylabel = 'Variable Gaussian Event State')
 
     # Example 4: Gaussian- increasing with time
     # For this we're using moving average. This is realistic because the further out from current time you get, 
@@ -140,8 +139,8 @@ def moving_avg(i):
     (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, **options)
 
     # Now lets plot the inputs and event_states
-    plot_timeseries(times, inputs, options={'ylabel': 'Moving Average Current (amps)'})
-    plot_timeseries(times, event_states, options={'ylabel': 'Moving Average Event State'})
+    inputs.plot(ylabel = 'Moving Average Current (amps)')
+    event_states.plot(ylabel = 'Moving Average Event State')
 
     # In this example future_loading.t has to be updated with current time before each prediction.
 
@@ -162,8 +161,8 @@ def future_loading(t, x=None):
     (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_loading, **options)
 
     # Now lets plot the inputs and event_states
-    plot_timeseries(times, inputs, options={'ylabel': 'Moving Average Current (amps)'})
-    plot_timeseries(times, event_states, options={'ylabel': 'Moving Average Event State'})
+    inputs.plot(ylabel = 'Moving Average Current (amps)')
+    event_states.plot(ylabel = 'Moving Average Event State')
 
     # In this example future_loading.t has to be updated with current time before each prediction.
 

examples/noise.py

@@ -30,15 +30,15 @@ def future_load(t, x=None):
     print('\t- impact time: {}s'.format(times[-1]))
 
     # Ex3: noise- more noise on position than velocity
-    process_noise = {'x': 0.25, 'v': 0.75}
+    process_noise = {'x': 0.25, 'v': 0.75, 'max_x': 0}
     m = ThrownObject(process_noise = process_noise) 
     print('\nExample with more noise on position than velocity')
     (times, _, states, outputs, _) = m.simulate_to_threshold(future_load, threshold_keys=[event], dt=0.005, save_freq=1)
     print('\t- states: {}'.format(['{}s: {}'.format(round(t,2), x) for (t,x) in zip(times, states)])) 
     print('\t- impact time: {}s'.format(times[-1]))
 
     # Ex4: noise- Ex3 but uniform
-    process_noise = {'x': 0.25, 'v': 0.75}
+    process_noise = {'x': 0.25, 'v': 0.75, 'max_x': 0}
     process_noise_dist = 'uniform'
     model_config = {'process_noise_dist': process_noise_dist, 'process_noise': process_noise}
     m = ThrownObject(**model_config) 
@@ -48,7 +48,7 @@ def future_load(t, x=None):
     print('\t- impact time: {}s'.format(times[-1]))
 
     # Ex5: noise- Ex3 but triangle
-    process_noise = {'x': 0.25, 'v': 0.75}
+    process_noise = {'x': 0.25, 'v': 0.75, 'max_x': 1e-9}
     process_noise_dist = 'triangular'
     model_config = {'process_noise_dist': process_noise_dist, 'process_noise': process_noise}
     m = ThrownObject(**model_config) 
@@ -76,7 +76,8 @@ def future_load(t, x=None):
     def apply_proportional_process_noise(self, x, dt = 1):
         return {
             'x': x['x'],  # No noise on state
-            'v': x['v'] - dt*0.5*x['v']
+            'v': x['v'] - dt*0.5*x['v'],
+            'max_x': x['max_x']  # No noise on max_x
         }
     model_config = {'process_noise': apply_proportional_process_noise}
     m = ThrownObject(**model_config)

examples/vectorized.py

@@ -18,7 +18,8 @@ def future_load(t, x=None):
     # For this example we are saying there are 4 throwers of various strengths and heights
     first_state = {
         'x': array([1.75, 1.8, 1.85, 1.9]),
-        'v': array([35, 39, 22, 47])
+        'v': array([35, 39, 22, 47]),
+        'max_x': array([1.75, 1.8, 1.85, 1.9])
     }
 
     # Step 3: Simulate to threshold

setup.py

@@ -11,7 +11,7 @@
 
 setup(
     name = 'prog_models',
-    version = '1.2.1', #pkg_resources.require("prog_models")[0].version,
+    version = '1.2.2', #pkg_resources.require("prog_models")[0].version,
     description = 'The NASA Prognostic Model Package is a python modeling framework focused on defining and building models for prognostics (computation of remaining useful life) of engineering systems, and provides a set of prognostics models for select components developed within this framework, suitable for use in prognostics applications for these components.',
     long_description=long_description,
     long_description_content_type='text/markdown',
@@ -47,6 +47,7 @@
     license = 'NOSA',
     project_urls={  # Optional
         'Bug Reports': 'https://github.com/nasa/prog_models/issues',
+        'Docs': 'https://nasa.github.io/prog_models',
         'Organization': 'https://prognostics.nasa.gov/',
         'Source': 'https://github.com/nasa/prog_models',
     },

sphinx_config/index.rst

@@ -29,15 +29,15 @@ Use the following to cite this repository:
 @misc{2021_nasa_prog_model,
   | author    = {Christopher Teubert and Chetan Kulkarni and Matteo Corbetta and Matthew Daigle},
   | title     = {Prognostics Model Python Package},
-  | month     = Nov,
+  | month     = Dec,
   | year      = 2021,
-  | version   = {1.2.0},
+  | version   = {1.2.2},
   | url       = {https://github.com/nasa/prog_models}
   | }
 
 The corresponding reference should look like this:
 
-C. Teubert, C. Kulkarni, M. Corbetta, M. Daigle, Prognostics Model Python Package, v1.2, Nov. 2021. URL https://github.com/nasa/prog_models.
+C. Teubert, C. Kulkarni, M. Corbetta, M. Daigle, Prognostics Model Python Package, v1.2.2, Dec. 2021. URL https://github.com/nasa/prog_models.
 
 Indices and tables
 -----------------------
@@ -51,4 +51,4 @@ Disclaimers
 
 No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND, EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION, IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER, CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS, RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND DISTRIBUTES IT "AS IS."
 
-Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.  IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES, EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON, OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
+Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.  IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES, EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON, OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.

src/prog_models/__init__.py

@@ -4,4 +4,4 @@
 from .prognostics_model import PrognosticsModel
 from .exceptions import ProgModelException, ProgModelInputException, ProgModelTypeError
 
-__version__ = '1.2.1'
+__version__ = '1.2.2'

src/prog_models/models/thrown_object.py

@@ -8,19 +8,54 @@
 class ThrownObject(PrognosticsModel):
     """
     Model that similates an object thrown into the air without air resistance
+
+    Events (2)
+        | falling: The object is falling
+        | impact: The object has hit the ground
+
+    Inputs/Loading: (0)
+
+    States: (2)
+        | x: Position in space (m)
+        | v: Velocity in space (m/s)
+
+    Outputs/Measurements: (1)
+        | x: Position in space (m)
+
+    Keyword Args
+    ------------
+        process_noise : Optional, float or Dict[Srt, float]
+          Process noise (applied at dx/next_state). 
+          Can be number (e.g., .2) applied to every state, a dictionary of values for each 
+          state (e.g., {'x1': 0.2, 'x2': 0.3}), or a function (x) -> x
+        process_noise_dist : Optional, String
+          distribution for process noise (e.g., normal, uniform, triangular)
+        measurement_noise : Optional, float or Dict[Srt, float]
+          Measurement noise (applied in output eqn).
+          Can be number (e.g., .2) applied to every output, a dictionary of values for each
+          output (e.g., {'z1': 0.2, 'z2': 0.3}), or a function (z) -> z
+        measurement_noise_dist : Optional, String
+          distribution for measurement noise (e.g., normal, uniform, triangular)
+        g : Optional, float
+            Acceleration due to gravity (m/s^2). Default is 9.81 m/s^2 (standard gravity)
+        thrower_height : Optional, float
+            Height of the thrower (m). Default is 1.83 m
+        throwing_speed : Optional, float
+            Speed at which the ball is thrown (m/s). Default is 40 m/s
     """
 
     inputs = []  # no inputs, no way to control
     states = [
-        'x',  # Position (m) 
-        'v'  # Velocity (m/s)
+        'x',    # Position (m) 
+        'v',    # Velocity (m/s)
+        'max_x' # Maximum state
         ]
     outputs = [
-        'x'  # Position (m)
+        'x'     # Position (m)
     ]
     events = [
         'falling',  # Event- object is falling
-        'impact'  # Event- object has impacted ground
+        'impact'    # Event- object has impacted ground
     ]
 
     is_vectorized = True
@@ -40,12 +75,15 @@ def __init__(self, **kwargs):
     def initialize(self, u=None, z=None):
         return {
             'x': self.parameters['thrower_height'],  # Thrown, so initial altitude is height of thrower
-            'v': self.parameters['throwing_speed']  # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
+            'v': self.parameters['throwing_speed'],  # Velocity at which the ball is thrown - this guy is a professional baseball pitcher
+            'max_x': self.parameters['thrower_height']
             }
 
-    def dx(self, x, u):
-        return {'x': x['v'],
-                'v': self.parameters['g']}  # Acceleration of gravity
+    def next_state(self, x, u, dt):
+        next_x =  x['x'] + x['v']*dt
+        return {'x': next_x,
+                'v': x['v'] + self.parameters['g']*dt,  # Acceleration of gravity
+                'max_x': maximum(x['max_x'], next_x)}
 
     def output(self, x):
         return {'x': x['x']}
@@ -59,8 +97,8 @@ def threshold_met(self, x):
         }
 
     def event_state(self, x): 
-        self.max_x = maximum(self.max_x, x['x'])  # Maximum altitude
+        x['max_x'] = maximum(x['max_x'], x['x'])  # Maximum altitude
         return {
             'falling': maximum(x['v']/self.parameters['throwing_speed'],0),  # Throwing speed is max speed
-            'impact': maximum(x['x']/self.max_x,0)  # 1 until falling begins, then it's fraction of height
+            'impact': maximum(x['x']/x['max_x'],0)  # 1 until falling begins, then it's fraction of height
         }

src/prog_models/prognostics_model.py

@@ -7,7 +7,7 @@
 import numpy as np
 from copy import deepcopy
 from warnings import warn
-from collections import UserDict
+from collections import UserDict, abc
 import types
 from array import array
 from .sim_result import SimResult, LazySimResult
@@ -249,6 +249,12 @@ def __init__(self, **kwargs):
         except Exception:
             raise ProgModelTypeError('Model noise poorly configured')
 
+    def __eq__(self, other):
+        """
+        Check if two models are equal
+        """
+        return self.__class__ == other.__class__ and self.parameters == other.parameters
+
     def __str__(self):
         return "{} Prognostics Model (Events: {})".format(type(self).__name__, self.events)
 
@@ -728,25 +734,32 @@ def simulate_to_threshold(self, future_loading_eqn, first_output = None, thresho
         ----------
         future_loading_eqn : callable
             Function of (t) -> z used to predict future loading (output) at a given time (t)
+
+        Keyword Arguments
+        -----------------
+        t0 : Number, optional
+            Starting time for simulation in seconds (default: 0.0) \n
+        dt : Number or function, optional
+            time step (s), e.g. dt = 0.1 or function (t, x) -> dt\n
+        save_freq : Number, optional
+            Frequency at which output is saved (s), e.g., save_freq = 10 \n
+        save_pts : List[Number], optional
+            Additional ordered list of custom times where output is saved (s), e.g., save_pts= [50, 75] \n
+        horizon : Number, optional
+            maximum time that the model will be simulated forward (s), e.g., horizon = 1000 \n
         first_output : dict, optional
             First measured output, needed to initialize state for some classes. Can be omitted for classes that dont use this
-        threshold_keys: [str], optional
+        threshold_keys: List[str] or str, optional
             Keys for events that will trigger the end of simulation.
             If blank, simulation will occur if any event will be met ()
-        options: keyword arguments, optional
-            Configuration options for the simulation \n
-            Note: configuration of the model is set through model.parameters.\n
-            Supported parameters:\n
-             * t0 (Number) : Starting time for simulation in seconds (default: 0.0) \n
-             * dt (Number or function): time step (s), e.g. dt = 0.1 or function (t, x) -> dt\n
-             * save_freq (Number): Frequency at which output is saved (s), e.g., save_freq = 10 \n
-             * save_pts (List[Number]): Additional ordered list of custom times where output is saved (s), e.g., save_pts= [50, 75] \n
-             * horizon (Number): maximum time that the model will be simulated forward (s), e.g., horizon = 1000 \n
-             * x (dict): initial state dict, e.g., x= {'x1': 10, 'x2': -5.3}\n
-             * thresholds_met_eqn (function/lambda): custom equation to indicate logic for when to stop sim f(thresholds_met) -> bool\n
-             * print (bool): toggle intermediate printing, e.g., print = True\n
+        x : dict, optional
+            initial state dict, e.g., x= {'x1': 10, 'x2': -5.3}\n
+        thresholds_met_eqn : function/lambda, optional
+            custom equation to indicate logic for when to stop sim f(thresholds_met) -> bool\n
+        print : bool, optional
+            toggle intermediate printing, e.g., print = True\n
             e.g., m.simulate_to_threshold(eqn, z, dt=0.1, save_pts=[1, 2])
-        
+    
         Returns
         -------
         times: Array[number]
@@ -778,13 +791,21 @@ def simulate_to_threshold(self, future_loading_eqn, first_output = None, thresho
         | first_output = {'o1': 3.2, 'o2': 1.2}
         | m = PrognosticsModel() # Replace with specific model being simulated
         | (times, inputs, states, outputs, event_states) = m.simulate_to_threshold(future_load_eqn, first_output)
+
+        Note
+        ----
+        configuration of the model is set through model.parameters.\n
         """
         # Input Validation
         if first_output and not all(key in first_output for key in self.outputs):
             raise ProgModelInputException("Missing key in 'first_output', must have every key in model.outputs")
 
         if not (callable(future_loading_eqn)):
             raise ProgModelInputException("'future_loading_eqn' must be callable f(t)")
+
+        if isinstance(threshold_keys, str):
+            # A single threshold key
+            threshold_keys = [threshold_keys]
 
         if threshold_keys and not all([key in self.events for key in threshold_keys]):
             raise ProgModelInputException("threshold_keys must be event names")
@@ -809,6 +830,8 @@ def simulate_to_threshold(self, future_loading_eqn, first_output = None, thresho
             raise ProgModelInputException("'save_freq' must be a number, was a {}".format(type(config['save_freq'])))
         if config['save_freq'] <= 0:
             raise ProgModelInputException("'save_freq' must be positive, was {}".format(config['save_freq']))
+        if not isinstance(config['save_pts'], abc.Iterable):
+            raise ProgModelInputException("'save_pts' must be list or array, was a {}".format(type(config['save_pts'])))
         if not isinstance(config['horizon'], Number):
             raise ProgModelInputException("'horizon' must be a number, was a {}".format(type(config['horizon'])))
         if config['horizon'] < 0: