Merge branch 'main' into feat/fix_material_descr_rect

.github/workflows/pre-commit.yml

@@ -0,0 +1,16 @@
+name: pre-commit
+
+on:
+  pull_request:
+  push:
+    branches: [main]
+
+jobs:
+  pre-commit:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v3
+    - uses: actions/setup-python@v3
+      with:
+        python-version: 3.10.14
+    - uses: pre-commit/[email protected]

.github/workflows/pytest.yml

@@ -0,0 +1,23 @@
+name: pytest
+
+on:
+  pull_request:
+  push:
+    branches: [main]
+
+jobs:
+  pytest:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v4
+    - name: Create Environment and run tests
+      shell: bash
+      run: |
+        wget -O Miniforge3.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh"
+        bash Miniforge3.sh -b -p "${HOME}/conda"
+        source "${HOME}/conda/etc/profile.d/conda.sh"
+        conda init --all
+        source "${HOME}/.bash_profile"
+        conda env create -f environment-dev.yml
+        conda activate boreholes-dev
+        pytest

src/stratigraphy/util/dataclasses.py

@@ -1,5 +1,7 @@
 """Dataclasses for stratigraphy module."""
 
+from __future__ import annotations
+
 from dataclasses import dataclass
 
 import numpy as np
@@ -16,6 +18,9 @@ class Point:
     def tuple(self) -> (float, float):
         return self.x, self.y
 
+    def distance_to(self, point: Point) -> float:
+        return np.sqrt((self.x - point.x) ** 2 + (self.y - point.y) ** 2)
+
 
 @dataclass
 class Line:
@@ -32,6 +37,7 @@ def __post_init__(self):
 
         self.slope = self.slope()
         self.intercept = self.intercept()
+        self.length = self.start.distance_to(self.end)
 
     def distance_to(self, point: Point) -> float:
         # Calculate the distance of the point to the line:

src/stratigraphy/util/geometric_line_utilities.py

@@ -237,26 +237,36 @@ def merge_parallel_lines(lines: list[Line], tol: int, angle_threshold: float) ->
         return merged_lines
 
 
-def _odr_regression(x: ArrayLike, y: ArrayLike) -> tuple:
-    """Perform orthogonal distance regression on the given data.
+def _odr_regression(x: ArrayLike, y: ArrayLike, weights: ArrayLike = None) -> tuple:
+    """Perform orthogonal distance regression (a.k.a. Deming regression) on the given data.
+
+    This algorithm minimises the quadratic distance from the given points to the resulting line, where each point
+    can be given a specific weight. If no weights are specified, then all points are weighted equally.
+
+    The implementation follow the paper "A tutorial on the total least squares method for fitting a straight line and
+    a plane" (https://www.researchgate.net/publication/272179120).
 
     Note: If the problem is ill-defined (i.e. denominator == nominator == 0),
     then the function will return phi=np.nan and r=np.nan.
 
     Args:
         x (ArrayLike): The x-coordinates of the data.
         y (ArrayLike): The y-coordinates of the data.
+        weights (ArrayLike, optional): The weight for each data point. Defaults to None.
 
     Returns:
         tuple: (phi, r), the best fit values for the line equation in normal form.
     """
-    x_mean = np.mean(x)
-    y_mean = np.mean(y)
-    nominator = -2 * np.sum((x - x_mean) * (y - y_mean))
-    denominator = np.sum((y - y_mean) ** 2 - (x - x_mean) ** 2)
+    if weights is None:
+        weights = np.ones((len(x),))
+
+    x_mean = np.mean(np.dot(weights, x)) / np.sum(weights)
+    y_mean = np.mean(np.dot(weights, y)) / np.sum(weights)
+    nominator = -2 * np.sum(np.dot(weights, (x - x_mean) * (y - y_mean)))
+    denominator = np.sum(np.dot(weights, (y - y_mean) ** 2 - (x - x_mean) ** 2))
     if nominator == 0 and denominator == 0:
         logger.warning(
-            "The problem is ill defined as both nominator and denominator for arctan are 0. "
+            "The line merging problem is ill defined as both nominator and denominator for arctan are 0. "
             "We return phi=np.nan and r=np.nan."
         )
         return np.nan, np.nan
@@ -289,7 +299,14 @@ def _merge_lines(line1: Line, line2: Line) -> Line | None:
     """
     x = np.array([line1.start.x, line1.end.x, line2.start.x, line2.end.x])
     y = np.array([line1.start.y, line1.end.y, line2.start.y, line2.end.y])
