Repository: matchms/matchms
Pull request #98 (2020-07-08 10:11:31)

Cosine extension
- fix wrong denominator in cosine score calculation
- add options to raise peaks to given power

jspaaks: What a magnificent PR, thank you!
jspaaks: Sure no problem. Thank you for reviewing. Should I merge now?

jspaaks: Please go ahead and merge.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 21d36370..cfe2d073 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,13 +9,18 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
-- Added similarity score based on comparing parent masses [#79](https://github.com/matchms/matchms/pull/79)
+- Peak weighting option for cosine score and modified cosine score [#96](https://github.com/matchms/matchms/issues/96)
+- Similarity score based on comparing parent masses [#79](https://github.com/matchms/matchms/pull/79)
 
 ### Changed
 
 - CosineGreedy function is now numba based [#86](https://github.com/matchms/matchms/pull/86)
 - Extended readthedocs documentation [#82](https://github.com/matchms/matchms/issues/82)
 
+### Fixed
+
+- Incorrect denominator for cosine score normalization [#98](https://github.com/matchms/matchms/pull/98)
+
 ## [0.4.0] - 2020-06-11
 
 ### Added
diff --git a/matchms/Scores.py b/matchms/Scores.py
index 963720d3..fdc38825 100644
--- a/matchms/Scores.py
+++ b/matchms/Scores.py
@@ -45,9 +45,9 @@ class Scores:
     .. testoutput::
 
         Cosine score between spectrum1 and spectrum3 is 0.00 with 0 matched peaks
-        Cosine score between spectrum1 and spectrum4 is 0.68 with 3 matched peaks
-        Cosine score between spectrum2 and spectrum3 is 0.10 with 1 matched peaks
-        Cosine score between spectrum2 and spectrum4 is 0.33 with 1 matched peaks
+        Cosine score between spectrum1 and spectrum4 is 0.80 with 3 matched peaks
+        Cosine score between spectrum2 and spectrum3 is 0.14 with 1 matched peaks
+        Cosine score between spectrum2 and spectrum4 is 0.61 with 1 matched peaks
     """
     def __init__(self, references: ReferencesType, queries: QueriesType, similarity_function: SimilarityFunction):
         """
diff --git a/matchms/Spectrum.py b/matchms/Spectrum.py
index 36bcb9d4..d4bcab55 100644
--- a/matchms/Spectrum.py
+++ b/matchms/Spectrum.py
@@ -117,7 +117,7 @@ class Spectrum:
             def make_stems():
                 """calculate where the stems of the spectrum peaks are going to be"""
                 x = numpy.empty([2, self.peaks.mz.size], dtype="float")
-                y = numpy.empty_like(x)
+                y = numpy.empty(x.shape)
                 for i, mz in enumerate(self.peaks.mz):
                     x[0:2, i] = [mz, mz]
                     y[0:2, i] = [0, self.peaks.intensities[i]]
diff --git a/matchms/calculate_scores.py b/matchms/calculate_scores.py
index d8cad33b..30036e02 100644
--- a/matchms/calculate_scores.py
+++ b/matchms/calculate_scores.py
@@ -36,8 +36,8 @@ def calculate_scores(references: ReferencesType,
     .. testoutput::
 
         Cosine score between spectrum1 and spectrum1 is 1.00 with 3 matched peaks
-        Cosine score between spectrum1 and spectrum2 is 0.52 with 1 matched peaks
-        Cosine score between spectrum2 and spectrum1 is 0.52 with 1 matched peaks
+        Cosine score between spectrum1 and spectrum2 is 0.83 with 1 matched peaks
+        Cosine score between spectrum2 and spectrum1 is 0.83 with 1 matched peaks
         Cosine score between spectrum2 and spectrum2 is 1.00 with 3 matched peaks
 
     Parameters
diff --git a/matchms/similarity/CosineGreedy.py b/matchms/similarity/CosineGreedy.py
index 85cdca92..b468c541 100644
--- a/matchms/similarity/CosineGreedy.py
+++ b/matchms/similarity/CosineGreedy.py
@@ -42,17 +42,25 @@ class CosineGreedy:
 
     .. testoutput::
 
-        Cosine score is 0.52 with 1 matched peaks
+        Cosine score is 0.83 with 1 matched peaks
 
     """
-    def __init__(self, tolerance: float = 0.1):
+    def __init__(self, tolerance: float = 0.1, mz_power: float = 0.0,
+                 intensity_power: float = 1.0):
         """
         Parameters
         ----------
-        tolerance
-            Peaks will be considered a match when <= tolerance apart.
+        tolerance:
+            Peaks will be considered a match when <= tolerance apart. Default is 0.1.
+        mz_power:
+            The power to raise m/z to in the cosine function. The default is 0, in which
+            case the peak intensity products will not depend on the m/z ratios.
+        intensity_power:
+            The power to raise intensity to in the cosine function. The default is 1.
         """
         self.tolerance = tolerance
+        self.mz_power = mz_power
+        self.intensity_power = intensity_power
 
     def __call__(self, spectrum1: SpectrumType, spectrum2: SpectrumType) -> Tuple[float, int]:
         """Calculate cosine score between two spectra.
@@ -71,11 +79,14 @@ class CosineGreedy:
         """
         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 = collect_peak_pairs(spec1, spec2, self.tolerance, shift=0.0,
+                                                mz_power=self.mz_power,
+                                                intensity_power=self.intensity_power)
             matching_pairs = sorted(matching_pairs, key=lambda x: x[2], reverse=True)
             return matching_pairs
 
         spec1 = get_peaks_array(spectrum1)
         spec2 = get_peaks_array(spectrum2)
         matching_pairs = get_matching_pairs()
-        return score_best_matches(matching_pairs, spec1, spec2)
+        return score_best_matches(matching_pairs, spec1, spec2,
+                                  self.mz_power, self.intensity_power)
diff --git a/matchms/similarity/ModifiedCosine.py b/matchms/similarity/ModifiedCosine.py
index 955a8288..77870669 100644
--- a/matchms/similarity/ModifiedCosine.py
+++ b/matchms/similarity/ModifiedCosine.py
@@ -43,17 +43,25 @@ class ModifiedCosine:
 
     .. testoutput::
 
-        Modified cosine score is 0.52 with 1 matched peaks
+        Modified cosine score is 0.83 with 1 matched peaks
 
     """
-    def __init__(self, tolerance: float = 0.1):
+    def __init__(self, tolerance: float = 0.1, mz_power: float = 0.0,
+                 intensity_power: float = 1.0):
         """
-        Args:
-        -----
-        tolerance
-            Peaks will be considered a match when <= tolerance a part. Default is 0.1.
+        Parameters
+        ----------
+        tolerance:
+            Peaks will be considered a match when <= tolerance apart. Default is 0.1.
+        mz_power:
+            The power to raise mz to in the cosine function. The default is 0, in which
+            case the peak intensity products will not depend on the m/z ratios.
+        intensity_power:
+            The power to raise intensity to in the cosine function. The default is 1.
         """
         self.tolerance = tolerance
+        self.mz_power = mz_power
+        self.intensity_power = intensity_power
 
     def __call__(self, spectrum1: SpectrumType, spectrum2: SpectrumType) -> Tuple[float, int]:
         """Calculate modified cosine score between two spectra.
@@ -71,15 +79,20 @@ class ModifiedCosine:
         """
         def get_matching_pairs():
             """Find all pairs of peaks that match within the given tolerance."""
-            zero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=0.0)
+            zero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=0.0,
+                                            mz_power=self.mz_power,
+                                            intensity_power=self.intensity_power)
             message = "Precursor_mz missing. Apply 'add_precursor_mz' filter first."
             assert spectrum1.get("precursor_mz") and spectrum2.get("precursor_mz"), message
             mass_shift = spectrum1.get("precursor_mz") - spectrum2.get("precursor_mz")
-            nonzero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=mass_shift)
+            nonzero_pairs = collect_peak_pairs(spec1, spec2, self.tolerance, shift=mass_shift,
+                                               mz_power=self.mz_power,
+                                               intensity_power=self.intensity_power)
             unsorted_matching_pairs = zero_pairs + nonzero_pairs
             return sorted(unsorted_matching_pairs, key=lambda x: x[2], reverse=True)
 
         spec1 = get_peaks_array(spectrum1)
         spec2 = get_peaks_array(spectrum2)
         matching_pairs = get_matching_pairs()
-        return score_best_matches(matching_pairs, spec1, spec2)
+        return score_best_matches(matching_pairs, spec1, spec2,
+                                  self.mz_power, self.intensity_power)
diff --git a/matchms/similarity/spectrum_similarity_functions.py b/matchms/similarity/spectrum_similarity_functions.py
index b1a99574..2b2b7dd7 100644
--- a/matchms/similarity/spectrum_similarity_functions.py
+++ b/matchms/similarity/spectrum_similarity_functions.py
@@ -5,7 +5,9 @@ from matchms.typing import SpectrumType
 
 
 @numba.njit
-def collect_peak_pairs(spec1, spec2, tolerance, shift=0):
+def collect_peak_pairs(spec1, spec2, tolerance, shift=0,
+                       mz_power=0.0, intensity_power=1.0):
+    # pylint: disable=too-many-arguments
     """Find matching pairs between two spectra.
 
     Args
@@ -18,6 +20,11 @@ def collect_peak_pairs(spec1, spec2, tolerance, shift=0):
         Peaks will be considered a match when <= tolerance appart.
     shift : float, optional
         Shift spectra peaks by shift. The default is 0.
+    mz_power: float, optional
+        The power to raise mz to in the cosine function. The default is 0, in which
+        case the peak intensity products will not depend on the m/z ratios.
+    intensity_power: float, optional
+        The power to raise intensity to in the cosine function. The default is 1.
 
     Returns
     -------
@@ -28,9 +35,12 @@ def collect_peak_pairs(spec1, spec2, tolerance, shift=0):
 
     for idx in range(len(spec1)):
         intensity = spec1[idx, 1]
+        mz = spec1[idx, 0]
         matches = numpy.where((numpy.abs(spec2[:, 0] - spec1[idx, 0] + shift) <= tolerance))[0]
         for match in matches:
-            matching_pairs.append((idx, match, intensity*spec2[match][1]))
+            power_prod_spec1 = ((mz ** mz_power) * (intensity ** intensity_power))
+            power_prod_spec2 = ((spec2[match][0] ** mz_power) * (spec2[match][1] ** intensity_power))
+            matching_pairs.append((idx, match, power_prod_spec1 * power_prod_spec2))
 
     return matching_pairs
 
@@ -44,7 +54,8 @@ def get_peaks_array(spectrum: SpectrumType) -> numpy.ndarray:
 
 
 def score_best_matches(matching_pairs: list, spec1: numpy.ndarray,
-                       spec2: numpy.ndarray) -> Tuple[float, int]:
+                       spec2: numpy.ndarray, mz_power: float = 0.0,
+                       intensity_power: float = 1.0) -> Tuple[float, int]:
     """Calculate cosine-like score by multiplying matches. Does recuire a sorted
     list of matching peaks (sorted by intensity product)."""
     used1 = set()
@@ -58,5 +69,8 @@ def score_best_matches(matching_pairs: list, spec1: numpy.ndarray,
             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))
+    spec1_power = numpy.power(spec1[:, 0], mz_power) * numpy.power(spec1[:, 1], intensity_power)
+    spec2_power = numpy.power(spec2[:, 0], mz_power) * numpy.power(spec2[:, 1], intensity_power)
+
+    score = score/(numpy.sqrt(numpy.sum(spec1_power**2)) * numpy.sqrt(numpy.sum(spec2_power**2)))
     return score, len(used_matches)
diff --git a/readthedocs/index.rst b/readthedocs/index.rst
index 06b2268f..f5bbfbed 100644
--- a/readthedocs/index.rst
+++ b/readthedocs/index.rst
@@ -77,14 +77,14 @@ Should output
     Added adduct M-H to metadata.
     ...
     Reference scan id: 1320
-    Query scan id: 2833
-    Score: 0.0137
-    Number of matching peaks: 20
-    ----------------------------
-    Reference scan id: 1320
     Query scan id: 1372
-    Score: 0.8379
+    Score: 0.9143
     Number of matching peaks: 25
+    ----------------------------
+    Reference scan id: 2833
+    Query scan id: 1320
+    Score: 0.0144
+    Number of matching peaks: 20
     ...
 
 Indices and tables


Test patch:
diff --git a/integration-tests/test_user_workflow.py b/integration-tests/test_user_workflow.py
index df457593..0f613c78 100644
--- a/integration-tests/test_user_workflow.py
+++ b/integration-tests/test_user_workflow.py
@@ -52,15 +52,15 @@ def test_user_workflow():
     actual_top10 = sorted_by_score[:10]
 
     expected_top10 = [
-        (references[48], queries[50], pytest.approx(0.9994510368270997, rel=1e-9), 25),
-        (references[50], queries[48], pytest.approx(0.9994510368270997, rel=1e-9), 25),
-        (references[46], queries[48], pytest.approx(0.9981252309590571, rel=1e-9), 27),
-        (references[48], queries[46], pytest.approx(0.9981252309590571, rel=1e-9), 27),
-        (references[46], queries[50], pytest.approx(0.9979632203390496, rel=1e-9), 22),
-        (references[50], queries[46], pytest.approx(0.9979632203390496, rel=1e-9), 22),
-        (references[73], queries[74], pytest.approx(0.9956795920716246, rel=1e-9), 23),
-        (references[74], queries[73], pytest.approx(0.9956795920716246, rel=1e-9), 23),
-        (references[57], queries[59], pytest.approx(0.9886557001269415, rel=1e-9), 46),
-        (references[59], queries[57], pytest.approx(0.9886557001269415, rel=1e-9), 46),
+        (references[48], queries[50], pytest.approx(0.9994783627790965, rel=1e-9), 25),
+        (references[50], queries[48], pytest.approx(0.9994783627790965, rel=1e-9), 25),
+        (references[46], queries[48], pytest.approx(0.9990141860269471, rel=1e-9), 27),
+        (references[48], queries[46], pytest.approx(0.9990141860269471, rel=1e-9), 27),
+        (references[46], queries[50], pytest.approx(0.9988793406908719, rel=1e-9), 22),
+        (references[50], queries[46], pytest.approx(0.9988793406908719, rel=1e-9), 22),
+        (references[57], queries[59], pytest.approx(0.9982171275552505, rel=1e-9), 46),
+        (references[59], queries[57], pytest.approx(0.9982171275552505, rel=1e-9), 46),
+        (references[73], queries[74], pytest.approx(0.9973823244169199, rel=1e-9), 23),
+        (references[74], queries[73], pytest.approx(0.9973823244169199, rel=1e-9), 23),
     ]
     assert actual_top10 == expected_top10
diff --git a/tests/test_cosine_greedy.py b/tests/test_cosine_greedy.py
index 9cadaccf..5dab196c 100644
--- a/tests/test_cosine_greedy.py
+++ b/tests/test_cosine_greedy.py
@@ -18,7 +18,9 @@ def test_cosine_greedy_without_parameters():
     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()
+    denominator = numpy.sqrt((spectrum_1.peaks.intensities ** 2).sum()) \
+        * numpy.sqrt((spectrum_2.peaks.intensities ** 2).sum())
+    expected_score = multiply_matching_intensities.sum() / denominator
 
     assert score == pytest.approx(expected_score, 0.0001), "Expected different cosine score."
     assert n_matches == len(expected_matches), "Expected different number of matching peaks."
@@ -38,7 +40,9 @@ def test_cosine_score_greedy_with_tolerance_0_2():
     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()
+    denominator = numpy.sqrt((spectrum_1.peaks.intensities ** 2).sum()) \
+        * numpy.sqrt((spectrum_2.peaks.intensities ** 2).sum())
+    expected_score = multiply_matching_intensities.sum() / denominator
 
     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."
@@ -58,7 +62,9 @@ def test_cosine_score_greedy_with_tolerance_2_0():
     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()
+    denominator = numpy.sqrt((spectrum_1.peaks.intensities ** 2).sum()) \
+        * numpy.sqrt((spectrum_2.peaks.intensities ** 2).sum())
+    expected_score = multiply_matching_intensities.sum() / denominator
 
     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."
@@ -80,3 +86,33 @@ def test_cosine_score_greedy_order_of_arguments():
 
     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_greedy_with_peak_powers():
+    """Compare output cosine score with own calculation on simple dummy spectrums.
+    Here testing the options to raise peak intensities to given powers.
+    """
+    mz_power = 0.5
+    intensity_power = 2.0
+    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 = CosineGreedy(tolerance=1.0, mz_power=mz_power, intensity_power=intensity_power)
+    score, n_matches = cosine_greedy(spectrum_1, spectrum_2)
+
+    # Derive expected cosine score
+    matches = [0, 1, 4]  # Those peaks have matching mz values (within given tolerance)
+    intensity1 = spectrum_1.peaks.intensities
+    mz1 = spectrum_1.peaks.mz
+    intensity2 = spectrum_2.peaks.intensities
+    mz2 = spectrum_2.peaks.mz
+    multiply_matching_intensities = (mz1[matches] ** mz_power) * (intensity1[matches] ** intensity_power) \
+        * (mz2[matches] ** mz_power) * (intensity2[matches] ** intensity_power)
+    denominator = numpy.sqrt((((mz1 ** mz_power) * (intensity1 ** intensity_power)) ** 2).sum()) \
+        * numpy.sqrt((((mz2 ** mz_power) * (intensity2 ** intensity_power)) ** 2).sum())
+    expected_score = multiply_matching_intensities.sum() / denominator
+
+    assert score == pytest.approx(expected_score, 0.0001), "Expected different cosine score."
+    assert n_matches == len(matches), "Expected different number of matching peaks."
diff --git a/tests/test_modified_cosine.py b/tests/test_modified_cosine.py
index bbaed356..5bff822d 100644
--- a/tests/test_modified_cosine.py
+++ b/tests/test_modified_cosine.py
@@ -58,7 +58,7 @@ def test_modified_cosine_with_mass_shift_5_tolerance_2():
     modified_cosine = ModifiedCosine(tolerance=2.0)
     score, n_matches = modified_cosine(norm_spectrum_1, norm_spectrum_2)
 
-    assert score == pytest.approx(0.903412, 0.0001), "Expected different modified cosine score."
+    assert score == pytest.approx(0.96788, 0.0001), "Expected different modified cosine score."
     assert n_matches == 6, "Expected 6 matching peaks."