Merge branch 'dev' into feature/particle_filter_fix

Author: teubert

examples/basic_example.py

@@ -20,9 +20,6 @@
 def run_example():
     # Step 1: Setup model & future loading
     m = ThrownObject(process_noise = 1)
-    def future_loading(t, x = None):
-        # No load for a thrown object
-        return m.InputContainer({})
     initial_state = m.initialize()
 
     # Step 2: Demonstrating state estimator
@@ -42,7 +39,7 @@ def future_loading(t, x = None):
     # Step 2c: Perform state estimation step, given some measurement, above what's expected
     example_measurements = m.OutputContainer({'x': 7.5})
     t = 0.1
-    u = future_loading(t)
+    u = m.InputContainer({})
     filt.estimate(t, u, example_measurements)  # Update state, given (example) sensor data
 
     # Step 2d: Print & Plot Resulting Posterior State
@@ -65,7 +62,7 @@ def future_loading(t, x = None):
     # Step 3b: Perform a prediction
     NUM_SAMPLES = 50
     STEP_SIZE = 0.01
-    mc_results = mc.predict(filt.x, future_loading, n_samples = NUM_SAMPLES, dt=STEP_SIZE, save_freq=STEP_SIZE)
+    mc_results = mc.predict(filt.x, n_samples = NUM_SAMPLES, dt=STEP_SIZE, save_freq=STEP_SIZE)
     print('Predicted time of event (ToE): ', mc_results.time_of_event.mean)
     # Here there are 2 events predicted, when the object starts falling, and when it impacts the ground.
 

examples/horizon.py

@@ -18,9 +18,7 @@
 
 def run_example():
     # Step 1: Setup model & future loading
-    def future_loading(t, x = None):
-        return {}
-    m = ThrownObject(process_noise = 0.25, measurement_noise = 0.2)
+    m = ThrownObject(process_noise=0.25, measurement_noise=0.2)
     initial_state = m.initialize()
 
     # Step 2: Demonstrating state estimator
@@ -53,7 +51,7 @@ def future_loading(t, x = None):
     PREDICTION_HORIZON = 7.75
     samples = filt.x  # Since we're using a particle filter, which is also sample-based, we can directly use the samples, without changes
     STEP_SIZE = 0.01
-    mc_results = mc.predict(samples, future_loading, dt=STEP_SIZE, horizon = PREDICTION_HORIZON)
+    mc_results = mc.predict(samples, dt=STEP_SIZE, horizon=PREDICTION_HORIZON)
     print("\nPredicted Time of Event:")
     metrics = mc_results.time_of_event.metrics()
     pprint(metrics)  # Note this takes some time

examples/predict_specific_event.py

@@ -12,8 +12,6 @@ def run_example():
     m = ThrownObject()
     initial_state = m.initialize()
     load = m.InputContainer({}) # Optimization - create once
-    def future_loading(t, x = None):
-        return load
 
     ## State Estimation - perform a single ukf state estimate step
     filt = state_estimators.UnscentedKalmanFilter(m, initial_state)
@@ -24,7 +22,7 @@ def future_loading(t, x = None):
     pred = predictors.UnscentedTransformPredictor(m)
 
     # Predict with a step size of 0.1
-    mc_results = pred.predict(filt.x, future_loading, dt=0.1, save_freq= 1, events=['impact'])
+    mc_results = pred.predict(filt.x, dt=0.1, save_freq= 1, events=['impact'])
 
     # Print Results
     for i, time in enumerate(mc_results.times):

examples/sensitivity.py

@@ -14,22 +14,18 @@ def run_example():
     # Step 1: Create instance of model
     m = ThrownObject()
 
-    # Step 2: Setup for simulation 
-    def future_load(t, x=None):
-        return m.InputContainer({})
-
-    # Step 3: Setup range on parameters considered
+    # Step 2: Setup range on parameters considered
     thrower_height_range = np.arange(1.2, 2.1, 0.1)
 
-    # Step 4: Sim for each 
+    # Step 3: Sim for each 
     event = 'impact'
     eods = np.empty(len(thrower_height_range))
     for (i, thrower_height) in zip(range(len(thrower_height_range)), thrower_height_range):
         m.parameters['thrower_height'] = thrower_height
-        simulated_results = m.simulate_to_threshold(future_load, threshold_keys=[event], dt =1e-3, save_freq =10)
+        simulated_results = m.simulate_to_threshold(threshold_keys=[event], dt =1e-3, save_freq =10)
         eods[i] = simulated_results.times[-1]
 