-    phi, r = _odr_regression(x, y)
+    # For the orthogonal distance regression, we weigh each point proportionally to the length of the line. The
+    # intuition behind this, is that we want the resulting merged line to be equal, regardless of whether we are
+    # merging a line AB with a line CD, or whether we are merging three lines AB, CE and ED, where E is the midpoint
+    # of CD. The current choice of weights only achieves this requirement if the lines are exactly parallel. Still,
+    # we use these weights also for non-parallel input lines, as it is still a good enough approximation for our
+    # purposes, and it keeps the mathematics reasonable simple.
+    weights = np.array([line1.length, line1.length, line2.length, line2.length])
+    phi, r = _odr_regression(x, y, weights)
     if np.isnan(phi) or np.isnan(r):
         return None
 

tests/test_geometric_line_utilities.py

@@ -14,15 +14,29 @@
 @pytest.fixture(
     params=[
         # Test case 1: horizontal line at y=1
-        (np.array([0, 1, 2, 3]), np.array([1, 1, 1, 1]), np.pi / 2, 1),
+        (np.array([0, 1, 2, 3]), np.array([1, 1, 1, 1]), np.array([1, 1, 1, 1]), np.pi / 2, 1),
         # Test case 2: 45 degrees through zero
-        (np.array([0, 1, 2, 3]), np.array([0, 1, 2, 3]), -np.pi / 4, 0),
+        (np.array([0, 1, 2, 3]), np.array([0, 1, 2, 3]), np.array([1, 1, 1, 1]), -np.pi / 4, 0),
         # Test case 3: vertical line at x=2
-        (np.array([2, 2, 2, 2]), np.array([0, 1, 2, 3]), 0, 2),
+        (np.array([2, 2, 2, 2]), np.array([0, 1, 2, 3]), np.array([1, 1, 1, 1]), 0, 2),
         # Test case 4: best fix is a horizontal line at y=0.5
-        (np.array([0, 0, 2, 2]), np.array([0, 1, 0, 1]), np.pi / 2, 0.5),
+        (np.array([0, 0, 2, 2]), np.array([0, 1, 0, 1]), np.array([1, 1, 1, 1]), np.pi / 2, 0.5),
         # Test case 4: best fix is a vertical line at x=0.5
-        (np.array([0, 1, 0, 1]), np.array([0, 0, 2, 2]), 0, 0.5),
+        (np.array([0, 1, 0, 1]), np.array([0, 0, 2, 2]), np.array([1, 1, 1, 1]), 0, 0.5),
+        (
+            np.array([0, 1, 2, 3]),
+            np.array([0, 1, 1, 0]),
+            np.array([3, 1, 1, 3]),
+            np.pi / 2,
+            1 / 4,
+        ),  # test impact of the weights (horizontal, parallel lines)
+        (
+            np.array([0, 0, 1]),
+            np.array([0.1, -0.1, 0]),
+            np.array([2, 1, 1]),
+            1.536249331404415,
+            0.03362011376179206,
+        ),  # test impact of the weights (three points)
     ]
 )
 def odr_regression_case(request):  # noqa: D103
@@ -31,8 +45,8 @@ def odr_regression_case(request):  # noqa: D103
 
 # Use the fixture in the test function
 def test_odr_regression(odr_regression_case):  # noqa: D103
-    x, y, expected_phi, expected_r = odr_regression_case
-    phi, r = _odr_regression(x, y)
+    x, y, weights, expected_phi, expected_r = odr_regression_case
+    phi, r = _odr_regression(x, y, weights)
     assert pytest.approx(expected_phi) == phi
     assert pytest.approx(expected_r) == r
 
@@ -49,7 +63,7 @@ def odr_regression_with_zeroes_case(request):  # noqa: D103
 
 def test_odr_regression_with_zeroes(odr_regression_with_zeroes_case):  # noqa: D103
     x, y = odr_regression_with_zeroes_case
-    phi, r = _odr_regression(x, y)
+    phi, r = _odr_regression(x, y, np.array([1, 1, 1, 1]))
     print(f"phi: {phi}, r: {r}")
     assert np.isnan(phi)  # Expected value
     assert np.isnan(r)  # Expected value