From e8aa6b7decd9f90205e164df422b01c7bfef4c7b Mon Sep 17 00:00:00 2001 From: Steven Hiscocks Date: Thu, 16 Mar 2023 16:20:14 +0000 Subject: [PATCH 1/3] Add log_weights to ParticleState This may have performance advantage by simply storing weights in log space, and converting to weights dynamically when requried (utilising a cache for the property). Changes here are mostly backwards compatible, other than an issue with `from_state` method. --- stonesoup/predictor/particle.py | 10 ++- stonesoup/resampler/particle.py | 9 +-- stonesoup/types/state.py | 94 ++++++++++++++++-------- stonesoup/types/tests/test_prediction.py | 2 +- stonesoup/updater/particle.py | 31 ++++---- 5 files changed, 93 insertions(+), 53 deletions(-) diff --git a/stonesoup/predictor/particle.py b/stonesoup/predictor/particle.py index 97f7ed205..5a9e04076 100644 --- a/stonesoup/predictor/particle.py +++ b/stonesoup/predictor/particle.py @@ -48,7 +48,7 @@ def predict(self, prior, timestamp=None, **kwargs): time_interval=time_interval, **kwargs) - return Prediction.from_state(prior, state_vector=new_state_vector, weight=prior.weight, + return Prediction.from_state(prior, state_vector=new_state_vector, weight=None, timestamp=timestamp, particle_list=None, transition_model=self.transition_model) @@ -93,9 +93,11 @@ def predict(self, prior, *args, **kwargs): *args, **kwargs) return Prediction.from_state(prior, state_vector=particle_prediction.state_vector, - weight=particle_prediction.weight, + weight=None, + log_weight=particle_prediction.log_weight, timestamp=particle_prediction.timestamp, - fixed_covar=kalman_prediction.covar, particle_list=None, + fixed_covar=kalman_prediction.covar, + particle_list=None, transition_model=self.transition_model) @@ -145,6 +147,7 @@ def predict(self, prior, timestamp=None, **kwargs): state_vector=copy.copy(prior.state_vector), parent=prior, particle_list=None, + weight=None, dynamic_model=copy.copy(prior.dynamic_model), timestamp=timestamp) @@ -213,6 +216,7 @@ def predict(self, prior, timestamp=None, **kwargs): new_particle_state = Prediction.from_state( prior, state_vector=copy.copy(prior.state_vector), + weight=None, parent=prior, particle_list=None, timestamp=timestamp) diff --git a/stonesoup/resampler/particle.py b/stonesoup/resampler/particle.py index ca242c32c..67c7696c4 100644 --- a/stonesoup/resampler/particle.py +++ b/stonesoup/resampler/particle.py @@ -2,7 +2,6 @@ from .base import Resampler from ..base import Property -from ..types.numeric import Probability from ..types.state import ParticleState @@ -25,9 +24,8 @@ def resample(self, particles): if not isinstance(particles, ParticleState): particles = ParticleState(None, particle_list=particles) n_particles = len(particles) - weight = Probability(1 / n_particles) - log_weights = np.asfarray(np.log(particles.weight)) + log_weights = particles.log_weight weight_order = np.argsort(log_weights, kind='stable') max_log_value = log_weights[weight_order[-1]] with np.errstate(divide='ignore'): @@ -43,7 +41,7 @@ def resample(self, particles): index = weight_order[np.searchsorted(cdf, np.log(u_j))] new_particles = particles[index] - new_particles.weight = np.full((n_particles, ), weight) + new_particles.log_weight = np.full((n_particles, ), np.log(1/n_particles)) return new_particles @@ -84,7 +82,8 @@ def resample(self, particles): particles = ParticleState(None, particle_list=particles) if self.threshold is None: self.threshold = len(particles) / 2 - if 1 / np.sum(np.square(particles.weight)) < self.threshold: # If ESS too small, resample + # If ESS too small, resample + if 1 / np.sum(np.exp(2*particles.log_weight)) < self.threshold: return self.resampler.resample(self.resampler, particles) else: return particles diff --git a/stonesoup/types/state.py b/stonesoup/types/state.py index 7f79b6cb6..3574edee9 100644 --- a/stonesoup/types/state.py +++ b/stonesoup/types/state.py @@ -73,7 +73,7 @@ def from_state(state: 'State', *args: Any, target_type: Optional[typing.Type] = new_kwargs = { name: getattr(state, name) for name in type(state).properties.keys() & target_type.properties.keys() - if name not in args_property_names} + if name not in args_property_names and name not in kwargs} new_kwargs.update(kwargs) @@ -617,6 +617,7 @@ class ParticleState(State): state_vector: StateVectors = Property(doc='State vectors.') weight: MutableSequence[Probability] = Property(default=None, doc='Weights of particles') + log_weight: np.ndarray = Property(default=None, doc='Log weights of particles') parent: 'ParticleState' = Property(default=None, doc='Parent particles') particle_list: MutableSequence[Particle] = Property(default=None, doc='List of Particle objects') @@ -625,6 +626,22 @@ class ParticleState(State): 'weighted sample covariance is then used.') def __init__(self, *args, **kwargs): + weight = next( + (val for name, val in zip(type(self).properties, args) if name == 'weight'), + kwargs.get('weight', None)) + log_weight, idx = next( + ((val, idx) for idx, (name, val) in enumerate(zip(type(self).properties, args)) + if name == 'log_weight'), + (kwargs.get('log_weight', None), None)) + + if weight is not None and log_weight is not None: + raise ValueError("Cannot provide both weight and log weight") + elif log_weight is None and weight is not None: + log_weight = np.log(np.asfarray(weight)) + if idx is not None: + args[idx] = log_weight + else: + kwargs['log_weight'] = log_weight super().__init__(*args, **kwargs) if (self.particle_list is not None) and \ @@ -650,8 +667,6 @@ def __init__(self, *args, **kwargs): if self.state_vector is not None and not isinstance(self.state_vector, StateVectors): self.state_vector = StateVectors(self.state_vector) - if self.weight is not None and not isinstance(self.weight, np.ndarray): - self.weight = np.array(self.weight) def __getitem__(self, item): if self.parent is not None: @@ -659,25 +674,26 @@ def __getitem__(self, item): else: parent = None - if self.weight is not None: - weight = self.weight[item] + if self.log_weight is not None: + log_weight = self.log_weight[item] else: - weight = None + log_weight = None if isinstance(item, int): result = Particle(state_vector=self.state_vector[:, item], - weight=weight, + weight=self.weight[item] if self.weight is not None else None, parent=parent) else: # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=weight, + weight=None, + log_weight=log_weight, parent=parent, particle_list=None) return result - @clearable_cached_property('state_vector', 'weight') + @clearable_cached_property('state_vector', 'log_weight') def particles(self): """Sequence of individual :class:`~.Particle` objects.""" if self.particle_list is not None: @@ -692,22 +708,38 @@ def ndim(self): """The number of dimensions represented by the state.""" return self.state_vector.shape[0] - @clearable_cached_property('state_vector', 'weight') + @clearable_cached_property('state_vector', 'log_weight') def mean(self): """Sample mean for particles""" if len(self) == 1: # No need to calculate mean return self.state_vector - return np.average(self.state_vector, axis=1, weights=np.asfarray(self.weight)) + return np.average(self.state_vector, axis=1, weights=np.exp(self.log_weight)) - @clearable_cached_property('state_vector', 'weight', 'fixed_covar') + @clearable_cached_property('state_vector', 'log_weight', 'fixed_covar') def covar(self): """Sample covariance matrix for particles""" if self.fixed_covar is not None: return self.fixed_covar - return np.cov(self.state_vector, ddof=0, aweights=np.asfarray(self.weight)) + return np.cov(self.state_vector, ddof=0, aweights=np.exp(self.log_weight)) + + @weight.setter + def weight(self, value): + self.log_weight = np.log(np.asfarray(value)) + @weight.getter + def weight(self): + try: + return self.__dict__['weight'] + except KeyError: + log_weight = self.log_weight + if log_weight is None: + return None + weight = Probability.from_log_ufunc(log_weight) + self.__dict__['weight'] = weight + return weight State.register(ParticleState) # noqa: E305 +ParticleState.log_weight._clear_cached.add('weight') class MultiModelParticleState(ParticleState): @@ -734,10 +766,10 @@ def __getitem__(self, item): else: parent = None - if self.weight is not None: - weight = self.weight[item] + if self.log_weight is not None: + log_weight = self.log_weight[item] else: - weight = None + log_weight = None if self.dynamic_model is not None: dynamic_model = self.dynamic_model[item] @@ -745,15 +777,17 @@ def __getitem__(self, item): dynamic_model = None if isinstance(item, int): - result = MultiModelParticle(state_vector=self.state_vector[:, item], - weight=weight, - parent=parent, - dynamic_model=dynamic_model) + result = MultiModelParticle( + state_vector=self.state_vector[:, item], + weight=self.weight[item] if self.weight is not None else None, + parent=parent, + dynamic_model=dynamic_model) else: # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=weight, + weight=None, + log_weight=log_weight, parent=parent, particle_list=None, dynamic_model=dynamic_model) @@ -780,10 +814,10 @@ def __getitem__(self, item): else: parent = None - if self.weight is not None: - weight = self.weight[item] + if self.log_weight is not None: + log_weight = self.log_weight[item] else: - weight = None + log_weight = None if self.model_probabilities is not None: model_probabilities = self.model_probabilities[:, item] @@ -791,15 +825,17 @@ def __getitem__(self, item): model_probabilities = None if isinstance(item, int): - result = RaoBlackwellisedParticle(state_vector=self.state_vector[:, item], - weight=weight, - parent=parent, - model_probabilities=model_probabilities) + result = RaoBlackwellisedParticle( + state_vector=self.state_vector[:, item], + weight=self.weight[item] if self.weight is not None else None, + parent=parent, + model_probabilities=model_probabilities) else: # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=weight, + weight=None, + log_weight=log_weight, parent=parent, particle_list=None, model_probabilities=model_probabilities) diff --git a/stonesoup/types/tests/test_prediction.py b/stonesoup/types/tests/test_prediction.py index e1b8a7f68..0370f5444 100644 --- a/stonesoup/types/tests/test_prediction.py +++ b/stonesoup/types/tests/test_prediction.py @@ -167,7 +167,7 @@ def test_from_state(prediction_type): assert prediction.tag == state.tag state = ParticleState(StateVectors([[1]]), weight=0.5, timestamp=datetime.datetime.now()) - prediction = prediction_type.from_state(state) + prediction = prediction_type.from_state(state, weight=None) if prediction_type is Prediction: assert isinstance(prediction, ParticleStatePrediction) else: diff --git a/stonesoup/updater/particle.py b/stonesoup/updater/particle.py index 4dc2f3855..7ca876138 100644 --- a/stonesoup/updater/particle.py +++ b/stonesoup/updater/particle.py @@ -11,7 +11,6 @@ from ..functions import cholesky_eps, sde_euler_maruyama_integration from ..predictor.particle import MultiModelPredictor, RaoBlackwellisedMultiModelPredictor from ..resampler import Resampler -from ..types.numeric import Probability from ..types.prediction import ( Prediction, ParticleMeasurementPrediction, GaussianStatePrediction, MeasurementPrediction) from ..types.update import ParticleStateUpdate, Update @@ -51,13 +50,13 @@ def update(self, hypothesis, **kwargs): else: measurement_model = hypothesis.measurement.measurement_model - new_weight = np.asfarray(np.log(predicted_state.weight)) + measurement_model.logpdf( + new_weight = predicted_state.log_weight + measurement_model.logpdf( hypothesis.measurement, predicted_state, **kwargs) # Normalise the weights new_weight -= logsumexp(new_weight) - predicted_state.weight = Probability.from_log_ufunc(new_weight) + predicted_state.log_weight = new_weight # Resample if self.resampler is not None: @@ -66,7 +65,8 @@ def update(self, hypothesis, **kwargs): return Update.from_state( state=hypothesis.prediction, state_vector=predicted_state.state_vector, - weight=predicted_state.weight, + weight=None, + log_weight=predicted_state.log_weight, hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, particle_list=None @@ -80,11 +80,10 @@ def predict_measurement(self, state_prediction, measurement_model=None, measurement_model = self.measurement_model new_state_vector = measurement_model.function(state_prediction, **kwargs) - new_weight = state_prediction.weight return MeasurementPrediction.from_state( state_prediction, state_vector=new_state_vector, timestamp=state_prediction.timestamp, - particle_list=None, weight=new_weight) + particle_list=None, weight=None) class GromovFlowParticleUpdater(Updater): @@ -253,6 +252,7 @@ def update(self, hypothesis, **kwargs): update = Update.from_state( hypothesis.prediction, + weight=None, hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, particle_list=None @@ -260,12 +260,12 @@ def update(self, hypothesis, **kwargs): transition_matrix = np.asanyarray(self.predictor.transition_matrix) - update.weight = update.weight \ - * measurement_model.pdf(hypothesis.measurement, update, **kwargs) \ - * transition_matrix[update.parent.dynamic_model, update.dynamic_model] + update.log_weight = update.log_weight \ + + measurement_model.logpdf(hypothesis.measurement, update, **kwargs) \ + + np.log(transition_matrix[update.parent.dynamic_model, update.dynamic_model]) # Normalise the weights - update.weight /= Probability.sum(update.weight) + update.log_weight -= logsumexp(update.log_weight) if self.resampler: update = self.resampler.resample(update) @@ -300,7 +300,8 @@ def update(self, hypothesis, **kwargs): update = Update.from_state( hypothesis.prediction, - weight=copy.copy(hypothesis.prediction.weight), + weight=None, + log_weight=copy.copy(hypothesis.prediction.log_weight), hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, particle_list=None @@ -309,12 +310,12 @@ def update(self, hypothesis, **kwargs): update.model_probabilities = self.calculate_model_probabilities( hypothesis.prediction, self.predictor) - update.weight = update.weight * measurement_model.pdf(hypothesis.measurement, - update, - **kwargs) + update.log_weight = update.log_weight + measurement_model.logpdf(hypothesis.measurement, + update, + **kwargs) # Normalise the weights - update.weight /= Probability.sum(update.weight) + update.log_weight -= logsumexp(update.log_weight) if self.resampler: update = self.resampler.resample(update) From 5ace9d533f9b94bd07299e60841c3a77714fe88a Mon Sep 17 00:00:00 2001 From: Steven Hiscocks Date: Fri, 17 Mar 2023 09:35:55 +0000 Subject: [PATCH 2/3] Customise from_state method on ParticleState This handles the fact that the weight and particle list should be ignored for created states, as the elements of these are converted/stored in other parts of an instance, and only used during instantiation of the class. --- stonesoup/predictor/particle.py | 14 +++----- stonesoup/types/state.py | 42 ++++++++++++++---------- stonesoup/types/tests/test_prediction.py | 2 +- stonesoup/updater/particle.py | 11 +------ 4 files changed, 32 insertions(+), 37 deletions(-) diff --git a/stonesoup/predictor/particle.py b/stonesoup/predictor/particle.py index 5a9e04076..57dafad25 100644 --- a/stonesoup/predictor/particle.py +++ b/stonesoup/predictor/particle.py @@ -48,8 +48,9 @@ def predict(self, prior, timestamp=None, **kwargs): time_interval=time_interval, **kwargs) - return Prediction.from_state(prior, state_vector=new_state_vector, weight=None, - timestamp=timestamp, particle_list=None, + return Prediction.from_state(prior, + state_vector=new_state_vector, + timestamp=timestamp, transition_model=self.transition_model) @@ -92,12 +93,11 @@ def predict(self, prior, *args, **kwargs): GaussianState(prior.state_vector, prior.covar, prior.timestamp), *args, **kwargs) - return Prediction.from_state(prior, state_vector=particle_prediction.state_vector, - weight=None, + return Prediction.from_state(prior, + state_vector=particle_prediction.state_vector, log_weight=particle_prediction.log_weight, timestamp=particle_prediction.timestamp, fixed_covar=kalman_prediction.covar, - particle_list=None, transition_model=self.transition_model) @@ -146,8 +146,6 @@ def predict(self, prior, timestamp=None, **kwargs): prior, state_vector=copy.copy(prior.state_vector), parent=prior, - particle_list=None, - weight=None, dynamic_model=copy.copy(prior.dynamic_model), timestamp=timestamp) @@ -216,9 +214,7 @@ def predict(self, prior, timestamp=None, **kwargs): new_particle_state = Prediction.from_state( prior, state_vector=copy.copy(prior.state_vector), - weight=None, parent=prior, - particle_list=None, timestamp=timestamp) # Change the value of the dynamic value randomly according to the diff --git a/stonesoup/types/state.py b/stonesoup/types/state.py index 3574edee9..bd1d42a7f 100644 --- a/stonesoup/types/state.py +++ b/stonesoup/types/state.py @@ -163,7 +163,7 @@ def from_state( if target_type is None: target_type = CreatableFromState.class_mapping[cls][state_type] - return State.from_state(state, *args, **kwargs, target_type=target_type) + return target_type.from_state(state, *args, **kwargs, target_type=target_type) class ASDState(Type): @@ -687,12 +687,32 @@ def __getitem__(self, item): # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=None, log_weight=log_weight, - parent=parent, - particle_list=None) + parent=parent) return result + @classmethod + def from_state(cls, state: 'State', *args: Any, target_type: Optional[typing.Type] = None, + **kwargs: Any) -> 'State': + + # Handle default presence of both particle_list and weight once class has been created by + # ignoring particle_list and weight (setting to None) if not provided. + particle_list, particle_list_idx = next( + ((val, idx) for idx, (name, val) in enumerate(zip(cls.properties, args)) + if name == 'particle_list'), + (kwargs.get('particle_list', None), None)) + if particle_list_idx is None: + kwargs['particle_list'] = particle_list + + weight, weight_idx = next( + ((val, idx) for idx, (name, val) in enumerate(zip(cls.properties, args)) + if name == 'weight'), + (kwargs.get('weight', None), None)) + if weight_idx is None: + kwargs['weight'] = weight + + return super().from_state(state, *args, target_type=target_type, **kwargs) + @clearable_cached_property('state_vector', 'log_weight') def particles(self): """Sequence of individual :class:`~.Particle` objects.""" @@ -786,10 +806,8 @@ def __getitem__(self, item): # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=None, log_weight=log_weight, parent=parent, - particle_list=None, dynamic_model=dynamic_model) return result @@ -834,15 +852,13 @@ def __getitem__(self, item): # Allow for Prediction/Update sub-types result = type(self).from_state(self, state_vector=self.state_vector[:, item], - weight=None, log_weight=log_weight, parent=parent, - particle_list=None, model_probabilities=model_probabilities) return result -class EnsembleState(Type): +class EnsembleState(State): r"""Ensemble State type This is an Ensemble state object which describes the system state as a @@ -925,11 +941,6 @@ def generate_ensemble(mean, covar, num_vectors): return ensemble - @property - def ndim(self): - """Number of dimensions in state vectors""" - return np.shape(self.state_vector)[0] - @property def num_vectors(self): """Number of columns in state ensemble""" @@ -953,9 +964,6 @@ def sqrt_covar(self): / np.sqrt(self.num_vectors - 1)) -State.register(EnsembleState) # noqa: E305 - - class CategoricalState(State): r"""CategoricalState type. diff --git a/stonesoup/types/tests/test_prediction.py b/stonesoup/types/tests/test_prediction.py index 0370f5444..e1b8a7f68 100644 --- a/stonesoup/types/tests/test_prediction.py +++ b/stonesoup/types/tests/test_prediction.py @@ -167,7 +167,7 @@ def test_from_state(prediction_type): assert prediction.tag == state.tag state = ParticleState(StateVectors([[1]]), weight=0.5, timestamp=datetime.datetime.now()) - prediction = prediction_type.from_state(state, weight=None) + prediction = prediction_type.from_state(state) if prediction_type is Prediction: assert isinstance(prediction, ParticleStatePrediction) else: diff --git a/stonesoup/updater/particle.py b/stonesoup/updater/particle.py index 7ca876138..c01b8337f 100644 --- a/stonesoup/updater/particle.py +++ b/stonesoup/updater/particle.py @@ -65,11 +65,9 @@ def update(self, hypothesis, **kwargs): return Update.from_state( state=hypothesis.prediction, state_vector=predicted_state.state_vector, - weight=None, log_weight=predicted_state.log_weight, hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, - particle_list=None ) @lru_cache() @@ -82,8 +80,7 @@ def predict_measurement(self, state_prediction, measurement_model=None, new_state_vector = measurement_model.function(state_prediction, **kwargs) return MeasurementPrediction.from_state( - state_prediction, state_vector=new_state_vector, timestamp=state_prediction.timestamp, - particle_list=None, weight=None) + state_prediction, state_vector=new_state_vector, timestamp=state_prediction.timestamp) class GromovFlowParticleUpdater(Updater): @@ -204,7 +201,6 @@ def update(self, hypothesis, **kwargs): hypothesis, weight=particle_update.weight, fixed_covar=kalman_update.covar, - particle_list=None, timestamp=particle_update.timestamp) def predict_measurement(self, state_prediction, *args, **kwargs): @@ -221,7 +217,6 @@ def predict_measurement(self, state_prediction, *args, **kwargs): state_vector=particle_prediction.state_vector, weight=state_prediction.weight, fixed_covar=kalman_prediction.covar, - particle_list=None, timestamp=particle_prediction.timestamp) @@ -252,10 +247,8 @@ def update(self, hypothesis, **kwargs): update = Update.from_state( hypothesis.prediction, - weight=None, hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, - particle_list=None ) transition_matrix = np.asanyarray(self.predictor.transition_matrix) @@ -300,11 +293,9 @@ def update(self, hypothesis, **kwargs): update = Update.from_state( hypothesis.prediction, - weight=None, log_weight=copy.copy(hypothesis.prediction.log_weight), hypothesis=hypothesis, timestamp=hypothesis.measurement.timestamp, - particle_list=None ) update.model_probabilities = self.calculate_model_probabilities( From b2c3da12582ce804744f770cb7dd7594be4e8af6 Mon Sep 17 00:00:00 2001 From: Steven Hiscocks Date: Tue, 11 Apr 2023 08:57:48 +0100 Subject: [PATCH 3/3] Cache weight on ParticleState when setting weight --- stonesoup/types/state.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/stonesoup/types/state.py b/stonesoup/types/state.py index bd1d42a7f..b6fc0ba07 100644 --- a/stonesoup/types/state.py +++ b/stonesoup/types/state.py @@ -744,7 +744,11 @@ def covar(self): @weight.setter def weight(self, value): - self.log_weight = np.log(np.asfarray(value)) + if value is None: + self.log_weight = None + else: + self.log_weight = np.log(np.asfarray(value)) + self.__dict__['weight'] = np.asanyarray(value) @weight.getter def weight(self):