-    # Step 5: Analysis
+    # Step 4: Analysis
     print('For a reasonable range of heights, impact time is between {} and {}'.format(round(eods[0],3), round(eods[-1],3)))
     sensitivity = (eods[-1]-eods[0])/(thrower_height_range[-1] - thrower_height_range[0])
     print('  - Average sensitivity: {} s per cm height'.format(round(sensitivity/100, 6)))
@@ -40,7 +36,7 @@ def future_load(t, x=None):
     eods = np.empty(len(throw_speed_range))
     for (i, throw_speed) in zip(range(len(throw_speed_range)), throw_speed_range):
         m.parameters['throwing_speed'] = throw_speed
-        simulated_results = m.simulate_to_threshold(future_load, threshold_keys=[event], options={'dt':1e-3, 'save_freq':10})
+        simulated_results = m.simulate_to_threshold(threshold_keys=[event], options={'dt':1e-3, 'save_freq':10})
         eods[i] = simulated_results.times[-1]
 
     print('\nFor a reasonable range of throwing speeds, impact time is between {} and {}'.format(round(eods[0],3), round(eods[-1],3)))

examples/state_limits.py

@@ -15,11 +15,7 @@ def run_example():
     # Step 1: Create instance of model (without drag)
     m = ThrownObject( cd = 0 )
 
-    # Step 2: Setup for simulation 
-    def future_load(t, x=None):
-        return {}
-
-    # add state limits
+    # Step 2: add state limits
     m.state_limits = {
         # object may not go below ground height
         'x': (0, inf),
@@ -30,7 +26,7 @@ def future_load(t, x=None):
 
     # Step 3: Simulate to impact
     event = 'impact'
-    simulated_results = m.simulate_to_threshold(future_load, threshold_keys=[event], dt=0.005, save_freq=1)
+    simulated_results = m.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1)
 
     # Print states
     print('Example 1')
@@ -42,7 +38,7 @@ def future_load(t, x=None):
     x0 = m.initialize(u = {}, z = {})
     x0['x'] = -1
 
-    simulated_results = m.simulate_to_threshold(future_load, threshold_keys=[event], dt=0.005, save_freq=1, x = x0)
+    simulated_results = m.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=1, x=x0)
 
     # Print states
     print('Example 2')
@@ -57,7 +53,7 @@ def future_load(t, x=None):
     m.parameters['g'] = -50000000
 
     print('Example 3')
-    simulated_results = m.simulate_to_threshold(future_load, threshold_keys=[event], dt=0.005, save_freq=0.3, x = x0, print = True, progress = False)
+    simulated_results = m.simulate_to_threshold(threshold_keys=[event], dt=0.005, save_freq=0.3, x=x0, print=True, progress=False)
 
     # Note that the limits can also be applied manually using the apply_limits function
     print('limiting states')

src/progpy/predictors/monte_carlo.py

@@ -19,27 +19,56 @@ class MonteCarlo(Predictor):
     
     Configuration Parameters
     ------------------------------
