Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 15 additions & 8 deletions stonesoup/measures/state.py
Original file line number Diff line number Diff line change
Expand Up @@ -188,32 +188,39 @@ def __call__(self, state1, state2):

Parameters
----------
state1 : :class:`~.State`
state2 : :class:`~.State`
state1 : :class:`~.State` whose `state_vector` has shape (n, 1)
state2 : :class:`~.State` whose `state_vector` has shape (n, m)

Returns
-------
float
Union(float, np.ndarray)
Squared Mahalanobis distance between a pair of input :class:`~.State`
objects
objects. Returns a 1D array if state2 has >1 length along axis 1.

"""
state_vector1 = getattr(state1, 'mean', state1.state_vector)
state_vector2 = getattr(state2, 'mean', state2.state_vector)

if len(state_vector1) != len(state_vector2):
raise ValueError(
f"Shape mismatch between state1 and state2 along axis 0, \
{len(state_vector1)} != {len(state_vector2)}."
)

if self.mapping is not None:
u = state_vector1[self.mapping, 0]
v = state_vector2[self.mapping2, 0]
v = state_vector2[self.mapping2, :]
# extract the mapped covariance data
vi = self._inv_cov(state1, tuple(self.mapping))
else:
u = state_vector1[:, 0]
v = state_vector2[:, 0]
v = state_vector2[:, :]
vi = self._inv_cov(state1)

delta = u - v
delta = -(v.T-u)

return np.dot(np.dot(delta, vi), delta)
# Return the diagonal elements of A@B@A.T
return np.einsum('ij,jk,ik->i', delta, vi, delta).squeeze()[()]

@staticmethod
def _inv_cov(state, mapping=None):
Expand Down
55 changes: 46 additions & 9 deletions stonesoup/measures/tests/test_state.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,11 @@
stateB_u = State(u, timestamp=t)

v = StateVector([[11.], [10.], [100.], [2.]])
vm = StateVectors(np.random.random((4, 10)))
vi = CovarianceMatrix(np.diag([20., 3., 7., 10.]))

state_v = GaussianState(v, vi, timestamp=t)
state_vm = State(vm)
stateB_v = State(v, timestamp=t)


Expand All @@ -48,9 +50,14 @@ def test_euclideanweighted():

def test_mahalanobis():
measure = measures.Mahalanobis()
assert measure(state_u, state_v) == distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui))
assert measure(state_u, state_v) == pytest.approx(distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui)))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[:, 0],
vec[:, 0],
np.linalg.inv(ui)))


def test_hellinger():
Expand Down Expand Up @@ -141,13 +148,23 @@ def test_hellinger_partial_mapping(mapping_type):
def test_mahalanobis_full_mapping(mapping_type):
mapping = mapping_type(np.arange(len(u)))
measure = measures.Mahalanobis(mapping=mapping)
assert measure(state_u, state_v) == distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui))
assert measure(state_u, state_v) == pytest.approx(distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui)))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[:, 0],
vec[:, 0],
np.linalg.inv(ui)))
measure = measures.Mahalanobis(mapping=mapping, mapping2=mapping)
assert measure(state_u, state_v) == distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui))
assert measure(state_u, state_v) == pytest.approx(distance.mahalanobis(u[:, 0],
v[:, 0],
np.linalg.inv(ui)))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[:, 0],
vec[:, 0],
np.linalg.inv(ui)))


def test_mahalanobis_partial_mapping(mapping_type):
Expand All @@ -157,23 +174,43 @@ def test_mahalanobis_partial_mapping(mapping_type):
assert measure(state_u, state_v) == \
distance.mahalanobis([10, 1],
[11, 10], np.linalg.inv(reduced_ui))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[mapping, 0],
vec[mapping, 0],
np.linalg.inv(reduced_ui)))
mapping = np.array([0, 3])
reduced_ui = CovarianceMatrix(np.diag([100, 10]))
measure = measures.Mahalanobis(mapping=mapping)
assert measure(state_u, state_v) == \
distance.mahalanobis([10, 1],
[11, 2], np.linalg.inv(reduced_ui))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[mapping, 0],
vec[mapping, 0],
np.linalg.inv(reduced_ui)))

mapping = mapping_type([0, 1])
measure = measures.Mahalanobis(mapping=mapping, mapping2=mapping)
assert measure(state_u, state_v) == \
distance.mahalanobis([10, 1],
[11, 10], np.linalg.inv(reduced_ui))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[mapping, 0],
vec[mapping, 0],
np.linalg.inv(reduced_ui)))
mapping = np.array([0, 3])
measure = measures.Mahalanobis(mapping=mapping, mapping2=mapping)
assert measure(state_u, state_v) == \
distance.mahalanobis([10, 1],
[11, 2], np.linalg.inv(reduced_ui))
result_nm = measure(state_u, state_vm)
for i, vec in enumerate(vm):
assert result_nm[i] == pytest.approx(distance.mahalanobis(u[mapping, 0],
vec[mapping, 0],
np.linalg.inv(reduced_ui)))

mapping = mapping_type([0, 1])
mapping2 = np.array([0, 3])
Expand Down