Repository: matchms/matchms
Pull request #40 (2020-05-20 15:17:36)

59 cosine hungarian
Originally here https://github.com/matchms/matchms-backup/pull/246

So far we implemented a greedy version of the cosine score (which is often used because it is faster and gives the same results in 99,9% (that's a guess) of the actual cases.

The Hungarian algorithm would be the fully correct solution of the optimization problem, so it is nice to have it as well.

Here an example where the greedy scores would fail:
```python
from matchms import Spectrum
from matchms.similarity import CosineGreedy, CosineGreedyNumba, CosineHungarian

test_spectrum1 = Spectrum(mz=np.array([100.005, 100.016]),
                         intensities=np.array([1.0, 0.9]),
                         metadata={})

test_spectrum2 = Spectrum(mz=np.array([100.005, 100.01]),
                         intensities=np.array([0.9, 1.0]),
                         metadata={})

similarity_measure = CosineGreedyNumba(tolerance=0.01)
print("CosineGreedy:", similarity_measure(test_spectrum1, test_spectrum2))

similarity_measure = CosineGreedy(tolerance=0.01)
print("CosineGreedyNumba:", similarity_measure(test_spectrum1, test_spectrum2))

similarity_measure = CosineHungarian(tolerance=0.01)
print("CosineHungarian:", similarity_measure(test_spectrum1, test_spectrum2))
```
``CosineGreedy: (0.5524861878453039, 1)``
``CosineGreedyNumba: (0.5524861878453039, 1)``
``CosineHungarian: (0.994475138121547, 2)``

CunliangGeng: I had a look at the test scripts and would skip these tests to solve this issue currently.
We need to create more issues to 
1. check, clean and format existing test cases before adding new tests in #4, which first needs a systematic design of tests.
2. format all source code based on pep8   
CunliangGeng: PR #43
fdiblen: @CunliangGeng thanks! Can you create an issue about linting?
CunliangGeng: What is the #32 about? 
fdiblen: Sorry, missed it. Could you update the issue to add some extra details (e.g. doing it in tests)?
CunliangGeng: Updated #32 and #4. The issue can be closed now.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5af93063..04281710 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,7 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
-
+- Cosine score based on Hungarian algorithm [#40](https://github.com/matchms/matchms/pull/40)
 - Modified cosine score [#26](https://github.com/matchms/matchms/issues/26)
 - Doc strings for many methods [#49](https://github.com/matchms/matchms/issues/49)
 - Examples in doc strings which are tested on CI
diff --git a/matchms/similarity/CosineHungarian.py b/matchms/similarity/CosineHungarian.py
new file mode 100644
index 00000000..f1e3f626
--- /dev/null
+++ b/matchms/similarity/CosineHungarian.py
@@ -0,0 +1,96 @@
+from typing import Tuple
+import numpy
+from scipy.optimize import linear_sum_assignment
+from matchms.similarity.collect_peak_pairs import collect_peak_pairs
+from matchms.typing import SpectrumType
+
+
+class CosineHungarian:
+    """Calculate 'cosine similarity score' between two spectra (using Hungarian algorithm).
+
+    The cosine score aims at quantifying the similarity between two mass spectra.
+    The score is calculated by finding best possible matches between peaks
+    of two spectra. Two peaks are considered a potential match if their
+    m/z ratios lie within the given 'tolerance'.
+    The underlying peak assignment problem is here solved using the Hungarian algorithm.
+    This can perform notably slower than the 'greedy' implementation in CosineGreedy, but
+    does represent a mathematically proper solution to the problem.
+
+    """
+    def __init__(self, tolerance=0.1):
+        """
+
+        Args:
+        ----
+        tolerance: float
+            Peaks will be considered a match when <= tolerance appart. Default is 0.1.
+        """
+        self.tolerance = tolerance
+
+    def __call__(self, spectrum1: SpectrumType, spectrum2: SpectrumType) -> Tuple[float, int]:
+        """Calculate cosine score between two spectra.
+
+        Args:
+        ----
+        spectrum1: SpectrumType
+            Input spectrum 1.
+        spectrum2: SpectrumType
+            Input spectrum 2.
+        """
+        def get_peaks_arrays():
+            """Get peaks mz and intensities as numpy array."""
+            spec1 = numpy.vstack((spectrum1.peaks.mz, spectrum1.peaks.intensities)).T
+            spec2 = numpy.vstack((spectrum2.peaks.mz, spectrum2.peaks.intensities)).T
+            assert max(spec1[:, 1]) <= 1, ("Input spectrum1 is not normalized. ",
+                                           "Apply 'normalize_intensities' filter first.")
+            assert max(spec2[:, 1]) <= 1, ("Input spectrum2 is not normalized. ",
+                                           "Apply 'normalize_intensities' filter first.")
+            return spec1, spec2
+
+        def get_matching_pairs():
+            """Get pairs of peaks that match within the given tolerance."""
+            matching_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=0.0)
+            return sorted(matching_pairs, key=lambda x: x[2], reverse=True)
+
+        def get_matching_pairs_matrix():
+            """Create matrix of multiplied intensities of all matching pairs
+            between spectrum1 and spectrum2.
+            Returns
+            paired_peaks1:
+                list of paired peaks of spectrum1
+            paired_peaks2:
+                list of paired peaks of spectrum2
+            matching_pairs_matrix:
+                Array of multiplied intensities between all matching peaks.
+            """
+            if len(matching_pairs) == 0:
+                return None, None, None
+            paired_peaks1 = list({x[0] for x in matching_pairs})
+            paired_peaks2 = list({x[1] for x in matching_pairs})
+            matrix_size = (len(paired_peaks1), len(paired_peaks2))
+            matching_pairs_matrix = numpy.ones(matrix_size)
+            for match in matching_pairs:
+                matching_pairs_matrix[paired_peaks1.index(match[0]),
+                                      paired_peaks2.index(match[1])] = 1 - match[2]
+            return paired_peaks1, paired_peaks2, matching_pairs_matrix
+
+        def solve_hungarian():
+            """Use hungarian agorithm to solve the linear sum assignment problem."""
+            row_ind, col_ind = linear_sum_assignment(matching_pairs_matrix)
+            score = len(row_ind) - matching_pairs_matrix[row_ind, col_ind].sum()
+            used_matches = [(paired_peaks1[x], paired_peaks2[y]) for (x, y) in zip(row_ind, col_ind)]
+            return score, used_matches
+
+        def calc_score():
+            """Calculate cosine similarity score."""
+            if matching_pairs_matrix is not None:
+                score, used_matches = solve_hungarian()
+                # Normalize score:
+                score = score/max(numpy.sum(spec1[:, 1]**2), numpy.sum(spec2[:, 1]**2))
+                return score, len(used_matches)
+            return 0.0, 0
+
+        spec1, spec2 = get_peaks_arrays()
+        matching_pairs = get_matching_pairs()
+        paired_peaks1, paired_peaks2, matching_pairs_matrix = get_matching_pairs_matrix()
+        return calc_score()
diff --git a/matchms/similarity/__init__.py b/matchms/similarity/__init__.py
index a0916968..33275a66 100644
--- a/matchms/similarity/__init__.py
+++ b/matchms/similarity/__init__.py
@@ -1,11 +1,13 @@
 """similarity module"""
 from .CosineGreedy import CosineGreedy
+from .CosineHungarian import CosineHungarian
 from .IntersectMz import IntersectMz
 from .ModifiedCosine import ModifiedCosine
 
 
 __all__ = [
     "CosineGreedy",
+    "CosineHungarian",
     "IntersectMz",
     "ModifiedCosine",
 ]
diff --git a/matchms/similarity/collect_peak_pairs.py b/matchms/similarity/collect_peak_pairs.py
index 9ba68a18..a425ab0b 100644
--- a/matchms/similarity/collect_peak_pairs.py
+++ b/matchms/similarity/collect_peak_pairs.py
@@ -13,7 +13,7 @@ def collect_peak_pairs(spec1, spec2, tolerance, shift=0):
     spec2: numpy array
         Spectrum peaks and intensities as numpy array.
     tolerance : float
-        Peaks will be considered a match when <= tolerance a part.
+        Peaks will be considered a match when <= tolerance appart.
     shift : float, optional
         Shift spectra peaks by shift. The default is 0.
 


Test patch:
diff --git a/tests/test_cosine_hungarian.py b/tests/test_cosine_hungarian.py
new file mode 100644
index 00000000..569df705
--- /dev/null
+++ b/tests/test_cosine_hungarian.py
@@ -0,0 +1,114 @@
+import numpy
+import pytest
+from matchms import Spectrum
+from matchms.similarity import CosineHungarian
+
+
+def test_cosine_hungarian_without_parameters():
+    """Compare output cosine score with own calculation on simple dummy spectrums."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 200, 300, 500, 510], dtype="float"),
+                          intensities=numpy.array([0.1, 0.2, 1.0, 0.3, 0.4], dtype="float"))
+
+    spectrum_2 = Spectrum(mz=numpy.array([100, 200, 290, 490, 510], dtype="float"),
+                          intensities=numpy.array([0.1, 0.2, 1.0, 0.3, 0.4], dtype="float"))
+    cosine_hungarian = CosineHungarian()
+    score, n_matches = cosine_hungarian(spectrum_1, spectrum_2)
+
+    # Derive expected cosine score
+    expected_matches = [0, 1, 4]  # Those peaks have matching mz values (within given tolerance)
+    multiply_matching_intensities = spectrum_1.peaks.intensities[expected_matches] \
+        * spectrum_2.peaks.intensities[expected_matches]
+    expected_score = multiply_matching_intensities.sum() / (spectrum_1.peaks.intensities ** 2).sum()
+
+    assert score == pytest.approx(expected_score, 0.0001), "Expected different cosine score."
+    assert n_matches == len(expected_matches), "Expected different number of matching peaks."
+
+
+def test_cosine_hungarian_with_tolerance_0_2():
+    """Compare output cosine score for tolerance 0.2 with own calculation on simple dummy spectrums."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 299, 300, 301, 510], dtype="float"),
+                          intensities=numpy.array([0.1, 1.0, 0.2, 0.3, 0.4], dtype="float"))
+
+    spectrum_2 = Spectrum(mz=numpy.array([100, 300, 301, 511], dtype="float"),
+                          intensities=numpy.array([0.1, 1.0, 0.3, 0.4], dtype="float"))
+    cosine_hungarian = CosineHungarian(tolerance=0.2)
+    score, n_matches = cosine_hungarian(spectrum_1, spectrum_2)
+
+    # Derive expected cosine score
+    expected_matches = [[0, 2, 3], [0, 1, 2]]  # Those peaks have matching mz values (within given tolerance)
+    multiply_matching_intensities = spectrum_1.peaks.intensities[expected_matches[0]] \
+        * spectrum_2.peaks.intensities[expected_matches[1]]
+    expected_score = multiply_matching_intensities.sum() / (spectrum_1.peaks.intensities ** 2).sum()
+
+    assert score == pytest.approx(expected_score, 0.0001), "Expected different cosine score."
+    assert n_matches == len(expected_matches[0]), "Expected different number of matching peaks."
+
+
+def test_cosine_hungarian_with_tolerance_2_0():
+    """Compare output cosine score for tolerance 2.0 with own calculation on simple dummy spectrums."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 299, 300, 301, 510], dtype="float"),
+                          intensities=numpy.array([0.1, 1.0, 0.2, 0.3, 0.4], dtype="float"))
+
+    spectrum_2 = Spectrum(mz=numpy.array([100, 300, 301, 511], dtype="float"),
+                          intensities=numpy.array([0.1, 1.0, 0.3, 0.4], dtype="float"))
+    cosine_hungarian = CosineHungarian(tolerance=2.0)
+    score, n_matches = cosine_hungarian(spectrum_1, spectrum_2)
+
+    # Derive expected cosine score
+    expected_matches = [[0, 1, 3, 4], [0, 1, 2, 3]]  # Those peaks have matching mz values (within given tolerance)
+    multiply_matching_intensities = spectrum_1.peaks.intensities[expected_matches[0]] \
+        * spectrum_2.peaks.intensities[expected_matches[1]]
+    expected_score = multiply_matching_intensities.sum() / (spectrum_1.peaks.intensities ** 2).sum()
+
+    assert score == pytest.approx(expected_score, 0.0001), "Expected different cosine score."
+    assert n_matches == len(expected_matches[0]), "Expected different number of matching peaks."
+
+
+def test_cosine_hungarian_order_of_arguments():
+    """Compare cosine scores for A,B versus B,A, which should give the same score."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 200, 299, 300, 301, 500, 510], dtype="float"),
+                          intensities=numpy.array([0.02, 0.02, 1.0, 0.2, 0.4, 0.04, 0.2], dtype="float"),
+                          metadata=dict())
+
+    spectrum_2 = Spectrum(mz=numpy.array([100, 200, 300, 301, 500, 512], dtype="float"),
+                          intensities=numpy.array([0.02, 0.02, 1.0, 0.2, 0.04, 0.2], dtype="float"),
+                          metadata=dict())
+
+    cosine_hungarian = CosineHungarian(tolerance=2.0)
+    score_1_2, n_matches_1_2 = cosine_hungarian(spectrum_1, spectrum_2)
+    score_2_1, n_matches_2_1 = cosine_hungarian(spectrum_2, spectrum_1)
+
+    assert score_1_2 == score_2_1, "Expected that the order of the arguments would not matter."
+    assert n_matches_1_2 == n_matches_2_1, "Expected that the order of the arguments would not matter."
+
+
+def test_cosine_hungarian_case_where_greedy_would_fail():
+    """Test case that would fail for cosine greedy implementations."""
+    spectrum_1 = Spectrum(mz=numpy.array([100.005, 100.016], dtype="float"),
+                          intensities=numpy.array([1.0, 0.9], dtype="float"),
+                          metadata={})
+
+    spectrum_2 = Spectrum(mz=numpy.array([100.005, 100.01], dtype="float"),
+                          intensities=numpy.array([0.9, 1.0], dtype="float"),
+                          metadata={})
+
+    cosine_hungarian = CosineHungarian(tolerance=0.01)
+    score, n_matches = cosine_hungarian(spectrum_1, spectrum_2)
+    assert score == pytest.approx(0.994475, 0.0001), "Expected different cosine score."
+    assert n_matches == 2, "Expected different number of matching peaks."
+
+
+def test_cosine_hungarian_case_without_matches():
+    """Test case for spectrums without any matching peaks."""
+    spectrum_1 = Spectrum(mz=numpy.array([100, 200], dtype="float"),
+                          intensities=numpy.array([1.0, 0.1], dtype="float"),
+                          metadata={})
+
+    spectrum_2 = Spectrum(mz=numpy.array([110, 210], dtype="float"),
+                          intensities=numpy.array([1.0, 0.1], dtype="float"),
+                          metadata={})
+
+    cosine_hungarian = CosineHungarian()
+    score, n_matches = cosine_hungarian(spectrum_1, spectrum_2)
+    assert score == 0.0, "Expected different cosine score."
+    assert n_matches == 0, "Expected different number of matching peaks."