Repository: matchms/matchms
Pull request #38 (2020-05-20 15:16:06)

238 cosine greedy numba
Originally here https://github.com/matchms/matchms-backup/pull/239


- Added ``CosineGreedyNumba`` to similarities
- Added unit tests.
- Added numba to environment.


jspaaks: Personally, I think this would be helpful. If we follow what' s described in #7 and #41 , a user with mass spec data can just do a ``conda install matchms`` (or later on, a ``pip install matchms``). They can then proceed with 

```python
from matchms.harmonization import harmfun1, harmfun2 
from matchms.similarity import simfun1, simfun2, simfun3

# get spectrum Spectrum from somewhere
# get reference_spectrums Spectrum[] from somewhere

similarity = Similarity(spectrum=spectrum, reference_spectrums=reference_spectrums)
similarity.add_harmonization_function(harmfun1, harmfun2)
similarity.add_similarity_function(simfun1, simfun2, simfun3) 
similarity.calculate()
```
In the future, when somebody wants to do this same analysis but for gene data, they can write their own ``harmfun1`` and ``harmfun2`` and import from that module like so:

```python
from otherpackage import harmfun1, harmfun2 
from matchms.similarity import simfun1, simfun2, simfun3

# get gene_sequence Spectrum from somewhere
# get reference_gene_sequences Spectrum[] from somewhere

similarity = Similarity(spectrum=gene_sequence, reference_spectrums=gene_sequences)
similarity.add_harmonization_function(harmfun1, harmfun2)
similarity.add_similarity_function(simfun1, simfun2, simfun3) 
similarity.calculate()
```

Additionally, a researcher who has an idea for a new-and-improved similarity measure, can do 

```python
from matchms.harmonization import harmfun1, harmfun2 
from matchms.similarity import simfun1, simfun2, simfun3
from yetanotherpackage import simfun4

# get spectrum Spectrum from somewhere
# get reference_spectrums Spectrum[] from somewhere

similarity = Similarity(spectrum=spectrum, reference_spectrums=reference_spectrums)
similarity.add_harmonization_function(harmfun1, harmfun2)
similarity.add_similarity_function(simfun1, simfun2, simfun3, simfun4) 
similarity.calculate()
```



florian-huber: Like this structure a lot!

Small remark to the scores.add_similarity_function(...) : Would be great to allow nested runs of multiple scores runs. So that a first sorting can be done by -say- simfun4 and only then run the other (computationally more expensive) scores on a subset.
jspaaks: Hi Florian,
Glad you like this proposed structure. The nested analysis that you mention would be possible under the proposed structure, like so:

```python
from matchms.harmonization import harmfun1, harmfun2 
from matchms.similarity import simfun1, simfun2, simfun3
from anotherpackage import simfun4
from matchms import Scores, Spectrum

# get spectrum Spectrum from somewhere
# get reference_spectrums Spectrum[] from somewhere

# first, calculate the similarity scores that are not so computationally intensive
scores = Scores(spectrum=spectrum, reference_spectrums=reference_spectrums)
scores.add_harmonization_function(harmfun1, harmfun2)
scores.add_similarity_function(simfun1, simfun2, simfun3) 
scores.calculate()
scores.sort_by("simfun1")

# based on the similarity scores, select a subset of interest for 
# doing a more computationally intensive analsysis, simfun4
reference_spectrums_subset = scores.top(1000).get_reference_spectrums()
scores_subset = Scores(spectrum=spectrum, reference_spectrums=reference_spectrums_subset)
scores_subset.add_harmonization_function(harmfun1, harmfun2)
scores_subset.add_similarity_function(simfun4) 
scores_subset.calculate()
```

florian-huber: That looks great and seems to fit quite well our workflow!
jspaaks: This issue can be closed.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 04281710..87e0bb32 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- Faster cosine score implementation using numba [#29](https://github.com/matchms/matchms/issues/29)
 - 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)
diff --git a/matchms/similarity/CosineGreedyNumba.py b/matchms/similarity/CosineGreedyNumba.py
new file mode 100644
index 00000000..df42d8d2
--- /dev/null
+++ b/matchms/similarity/CosineGreedyNumba.py
@@ -0,0 +1,72 @@
+from typing import Tuple
+import numpy
+from matchms.similarity.collect_peak_pairs import collect_peak_pairs
+from matchms.typing import SpectrumType
+
+
+class CosineGreedyNumba:
+    """Calculate 'cosine similarity score' between two spectra.
+
+    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 sovled in a 'greedy' way.
+    This can perform notably faster, but does occationally deviate slightly from
+    a fully correct solution (as with the Hungarian algorithm). In practice this
+    will rarely affect similarity scores notably, in particular for smaller
+    tolerances.
+    """
+    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)
+            matching_pairs = sorted(matching_pairs, key=lambda x: x[2], reverse=True)
+            return matching_pairs
+
+        def calc_score():
+            """Calculate cosine similarity score."""
+            used1 = set()
+            used2 = set()
+            score = 0.0
+            used_matches = []
+            for match in matching_pairs:
+                if not match[0] in used1 and not match[1] in used2:
+                    score += match[2]
+                    used1.add(match[0])  # Every peak can only be paired once
+                    used2.add(match[1])  # Every peak can only be paired once
+                    used_matches.append(match)
+            # Normalize score:
+            score = score/max(numpy.sum(spec1[:, 1]**2), numpy.sum(spec2[:, 1]**2))
+            return score, len(used_matches)
+
+        spec1, spec2 = get_peaks_arrays()
+        matching_pairs = get_matching_pairs()
+        return calc_score()
diff --git a/matchms/similarity/__init__.py b/matchms/similarity/__init__.py
index 33275a66..c9e57f14 100644
--- a/matchms/similarity/__init__.py
+++ b/matchms/similarity/__init__.py
@@ -1,5 +1,6 @@
 """similarity module"""
 from .CosineGreedy import CosineGreedy
+from .CosineGreedyNumba import CosineGreedyNumba
 from .CosineHungarian import CosineHungarian
 from .IntersectMz import IntersectMz
 from .ModifiedCosine import ModifiedCosine
@@ -7,6 +8,7 @@ from .ModifiedCosine import ModifiedCosine
 
 __all__ = [
     "CosineGreedy",
+    "CosineGreedyNumba",
     "CosineHungarian",
     "IntersectMz",
     "ModifiedCosine",


Test patch:
diff --git a/tests/test_cosine_greedy_numba.py b/tests/test_cosine_greedy_numba.py
new file mode 100644
index 00000000..394ea274
--- /dev/null
+++ b/tests/test_cosine_greedy_numba.py
@@ -0,0 +1,82 @@
+import numpy
+import pytest
+from matchms import Spectrum
+from matchms.similarity import CosineGreedyNumba
+
+
+def test_cosine_greedy_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_greedy_numba = CosineGreedyNumba()
+    score, n_matches = cosine_greedy_numba(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_score_greedy_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_greedy_numba = CosineGreedyNumba(tolerance=0.2)
+    score, n_matches = cosine_greedy_numba(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_score_greedy_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_greedy_numba = CosineGreedyNumba(tolerance=2.0)
+    score, n_matches = cosine_greedy_numba(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_score_greedy_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_greedy_numba = CosineGreedyNumba(tolerance=2.0)
+    score_1_2, n_matches_1_2 = cosine_greedy_numba(spectrum_1, spectrum_2)
+    score_2_1, n_matches_2_1 = cosine_greedy_numba(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."