AdrianBenson/issue36

.gitignore

@@ -186,3 +186,4 @@ docs/_build
 docs/source/api/generated
 docs/source/api/api/*
 docs/zz/*
+ARCHITECTURE.md

AUTHORS.md

@@ -7,4 +7,4 @@ order by last name) and are considered "The gSolve Developers":
 * Alison Kirkby - Earth Sciences New Zealand, New Zealand
 * Craig Miller - Earth Sciences New Zealand, New Zealand
 * Vaughan Stagpoole - Earth Sciences New Zealand, New Zealand
-* Aleksandr Spesivtsev - Earth Sciences New Zealand, New Zealand
+* Aleksandr Spesivtsev - Earth Sciences New Zealand, New Zealand

CHANGELOG.md

@@ -11,7 +11,6 @@ New features:
 * Added GPLv3 licence
 * PyGTide added for ETERNA tide correction
 
-Maintenance: 
+Maintenance:
 
 * Bug fix for terrain corrections
-

README.md

@@ -1,17 +1,15 @@
 ![gSolve logo](docs/source/_static/gsolve_logo.png)
-    
+
 [![codecov](https://codecov.io/gh/GNS-Science/gsolve/branch/main/graph/badge.svg)](https://codecov.io/gh/GNS-Science/gsolve)
 [![GitHub License](https://img.shields.io/github/license/GNS-Science/gsolve)](https://github.com/GNS-Science/gsolve/blob/main/LICENSE)
 ![Publish to PyPI](https://github.com/GNS-Science/gsolve/actions/workflows/publish.yml/badge.svg)
 [![Pypi version](https://img.shields.io/pypi/v/gsolve)](https://pypi.org/project/gsolve/)
 
-
-
 # gSolve
 
-gSolve, a Python computer library by Earth Sciences New Zealand (formerly GNS Science) to transform relative gravity survey measurements to absolute gravity values and gravity anomalies and disturbances.  
+gSolve, a Python computer library by Earth Sciences New Zealand (formerly GNS Science) to transform relative gravity survey measurements to absolute gravity values and gravity anomalies and disturbances.
 
-It is suitable for time varying gravity as well as Bouguer gravity.  
+It is suitable for time varying gravity as well as Bouguer gravity.
 
 This version is a substantial re-write of the previous python version to remove the limitation of a graphical user interface and to update to python 3 with modern software management.
 
@@ -30,8 +28,10 @@ Process gravity data for time change microgravity and Bouguer surveys.
 * option to correct for ocean loading using [pyhardisp](https://github.com/craigmillernz/pyhardisp).
 * correct for drift across loops or whole survey.
 * network adjustment with three different network adjustment algorithms depending on user requirements.
-# calibrate meters to absolute values.
-* residuals can be filtered using a percentile cut filter.  
+
+# calibrate meters to absolute values
+
+* residuals can be filtered using a percentile cut filter.
 
 ## Corrections
 
@@ -42,7 +42,7 @@ Process gravity data for time change microgravity and Bouguer surveys.
 ## Plotting
 
 * plot raw observations
-* residual cumuluative probability density functions, CDF.
+* residual cumulative probability density functions, CDF.
 * drift curve
 * network map
 
@@ -65,14 +65,14 @@ Full documentation is available here. [gSolve](https://gns-science.github.io/gso
 
 # Authors and acknowledgment
 
-gSolve builds on many previous authors.  
+gSolve builds on many previous authors.
 
 The current author team is Adrian Benson, Alison Kirkby, Craig Miller, Aleksandr Spesivtsev, Vaughan Stagpoole.
 
 This version supersedes previous Gsolve versions e.g.
-McCubbine, J., Tontini, F. C., Stagpoole, V., Smith, E., & O’Brien, G. (2018). Gsolve, a Python computer program with a graphical user interface to transform relative gravity survey measurements to absolute gravity values and gravity anomalies. SoftwareX, 7, 129–137. 
+McCubbine, J., Tontini, F. C., Stagpoole, V., Smith, E., & O’Brien, G. (2018). Gsolve, a Python computer program with a graphical user interface to transform relative gravity survey measurements to absolute gravity values and gravity anomalies. SoftwareX, 7, 129–137.
 
-# How to cite gSolve  
+# How to cite gSolve
 
 Link to JOSS paper here when it is ready.
 
@@ -82,4 +82,4 @@ Python 3.12+
 
 # License
 
-Licenced with GPLv3.
+Licensed with GPLv3.

cspell.json

@@ -15,5 +15,9 @@
         "data-science-tools",
         "git",
         "en-au"
+    ],
+    "ignorePaths": [
+        ".git/**",
+        ".gitignore"
     ]
 }

project-words.txt

@@ -1,5 +1,6 @@
 addopts
 afactor
+aimport
 airgap
 Aleksandr
 alisonk

pyproject.toml

@@ -73,6 +73,7 @@ dev = [
     "burocrata>=0.3.0",
     "ipykernel>=7.2.0",
     "ipywidgets>=8.1.8",
+    "marimo>=0.25.0",
     "pandas-stubs==2.3.3.260113",
     "pandas-vet>=2023.8.2",
     "pyrefly>=1.3.1",

src/gsolve/core/data.py

@@ -21,7 +21,6 @@
 import abc
 import copy
 import dataclasses
-import warnings
 from collections.abc import Callable
 from copy import deepcopy
 from types import MappingProxyType
@@ -42,6 +41,7 @@
     write_excel_worksheet,
 )
 from gsolve.core.utils import (
+    GSolveUserWarning,
     merge_datetime_columns,
     normalize_field_names,
     prepare_writable_df,
@@ -104,7 +104,9 @@ class DataFieldSpecification:
     DataFieldSpecification("absolute_gravity", float, default=np.nan),
 ]
 # TODO: make this a class?
-COMMON_FIELDS: dict[str, DataFieldSpecification] = {f.name: f for f in _COMMON_FIELDS}
+COMMON_FIELDS: MappingProxyType[str, DataFieldSpecification] = MappingProxyType(
+    {f.name: f for f in _COMMON_FIELDS}
+)
 
 
 @dataclasses.dataclass
@@ -503,18 +505,19 @@ def set_column(
 
     def _data_ok(self, warn: bool = True) -> bool:
         """Test whether data are complete according to specifications in ``obj._known_fields``."""
-        rval = True
+        all_ok = True
+        warner = GSolveUserWarning()
         for f in self.required_fields():
             if f not in self.data.columns:
                 if warn:
-                    warnings.warn(f"Missing required field: '{f}'")
-                rval = False
+                    warner.warn(f"Missing required field: '{f}'")
+                all_ok = False
             if self.data[f].isna().any() or self.data[f].eq("").any():
                 if warn:
-                    warnings.warn(f"Required field has empty records: '{f}'")
-                rval = False
+                    warner.warn(f"Required field has empty records: '{f}'")
+                all_ok = False
 
-        return rval
+        return all_ok
 
     @classmethod
     def from_dataframe(

src/gsolve/core/utils.py

@@ -22,7 +22,9 @@
 
 import itertools
 import sys
+import warnings
 from collections.abc import Sequence
+from pathlib import Path
 from typing import Any, Literal, TypeAliasType, get_args, get_origin, overload
 
 import numpy as np
@@ -49,6 +51,7 @@
 __all__ = [
     "DEFAULT_TIMESTAMP_COLUMNS",
     "GSolveDataWarning",
+    "GSolveUserWarning",
     "check_duplicate_index",
     "columns_to_timestamp",
     "expand_datetime_column",
@@ -926,7 +929,7 @@ def _display(self, msg: str) -> None:
     def print_msgs(self) -> None:
         """Print all stored warning messages."""
         for msg in self.messages:
-            print(f"{self.prefix}: {msg}")  # ruff: ignore[print]
+            print(f"{self.prefix}: {msg}", file=sys.stderr)  # ruff: ignore[print]
 
     def final_msg(self) -> None:
         """Print closing summary message."""
@@ -1158,3 +1161,71 @@ def convert_single_timestamp_arg(
         raise TypeError(msg)
 
     return t_
+
+
+class GSolveUserWarning:
+    """Simple class to make python warnings friendlier.
+
+    For cases where the issue is with the user's code or data.
+
+      - Ensures that a warning's trace points to the users code.
+      - Maintains a count and cache of warnings issued
+
+
+    Parameters
+    ----------
+    prefix : str, optional
+        Prepend prefix to all messages.
+    default_category : Warning, default is UserWarning
+        The category to assign warnings.
+    """
+
+    prefix: str | None
+    default_category: Warning | None
+    _count: int
+    skip_file_prefixes: tuple[str, ...]
+    _messages: list[str]
+
+    def __init__(
+        self, prefix: str | None = None, default_category: Warning = UserWarning
+    ):
+        self.prefix = prefix
+        self._count = 0
+        self.default_category = default_category
+        self.skip_file_prefixes = (str(Path(__file__).parents[1]),)
+        self._messages = []
+
+    def warn(
+        self, message: str, show_prefix: bool = True, category: Warning | None = None
+    ):
+        """
+        Issue a python warning for user errors.
+
+        Parameters
+        ----------
+        message : str
+            The warning message to issue.
+        show_prefix : bool, default True
+            Prepend message with the ``default_prefix`` set at initialisation.
+        category : Warning | None, optional
+            If specified, use ``category`` for the warning category, otherwise use
+            the category specified at initialisation (default=UserWarning).
+        """
+        if not self.prefix is None and show_prefix:
+            msg = f"{self.prefix} {message}"
+        else:
+            msg = f"{message}"
+
+        if category is None:
+            category = self.default_category
+
+        warnings.warn(
+            message=msg, category=category, skip_file_prefixes=self.skip_file_prefixes
+        )
+        self._count += 1
+        self._messages.append([msg, category])
+
+    @property
+    def count(self) -> int:
+        """The number of warnings issued by this object."""
+        return self._count

src/gsolve/gsolve_algorithms.py

@@ -23,17 +23,23 @@
 import pandas as pd
 
 from gsolve.core._typing import GSolveSolverMethod, GSolveSolverReturn
+from gsolve.core.utils import GSolveUserWarning
 from gsolve.gsolve_outputs import GSolveResults
 
 __all__ = ["GSolveSolverMethod", "call_gsolve_calibration", "call_gsolve_lstsq"]
 
+
 _GSOLVE_SOLVER_METHODS: dict[int, str] = {
     1: "Unconstrained least squares",
     2: "Partially constrained least squares",
     3: "Constrained least squares",
 }
 
 
+class GSolveSolverWarning(UserWarning):
+    """Raised when gsolve may produce unreliable results."""
+
+
 def call_gsolve_lstsq(
     obs: pd.DataFrame,
     ref_sites: pd.DataFrame,
@@ -387,13 +393,20 @@ def g_solver_lstsq(  # ruff: ignore[too-many-positional-arguments]
         residuals = b - np.dot(A, solution[:n_parameters])
 
         # Define percentile clipping interval
+        # TODO: Should this be done using residual**2 and clip upper only?
         perc = (100.0 - percentile_clipping) / 2
         ci_l = np.percentile(residuals[:n_obs], perc)
         ci_h = np.percentile(residuals[:n_obs], 100.0 - perc)
 
         # Build mask of outliers
         mask = ((residuals[:n_obs] > ci_l) & (residuals[:n_obs] < ci_h)).flatten()
-
+        _check_post_clip_data_are_ok(
+            mask=mask,
+            obs_site_id=obs_site_id,
+            ties_site_id=ties_site_id,
+            obs_loop=obs_loop,
+            use_loops=use_loops,
+        )
         # Mask outliers
         A[:n_obs, :][~mask] = 0
         b[:n_obs, :][~mask] = 0
@@ -446,3 +459,80 @@ def g_solver_lstsq(  # ruff: ignore[too-many-positional-arguments]
         calibration_factor = None
 
     return gravity, residuals, gravity_var, drift, baseline, calibration_factor, mask
+
+
+def _check_post_clip_data_are_ok(
+    mask: np.ndarray[tuple[int]],
+    obs_site_id: np.ndarray[tuple[int]],
+    ties_site_id: np.ndarray[tuple[int]],
+    obs_loop: np.ndarray[tuple[int]],
+    use_loops: bool,
+) -> bool:
+    """Check if percentile clipping does damage."""
+    warner = GSolveUserWarning(
+        prefix="After percentile clipping:", default_category=GSolveSolverWarning
+    )
+    obs_site_id_remain = obs_site_id[mask]
+    obs_loop_remain = obs_loop[mask]
+
+    # warn about site removal
+    if len(dropped_sites := np.setdiff1d(obs_site_id, obs_site_id_remain)) > 0:
+        warner.warn(f"sites were completely removed: {dropped_sites}")
+
+    # ensure all tie sites made it
+    dropped_ties = np.intersect1d(ties_site_id, obs_site_id_remain)
+    if len(dropped_ties) == 0:
+        msg = (
+            "All reference sites were completely removed due to "
+            f"percentile clipping: {dropped_ties}"
+        )
+        raise ValueError(msg)
+
+    if len(dropped_ties) < len(ties_site_id):
+        warner.warn(
+            f"{len(dropped_ties)} of {len(ties_site_id)} reference sites "
+            f"were completely: {dropped_ties})"
+        )
+
+    # if using loops and there was actually more than 1 loop
+    # - check that loops were not completely removed
+    # - check that loops have common stations or include a ref site
+
+    # if use_loops and have > 1 loops, then check_loops
+    check_loops = use_loops and not (obs_loop[0] == obs_loop).all()
+    if check_loops:
+        if len(dropped_loops := np.setdiff1d(obs_loop, obs_loop_remain)) > 0:
+            warner.warn(f"some loops were completely removed: {dropped_loops}")
+
+        # check if loops are now isolated from rest of survey
+        for loop_id in np.unique(obs_loop_remain):
+            m_pre = obs_loop == loop_id
+            m_post = obs_loop_remain == loop_id
+
+            in_other_loops_pre = np.intersect1d(obs_site_id[m_pre], obs_site_id[~m_pre])
+            in_other_loops_pre = len(in_other_loops_pre) > 0
+            in_other_loops_post = np.intersect1d(
+                obs_site_id_remain[m_post], obs_site_id_remain[~m_post]
+            )
+            in_other_loops_post = len(in_other_loops_post) > 0
+
+            has_ref_gravity_pre = np.isin(ties_site_id, obs_site_id[m_pre]).any()
+            has_ref_gravity_post = np.isin(
+                ties_site_id, obs_site_id_remain[m_post]
+            ).any()
+
+            isolated_loop_pre = not in_other_loops_pre and not has_ref_gravity_pre
+            isolated_loop_post = not in_other_loops_post and not has_ref_gravity_post
+
+            if isolated_loop_pre:
+                # Already isolated so, nothing will have changed
+                continue
+            if isolated_loop_post:
+                warner.warn(
+                    f"loop '{loop_id}' has no sites in common with the rest of survey"
+                )
+
+            if has_ref_gravity_pre and not has_ref_gravity_post:
+                warner.warn(f"all reference sites removed from loop '{loop_id}'")
+
+    return warner.count == 0