-    t0 : float
-        Initial time at which prediction begins, e.g., 0
-    dt : float
-        Simulation step size (s), e.g., 0.1
-    events : list[str]
-        Events to predict (subset of model.events) e.g., ['event1', 'event2']
-    horizon : float
-        Prediction horizon (s)
-    n_samples : int
-        Number of samples to use. If not specified, a default value is used. If state is type UnweightedSamples and n_samples is not provided, the provided unweighted samples will be used directly.
-    save_freq : float
-        Frequency at which results are saved (s)
-    save_pts : list[float]
-        Any additional savepoints (s) e.g., [10.1, 22.5]
+    n_samples : int, optional
+        Default number of samples to use. If not specified, a default value is used. If state is type UnweightedSamples and n_samples is not provided, the provided unweighted samples will be used directly.
+    save_freq : float, optional
+        Default frequency at which results are saved (s).
     """
 
+    __DEFAULT_N_SAMPLES = 100 # Default number of samples to use, if none specified and not UncertainData
+
     default_parameters = { 
-        'n_samples': 100  # Default number of samples to use, if none specified
+        'n_samples': None
     }
 
-    def predict(self, state: UncertainData, future_loading_eqn: Callable, **kwargs) -> PredictionResults:
+    def predict(self, state: UncertainData, future_loading_eqn: Callable = None, **kwargs) -> PredictionResults:
+        """
+        Perform a single prediction
+
+        Parameters
+        ----------
+        state : UncertainData 
+            Distribution representing current state of the system
+        future_loading_eqn : function (t, x=None) -> z, optional
+            Function to generate an estimate of loading at future time t, and state x
+
+        Keyword Arguments
+        ------------------
+        t0 : float, optional
+            Initial time at which prediction begins, e.g., 0
+        dt : float, optional
+            Simulation step size (s), e.g., 0.1
+        events : list[str], optional
+            Events to predict (subset of model.events) e.g., ['event1', 'event2']
+        horizon : float, optional
+            Prediction horizon (s)
+        n_samples : int, optional
+            Number of samples to use. If not specified, a default value is used. If state is type UnweightedSamples and n_samples is not provided, the provided unweighted samples will be used directly.
+        save_freq : float, optional
+            Frequency at which results are saved (s)
+        save_pts : list[float], optional
+            Any additional savepoints (s) e.g., [10.1, 22.5]
+
+        Return
+        ----------
+        result from prediction, including: NameTuple
+            * times (List[float]): Times for each savepoint such that inputs.snapshot(i), states.snapshot(i), outputs.snapshot(i), and event_states.snapshot(i) are all at times[i]            
+            * inputs (Prediction): Inputs at each savepoint such that inputs.snapshot(i) is the input distribution (type UncertainData) at times[i]
+            * states (Prediction): States at each savepoint such that states.snapshot(i) is the state distribution (type UncertainData) at times[i]
+            * outputs (Prediction): Outputs at each savepoint such that outputs.snapshot(i) is the output distribution (type UncertainData) at times[i]
+            * event_states (Prediction): Event states at each savepoint such that event_states.snapshot(i) is the event state distribution (type UncertainData) at times[i]
+            * time_of_event (UncertainData): Distribution of predicted Time of Event (ToE) for each predicted event, represented by some subclass of UncertaintData (e.g., MultivariateNormalDist)
+        """
         if isinstance(state, dict) or isinstance(state, self.model.StateContainer):
             from progpy.uncertain_data import ScalarData
             state = ScalarData(state, _type = self.model.StateContainer)
@@ -48,16 +77,28 @@ def predict(self, state: UncertainData, future_loading_eqn: Callable, **kwargs)
         else:
             raise TypeError("state must be UncertainData, dict, or StateContainer")
 
+        if future_loading_eqn is None:
+            future_loading_eqn = lambda t, x=None: self.model.InputContainer({})
+
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
         params['print'] = False
         params['progress'] = False
 
+        if not isinstance(state, UnweightedSamples) and params['n_samples'] is None:
+            # if not unweighted samples, some sample number is required, so set to default.
+            params['n_samples'] = MonteCarlo.__DEFAULT_N_SAMPLES
+        elif isinstance(state, UnweightedSamples) and params['n_samples'] is None:
+            params['n_samples'] = len(state)  # number of samples is from provided state
+
         if len(params['events']) == 0 and 'horizon' not in params:
             raise ValueError("If specifying no event (i.e., simulate to time), must specify horizon")
 
-        # Sample from state if n_samples specified or state is not UnweightedSamples
-        if not isinstance(state, UnweightedSamples) or len(state) != params['n_samples']:
+        # Sample from state if n_samples specified or state is not UnweightedSamples (Case 2)
+        # Or if is Unweighted samples, but there are the wrong number of samples (Case 1)
+        if (
+            (isinstance(state, UnweightedSamples) and len(state) != params['n_samples'])  # Case 1
+            or not isinstance(state, UnweightedSamples)):  # Case 2
             state = state.sample(params['n_samples'])
 
         es_eqn = self.model.event_state

src/progpy/predictors/unscented_transform.py

@@ -123,14 +123,14 @@ def state_transition(x, dt):
         self.filter = kalman.UnscentedKalmanFilter(num_states, num_measurements, self.parameters['dt'], measure, state_transition, self.sigma_points)
         self.filter.Q = self.parameters['Q']
 
-    def predict(self, state, future_loading_eqn: Callable, **kwargs) -> PredictionResults:
+    def predict(self, state, future_loading_eqn: Callable = None, **kwargs) -> PredictionResults:
         """
         Perform a single prediction
 
         Parameters
         ----------
         state (UncertaintData): Distribution of states
-        future_loading_eqn : function (t, x={}) -> z
+        future_loading_eqn: function (t, x=None) -> z, optional
             Function to generate an estimate of loading at future time t
         options (optional, kwargs): configuration options\n
         Any additional configuration values. Note: These parameters can also be specified in the predictor constructor. The following configuration parameters are supported: \n
