Resolve "Modematching analyzing script"

.gitignore

@@ -23,3 +23,4 @@ src/tests/test_io/gzip_data_files/
 *.gz
 
 \.dmypy\.json
+TestCSV*

.gitlab-ci.yml

@@ -53,7 +53,7 @@ test:
       - src/tests/htmlcov
     expire_in: 30 days
 
-codequality:
+.codequality:
   image: lasnq/openqlab
   stage: test
   allow_failure: true

.pylintrc

@@ -313,7 +313,7 @@ indent-after-paren=4
 indent-string='    '
 
 # Maximum number of characters on a single line.
-max-line-length=160
+max-line-length=100
 
 # Maximum number of lines in a module.
 max-module-lines=1000

doc/analysis.rst

@@ -16,6 +16,14 @@ Fit beam data obtained from a beam analyzer to the gaussian beam model using non
 .. automodule:: openqlab.analysis.gaussian_beam
     :members:
 
+Mode Matching
+-------------
+
+Mode matching calculations
+
+.. automodule:: openqlab.analysis.modematching
+    :members:
+
 Phase
 -----
 

setup.cfg

@@ -89,6 +89,6 @@ version = 3.0.3
 package = openqlab
 
 [pycodestyle]
-max-line-length = 160
+max-line-length = 100
 statistics = True
 exclude = OldImporters

src/openqlab/analysis/modematching.py

+"""Automatically calculate modematchings with data taken from an oscilloscope."""
+
+from typing import Optional
+import logging as log
+from pandas import Series
+from scipy.signal import find_peaks
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def calculate(data: Series,  # pylint: disable=invalid-name
+              plot: bool = False,
+              U_max: Optional[float] = None,
+              offset: Optional[float] = None, rel_prominence: float = .02) -> float:
+    """Calculate the mode matching.
+
+    It assumes a cavity scan bounded by two peaks of the main mode.
+    The method looks for the smaller peaks where the detection threshold
+    can be adjusted with :obj:`rel_prominence`.
+
+    Offset
+        The default method to find out the offset is by calculating the mean value.
+        If you have measured it more precisely, use the parameter :obj:`offset`.
+
+    Improve precision
+        To get a better resolution for small peaks there is an option to take data
+        with a clipped main mode. Use the parameter :obj:`U_max` to manually set
+        the measured maximum value.
+
+    Parameters
+    ----------
+
+    data : Series
+        Measured data (just one column).
+    plot : bool
+        Make a plot to see, if the correct peaks where detected.
+    U_max : Optional[float]
+        U_max is the parameter to set the peak voltage of the clipped main peak.
+    rel_prominence : float
+        rel_prominence is the parameter to adjust the threshold for the detection
+        of small peaks.
+
+    Returns
+    -------
+
+    float
+        Calculated mode matching value.
+    """
+    if len(data.shape) != 1:
+        raise ValueError('The DataFrame should only contain one single column.')
+    data = data.dropna()
+
+    # Adjust offset
+    if offset is None:
+        offset = np.median(data)
+    data -= offset
+
+    # Make peaks positive if necessary
+    _adjust_peak_sign(data)
+
+    # Find highest value
+    max_value = np.max(data)
+    if U_max is None:
+        U_max = max_value
+    else:
+        U_max = abs(U_max - offset)
+
+    # Find peaks
+    peaks, peak_dict = find_peaks(data, prominence=U_max*rel_prominence)
+
+    # Find occurences of the main mode. Two are assumed
+    main_mode = np.where(peak_dict['prominences'] >= max_value*.9)[0]
+    if len(main_mode) != 2:
+        raise ValueError('The main mode must occur exactly two times for the algorithm to work, '
+                         f'but it found {len(main_mode)} main modes.')
+
+    # Main peak voltage
+    log.info(f'U_max: {U_max}')
+    # Sum of all different modes (excluding 2nd main mode)
+    U_sum = sum(data.iloc[peaks[main_mode[0]+1:main_mode[1]]], U_max)  # pylint: disable=invalid-name
+    # This version with U_max makes it possible to manually
+    # include a clipped value for the main peak
+    log.info(f'U_sum: {U_sum}')
+
+    # Mode matching
+    modematching = U_max / U_sum
+
+    # Plotting
+    if plot:
+        index_first, index_last = peaks[main_mode]
+        axes = data.plot()
+        plt.axvline(x=data.index[index_first], color='gray')
+        plt.axvline(x=data.index[index_last], color='gray')
+        plt.axvspan(data.index[0], data.index[index_first], color='gray', alpha=0.5)
+        plt.axvspan(data.index[index_last], data.index[-1], color='gray', alpha=0.5)
+        data.iloc[peaks].plot(style='.')
+        data.iloc[peaks[main_mode]].plot(style='o')
+
+        axes.set_xlim(data.index[0], data.index[-1])
+
+    return modematching
+
+
+def _adjust_peak_sign(data: Series):
+    minimum = np.min(data)
+    maximum = np.max(data)
+
+    if abs(minimum) > abs(maximum):
+        data *= -1

src/tests/test_modematching.py

+import unittest
+import matplotlib.pyplot as plt
+from openqlab import analysis, io
+from openqlab.analysis.modematching import calculate
+
+
+class TestModematching(unittest.TestCase):
+    def setUp(self):
+        self.data = io.read('modematching/opo_scan.csv')
+
+    def test_default_inverted(self):
+        result = calculate(self.data['opo_scan_2'], plot=True)
+        # plt.show()
+        self.assertAlmostEqual(result, .955, places=2)
+
+    def test_default_not_inverted(self):
+        self.data *= -1
+        result = calculate(self.data['opo_scan_2'], plot=True)
+        self.assertAlmostEqual(result, .955, places=2)
+
+    def test_to_many_columns(self):
+        with self.assertRaisesRegex(ValueError, 'only contain one single'):
+            calculate(self.data)
+
+    def test_with_clipped_data(self):
+        self.data.clip(lower=-1, upper=1, inplace=True)
+        result = calculate(self.data['opo_scan_2'], U_max=6.668)
+        self.assertAlmostEqual(result, .955, places=2)
+
+    def test_fail_with_only_one_peak(self):
+        with self.assertRaises(ValueError) as err:
+            calculate(self.data['opo_scan_2'].loc[:.01])
+        self.assertEqual(err.exception.args[0],
+                'The main mode must occur exactly two times for the algorithm to work, but it found 1 main modes.')
+
+    def test_with_much_offset(self):
+        self.data *= -1
+        self.data -= 2
+        result = calculate(self.data['opo_scan_2'], U_max=6.668-2, offset=-1.98, plot=True)
+        self.assertAlmostEqual(result, .955, places=2)

watch_tests.local.sh

@@ -2,5 +2,5 @@
 while true
 do
   inotifywait -e modify src/**/*.py
-  make mypy && make test && make pylint
+  make mypy && make test
 done