@@ -169,6 +169,9 @@ def predict(self, state, future_loading_eqn: Callable, **kwargs) -> PredictionRe
         else:
             raise TypeError("state must be UncertainData, dict, or StateContainer")
 
+        if future_loading_eqn is None:
+            future_loading_eqn = lambda t, x=None: self.model.InputContainer({})
+
         params = deepcopy(self.parameters) # copy parameters
         params.update(kwargs) # update for specific run
 

tests/test_predictors.py

@@ -66,10 +66,9 @@ def test_UTP_ThrownObject(self):
         m = ThrownObject()
         pred = UnscentedTransformPredictor(m)
         samples = MultivariateNormalDist(['x', 'v'], [1.83, 40], [[0.1, 0.01], [0.01, 0.1]])
-        def future_loading(t, x={}):
-            return {}
 
-        mc_results = pred.predict(samples, future_loading, dt=0.01, save_freq=1)
+        # No future loading (i.e., no load)
+        mc_results = pred.predict(samples, dt=0.01, save_freq=1)
         self.assertAlmostEqual(mc_results.time_of_event.mean['impact'], 8.21, 0)
         self.assertAlmostEqual(mc_results.time_of_event.mean['falling'], 4.15, 0)
         # self.assertAlmostEqual(mc_results.times[-1], 9, 1)  # Saving every second, last time should be around the 1s after impact event (because one of the sigma points fails afterwards)
@@ -126,10 +125,9 @@ def future_loading(t, x=None):
     def test_MC(self):
         m = ThrownObject()
         mc = MonteCarlo(m)
-        def future_loading(t=None, x=None):
-            return {}
-            
-        mc.predict(m.initialize(), future_loading, dt=0.2, num_samples=3, save_freq=1)
+        
+        # Test with empty future loading (i.e., no load)
+        mc.predict(m.initialize(), dt=0.2, num_samples=3, save_freq=1)
 
     def test_prediction_mvnormaldist(self):
         times = list(range(10))
@@ -385,6 +383,48 @@ def test_utp_surrogate(self):
 
     def test_mc_surrogate(self):
         self._test_surrogate_pred(MonteCarlo)
+    
+    def test_mc_num_samples(self):
+        """
+        This test confirms that monte carlos sampling logic works as expected
+        """
+        m = ThrownObject()
+        def future_load(t, x=None):
+            return m.InputContainer({})
+
+        pred = MonteCarlo(m)
+
+        # First test- scalar input
+        x_scalar = ScalarData({'x': 10, 'v': 0})
+        # Should default to 100 samples
+        result = pred.predict(x_scalar, future_load)
+        self.assertEqual(len(result.time_of_event), 100)
+        # Repeat with less samples
+        result = pred.predict(x_scalar, future_load, n_samples=10)
+        self.assertEqual(len(result.time_of_event), 10)
+        
+        # Second test- Same, but with multivariate normal input
+        # Behavior should be the same
+        x_mvnormal = MultivariateNormalDist(['x', 'v'], [10, 0], [[0.1, 0], [0, 0.1]])
+        # Should default to 100 samples
+        result = pred.predict(x_mvnormal, future_load)
+        self.assertEqual(len(result.time_of_event), 100)
+        # Repeat with less samples
+        result = pred.predict(x_mvnormal, future_load, n_samples=10)
+        self.assertEqual(len(result.time_of_event), 10)
+
+        # Third test- UnweightedSamples input
+        x_uwsamples = UnweightedSamples([{'x': 10, 'v': 0}, {'x': 9.9, 'v': 0.1}, {'x': 10.1, 'v': -0.1}])
+        # Should default to same as x_uwsamples - HERE IS THE DIFFERENCE FROM OTHER TYPES
+        result = pred.predict(x_uwsamples, future_load)
+        self.assertEqual(len(result.time_of_event), 3)
+        # Should be exact same data, in the same order
+        for i in range(3):
+            self.assertEqual(result.states[i][0]['x'], x_uwsamples[i]['x'])
+            self.assertEqual(result.states[i][0]['v'], x_uwsamples[i]['v'])
+        # Repeat with more samples
+        result = pred.predict(x_uwsamples, future_load, n_samples=10)
+        self.assertEqual(len(result.time_of_event), 10)
 
 # This allows the module to be executed directly    
 def main():