From 4a6b30fd2246cf66429380a1eafe9ec312435a47 Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Wed, 9 Sep 2026 18:21:21 +0300 Subject: [PATCH 1/4] Add release and validation workflows, and a release guide Set up the GitHub Actions and documentation needed to run releases through GitHub and Zenodo automatically, following the same pattern as LayTracer. PyNetDesign is not published to PyPI, so the packaging jobs stop after building and checking the distributions. - add validate-zenodo.yml, which validates .zenodo.json on every push and pull request to main and dev, and is reusable from other workflows - add validate-tag.yml, a reusable workflow enforcing strict vMAJOR.MINOR.PATCH tags that point at a commit contained in main - add pytest.yml, running pytest on Python 3.9 through 3.12 - add release-on-published.yml, which validates the tag and the Zenodo metadata, then builds the distributions and asserts that the built version matches the tag - add RELEASE.md documenting the versioning rules, the release flow, the Zenodo integration and how the pre-integration v1.0.0 record is cross-linked to the new concept record --- .github/workflows/pytest.yml | 35 +++++++ .github/workflows/release-on-published.yml | 75 ++++++++++++++ .github/workflows/validate-tag.yml | 68 +++++++++++++ .github/workflows/validate-zenodo.yml | 35 +++++++ RELEASE.md | 111 +++++++++++++++++++++ 5 files changed, 324 insertions(+) create mode 100644 .github/workflows/pytest.yml create mode 100644 .github/workflows/release-on-published.yml create mode 100644 .github/workflows/validate-tag.yml create mode 100644 .github/workflows/validate-zenodo.yml create mode 100644 RELEASE.md diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml new file mode 100644 index 0000000..521d50a --- /dev/null +++ b/.github/workflows/pytest.yml @@ -0,0 +1,35 @@ +name: Pytest + +on: + push: + branches: [ "main", "dev" ] + pull_request: + branches: [ "main", "dev" ] + +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + python-version: ["3.9", "3.10", "3.11", "3.12"] + + steps: + - uses: actions/checkout@v6 + with: + fetch-depth: 0 + fetch-tags: true + + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install setuptools wheel pytest + if [ -f pyproject.toml ]; then pip install -e . ; fi + + - name: Test with pytest + run: | + pytest pytests/ diff --git a/.github/workflows/release-on-published.yml b/.github/workflows/release-on-published.yml new file mode 100644 index 0000000..1893ac2 --- /dev/null +++ b/.github/workflows/release-on-published.yml @@ -0,0 +1,75 @@ +name: Release On Published + +on: + release: + types: [published] + workflow_dispatch: + inputs: + tag: + description: "Release tag to publish" + required: true + type: string + +permissions: + contents: read + +jobs: + validate-tag: + uses: ./.github/workflows/validate-tag.yml + with: + tag: ${{ github.event.release.tag_name || inputs.tag }} + permissions: + contents: read + + validate-zenodo: + needs: validate-tag + uses: ./.github/workflows/validate-zenodo.yml + with: + ref: ${{ needs.validate-tag.outputs.ref }} + permissions: + contents: read + + build: + runs-on: ubuntu-latest + needs: [validate-tag, validate-zenodo] + + steps: + - name: Checkout + uses: actions/checkout@v6 + with: + ref: ${{ needs.validate-tag.outputs.ref }} + fetch-depth: 0 + fetch-tags: true + + - name: Setup Python + uses: actions/setup-python@v6 + with: + python-version: "3.11" + + - name: Build distributions + run: | + python -m pip install --upgrade pip + python -m pip install build twine + python -m build + python -m twine check dist/* + + - name: Verify distribution version + shell: bash + run: | + expected="${{ needs.validate-tag.outputs.version }}" + test -f "dist/pynetdesign-${expected}.tar.gz" || { + echo "Expected sdist dist/pynetdesign-${expected}.tar.gz was not built" + ls -la dist + exit 1 + } + compgen -G "dist/pynetdesign-${expected}-*.whl" > /dev/null || { + echo "Expected wheel dist/pynetdesign-${expected}-*.whl was not built" + ls -la dist + exit 1 + } + + - name: Upload built artifacts + uses: actions/upload-artifact@v7 + with: + name: dist-release-published + path: dist/ diff --git a/.github/workflows/validate-tag.yml b/.github/workflows/validate-tag.yml new file mode 100644 index 0000000..b7125bb --- /dev/null +++ b/.github/workflows/validate-tag.yml @@ -0,0 +1,68 @@ +name: Validate Tag + +on: + workflow_call: + inputs: + tag: + required: true + type: string + outputs: + tag: + description: "Validated tag name" + value: ${{ jobs.validate.outputs.tag }} + ref: + description: "Validated tag ref" + value: ${{ jobs.validate.outputs.ref }} + version: + description: "Validated version without leading v" + value: ${{ jobs.validate.outputs.version }} + +permissions: + contents: read + +jobs: + validate: + runs-on: ubuntu-latest + outputs: + tag: ${{ steps.tag.outputs.tag }} + ref: ${{ steps.tag.outputs.ref }} + version: ${{ steps.tag.outputs.version }} + + steps: + - name: Checkout tag + uses: actions/checkout@v6 + with: + ref: ${{ inputs.tag }} + fetch-depth: 0 + fetch-tags: true + + - name: Validate strict semver tag + id: tag + shell: bash + run: | + TAG="${{ inputs.tag }}" + if [[ "$TAG" == refs/tags/* ]]; then + TAG="${TAG##*/}" + fi + if [[ ! "$TAG" =~ ^v[0-9]+\.[0-9]+\.[0-9]+$ ]]; then + echo "Invalid tag '$TAG'. Expected strict semver format: vMAJOR.MINOR.PATCH" + exit 1 + fi + echo "tag=$TAG" >> "$GITHUB_OUTPUT" + echo "ref=refs/tags/$TAG" >> "$GITHUB_OUTPUT" + echo "version=${TAG#v}" >> "$GITHUB_OUTPUT" + + - name: Ensure tag points to main history + shell: bash + run: | + TAG="${{ steps.tag.outputs.tag }}" + git fetch origin main + TAG_COMMIT="$(git rev-list -n 1 "$TAG")" + if [[ -z "$TAG_COMMIT" ]]; then + echo "Unable to resolve commit for tag '$TAG'" + exit 1 + fi + if ! git merge-base --is-ancestor "$TAG_COMMIT" origin/main; then + echo "Tag '$TAG' does not point to a commit contained in main" + exit 1 + fi diff --git a/.github/workflows/validate-zenodo.yml b/.github/workflows/validate-zenodo.yml new file mode 100644 index 0000000..db40e37 --- /dev/null +++ b/.github/workflows/validate-zenodo.yml @@ -0,0 +1,35 @@ +name: Validate Zenodo Metadata + +on: + push: + branches: ["main", "dev"] + pull_request: + branches: ["main", "dev"] + workflow_call: + inputs: + ref: + required: false + type: string + workflow_dispatch: + inputs: + ref: + description: "Optional branch, tag, or SHA to validate" + required: false + type: string + +permissions: + contents: read + +jobs: + validate: + runs-on: ubuntu-latest + steps: + - name: Checkout code + uses: actions/checkout@v6 + with: + ref: ${{ inputs.ref || github.ref }} + + - name: Validate Zenodo Metadata + uses: vsoch/zenodo-validator@main + with: + path: ".zenodo.json" diff --git a/RELEASE.md b/RELEASE.md new file mode 100644 index 0000000..d0cac1d --- /dev/null +++ b/RELEASE.md @@ -0,0 +1,111 @@ +# Releasing PyNetDesign + +This project uses [`setuptools_scm`](https://setuptools-scm.readthedocs.io), so the package +version is derived from Git tags. Use strict semantic version tags of the form `vX.Y.Z` +(for example `v1.0.1`). Any commit that is not exactly on a tag resolves to a development +version such as `1.0.1.dev3+g1a2b3c4`. + +## Versioning rules + +- Tags must match `vMAJOR.MINOR.PATCH` exactly. Pre-release and metadata suffixes are + rejected by CI. +- A release tag must point at a commit contained in `main`. Tags created on other branches + are rejected during validation. +- Development happens on `dev`; `dev` is merged into `main` through a pull request, and the + tag is placed on the resulting merge commit on `main`. + +## Standard release flow + +1. Make sure `dev` is green and up to date with `main`. +2. Update `CHANGELOG.md`: move the relevant `## [Unreleased]` entries into a new + `## [vX.Y.Z] - YYYY-MM-DD` section and add the comparison links at the bottom. +3. Update `.zenodo.json`: set `version` to the new release and refresh `description` with + the release highlights. +4. Open and merge the `dev` → `main` pull request. +5. Tag the merge commit on `main` and push the tag: + + ```bash + git checkout main + git pull + git tag -a vX.Y.Z -m "PyNetDesign vX.Y.Z" + git push origin vX.Y.Z + ``` + +6. Verify the version resolves correctly before publishing: + + ```bash + git describe --tags # -> vX.Y.Z + python -m build # -> dist/pynetdesign-X.Y.Z.tar.gz + ``` + +7. Publish the GitHub Release for that tag, using the new `CHANGELOG.md` section as the + release body. +8. Merge `main` back into `dev` so that `dev` contains the tag in its history and + `setuptools_scm` resolves a sensible development version there. + +## Continuous integration + +| Workflow | Trigger | What it does | +| --- | --- | --- | +| `pytest.yml` | push / PR on `main`, `dev` | Runs `pytest pytests/` on Python 3.9-3.12 | +| `validate-zenodo.yml` | push / PR on `main`, `dev`; reusable | Validates `.zenodo.json` against the Zenodo schema | +| `validate-tag.yml` | reusable only | Enforces strict semver and that the tag is contained in `main` | +| `release-on-published.yml` | GitHub Release `published` | Validates the tag, validates the Zenodo metadata, then builds and checks the distributions | + +`release-on-published.yml` also accepts a `workflow_dispatch` with a `tag` input, so a +release can be re-validated without republishing it. + +Notes: + +- Zenodo metadata validation must pass before the package build starts. +- Tag validation is shared through `validate-tag.yml`, so the same rules apply to every + workflow that acts on a tag. +- PyNetDesign is not published to PyPI. Installation is from source; see `README.md`. + +## Zenodo + +The repository is archived on Zenodo through the GitHub integration. + +- GitHub repository ID: `958300487` +- Original manual deposit (v1.0.0): [10.5281/zenodo.14945917](https://doi.org/10.5281/zenodo.14945917) + +### One-time integration setup + +1. Sign in to (test the flow on first if + you prefer). +2. In Zenodo's GitHub settings, enable the `PyNetDesign` repository. +3. Confirm the repository has a valid `.zenodo.json` — `validate-zenodo.yml` checks this on + every push. + +> **Important:** v1.0.0 was deposited manually before the GitHub integration existed. +> The integration was therefore enabled only *after* the `v1.0.0` GitHub release was +> published, so that Zenodo would not mint a second, duplicate DOI for that version. The +> first integration-created deposit is `v1.0.1`, which starts a new concept DOI. + +### Per-release checklist + +1. **Before tagging** — update `.zenodo.json` (`version`, `description`, `keywords`, + `creators`, `related_identifiers` as needed) and make sure `LICENSE` and `README.md` + are consistent with the release notes. +2. **Create the release** — push the `vX.Y.Z` tag, then publish the GitHub Release with + the highlights from `CHANGELOG.md`. +3. **Validate the deposition** — confirm Zenodo auto-created a new versioned deposition + from the GitHub release, and check the title, version, creators, ORCID, license and + description in the Zenodo UI. Verify the source archive is attached. +4. **Communicate the DOIs** — copy the version DOI and the concept DOI into the release + notes, `README.md` and `docs/source/citing.rst`. The concept DOI always resolves to the + latest version; the version DOI pins one release. +5. **Post-release verification** — open the DOI link, check that the citation metadata + resolves, and confirm the new version appears under the concept record. + +### Linking the pre-integration record + +Because the manual v1.0.0 deposit and the integration-created records are separate Zenodo +lineages, they must be cross-linked by hand in both directions: + +- **Forward**, in the repository: `.zenodo.json` carries a `related_identifiers` entry with + relation `isNewVersionOf` pointing at `10.5281/zenodo.14945917`. Zenodo honours this on + every deposit, so it persists across releases. +- **Backward**, on Zenodo: open record `14945917`, choose *Edit*, add a related identifier + with relation `isPreviousVersionOf` pointing at the new version DOI, and publish. This is + a metadata-only edit and preserves the existing DOI. From 51f3db77b5f21ba266939c337e1e1d2fd3c4f1d5 Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Wed, 9 Sep 2026 18:27:34 +0300 Subject: [PATCH 2/4] Update gitignore --- .gitignore | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/.gitignore b/.gitignore index df4c776..6e38af4 100644 --- a/.gitignore +++ b/.gitignore @@ -204,3 +204,11 @@ cython_debug/ # PyPI configuration file .pypirc + + +# Agents +.agent/ +.codex* +AGENTS.md +CLAUDE.md +.claude/ \ No newline at end of file From 7e11b0194415fae20a128548c174503dd7efc7d7 Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Wed, 9 Sep 2026 23:03:25 +0300 Subject: [PATCH 3/4] Fix seismic moment underestimated by a factor of 2*pi*f calculate_M0 applied the amplitude-to-displacement spectral conversion a second time, after retrieve_min_amps had already applied it through get_scaling_displacement. The seismic moment was therefore too small by 2*pi*f_r and all minimum detectable moment magnitudes were too low by dMw = (2/3)*log10(2*pi*f_r), between about 1.0 and 1.9 magnitude units. Because f_r = min(v*Q/(pi*r), f_c) depends on the source-receiver distance the error varies across the grid, so it changed the shape of the sensitivity maps and not only their colour scale. Verified against the reference implementation on test case 1 with the free surface correction disabled in both: over all 81000 grid points the corrected output now reproduces the reference to within 1e-15, up to the known constant 2/3*log10(0.52/sqrt(4/15)) = 0.0020 that comes from the different radiation pattern constants. Fixes #3 - drop the spurious 1/(2*pi*f) factor in calculate_M0 and document the relation M0 = 4*pi*rho*v^3*r/|R| * |U(f)| * exp(pi*f*t*) that it evaluates, stating that amps is the displacement amplitude spectrum and giving units for the other arguments - add pytests/test_magnitude.py covering the seismic moment relation, the moment magnitude conversion and its inverse, the amplitude-type scaling coefficients, the peak frequency and its corner clamp, and the DAS strain and strain-rate thresholds; the closed-form checks fail on the previous formula - remove loc_uncertainty_grid, which called calculate_pdf and get_uncertainty_from_pdf; neither has ever been part of the package, so every call raised NameError - replace the match statement in plot_sensitivity_slice with if/elif so the visualization module imports on Python 3.9, which the CI matrix covers - record the release in CHANGELOG.md and .zenodo.json, and link the new Zenodo deposit to the manual v1.0.0 record with isNewVersionOf --- .zenodo.json | 11 +- CHANGELOG.md | 24 +- pynetdesign/modelling/homo.py | 287 ++--------------------- pynetdesign/visualization/sensitivity.py | 19 +- pytests/test_magnitude.py | 226 ++++++++++++++++++ 5 files changed, 283 insertions(+), 284 deletions(-) create mode 100644 pytests/test_magnitude.py diff --git a/.zenodo.json b/.zenodo.json index dd0725a..3952c93 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,9 +1,9 @@ { "title": "PyNetDesign: a Python framework for microseismic monitoring network design", - "version": "1.0.0", + "version": "1.0.1", "license": "MIT", "upload_type": "software", - "description": "PyNetDesign is an open-source Python library that provides a simple and functional framework for designing and testing microseismic monitoring networks.", + "description": "PyNetDesign is an open-source Python library that provides a simple and functional framework for designing and testing microseismic monitoring networks. Version 1.0.1 is a bugfix release. It corrects the seismic moment calculation, which underestimated the moment by a factor of 2*pi*f because the amplitude-to-displacement spectral conversion was applied twice. Minimum detectable moment magnitudes reported by version 1.0.0 were too low by roughly 1.0 to 1.9 magnitude units, by an amount that varies with the source-receiver distance, so results produced with version 1.0.0 should be recomputed. The release also restores Python 3.9 compatibility of the visualization module, removes a non-functional location uncertainty routine that referenced functions absent from the package, and adds a test module covering the magnitude relations.", "keywords": [ "earthquakes", "microseismicity", @@ -18,6 +18,13 @@ ], "access_right": "open", "language": "eng", + "related_identifiers": [ + { + "identifier": "10.5281/zenodo.14945917", + "relation": "isNewVersionOf", + "resource_type": "software" + } + ], "creators": [ { "affiliation": "Center for Integrative Petroleum Research, King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia", diff --git a/CHANGELOG.md b/CHANGELOG.md index 2961de5..475e703 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [v1.0.1] - 2026-09-09 + +Bugfix release correcting the seismic moment calculation. + +### Fixed + +- the seismic moment is no longer underestimated by a factor of `2*pi*f` (#3). `calculate_M0` applied the amplitude-to-displacement spectral conversion a second time, after `retrieve_min_amps` had already applied it through `get_scaling_displacement`. Minimum detectable moment magnitudes reported by v1.0.0 were therefore too low by `dMw = (2/3)*log10(2*pi*f_r)`, which is 1.20 at a representative frequency of 10 Hz and 1.87 at 100 Hz. Because the representative frequency `f_r = min(v*Q/(pi*r), f_c)` depends on the source-receiver distance, the correction varies across the grid and changes the shape of the sensitivity maps, not only their colour scale +- `pynetdesign.visualization` can be imported on Python 3.9 again; `plot_sensitivity_slice` used a `match` statement, which is syntax available only from Python 3.10 onwards + +### Changed + +- `calculate_M0` now documents the seismic moment relation it evaluates and states that the amplitudes it takes are the displacement amplitude spectrum, with units given for the remaining arguments + +### Removed + +- `loc_uncertainty_grid`, which called `calculate_pdf` and `get_uncertainty_from_pdf`. Neither function has ever been part of the package, so every call raised `NameError`. Location uncertainty will be reintroduced later together with the functions it depends on + +### Added + +- `pytests/test_magnitude.py`, covering the seismic moment relation, the moment magnitude conversion and its inverse, the amplitude-type scaling coefficients, the peak frequency and its corner-frequency clamp, and the DAS strain and strain-rate detection thresholds. The closed-form checks fail on the v1.0.0 formula and so guard against a recurrence of #3 + ## [v1.0.0] - 2025-04-01 Initial public release, archived on Zenodo as @@ -35,5 +56,6 @@ code itself is unchanged from the deposited version. - conda environment files for user and development installations, and Windows installation scripts - MIT license -[Unreleased]: https://github.com/danikiev/PyNetDesign/compare/v1.0.0...HEAD +[Unreleased]: https://github.com/danikiev/PyNetDesign/compare/v1.0.1...HEAD +[v1.0.1]: https://github.com/danikiev/PyNetDesign/compare/v1.0.0...v1.0.1 [v1.0.0]: https://github.com/danikiev/PyNetDesign/releases/tag/v1.0.0 diff --git a/pynetdesign/modelling/homo.py b/pynetdesign/modelling/homo.py index 34bb4ea..a1a377c 100644 --- a/pynetdesign/modelling/homo.py +++ b/pynetdesign/modelling/homo.py @@ -243,262 +243,6 @@ def calculate_Mw_min(f, fcorner, density, v, Q, rad_pattern, min_amps, fs_coef, return Mw_min_grid -def loc_uncertainty_grid(grid_coords: np.ndarray, - velocity_df: pd.DataFrame, - loc_gx: np.ndarray, - loc_gy: np.ndarray, - loc_gz: np.ndarray, - source_mw: float, - geometry_df: pd.DataFrame = None, - station_coords: np.ndarray = None, - min_amps_p: np.ndarray = None, - min_amps_s: np.ndarray = None, - rad_pattern_p: float = 0.52, - rad_pattern_s: float = 0.63, - sigma_p: float = None, - sigma_s: float = None, - min_sigma: float = None, - wave_mode: str = 'PS', - use_free_surface: bool = False, - use_station_directionality: bool = False, - output_pdf: bool = False, - strict_nan_check: bool = False): - - r""" - Computes location uncertainty of a microseismic event using P and/or S waves - assuming 3D elastic homogeneous medium with attenuation. - It takes into account picking uncertainties for P and/or S waves. - - Parameters - ---------- - grid_coords : :obj:`numpy.ndarray` - Array of imaging grid coordinates [[xs1, ys1, zs1], [xs2, ys2, zs2], ...] - velocity_df : :obj:`pandas.DataFrame` - Velocity model DataFrame - loc_gx : :obj:`numpy.ndarray` - Vector of x-coordinates of the location grid - loc_gy : :obj:`numpy.ndarray` - Vector of y-coordinates of the location grid - loc_gz : :obj:`numpy.ndarray` - Vector of z-coordinates of the location grid - source_mw : :obj:`float` - Moment magnitude of the synthetic source event - geometry_df : :obj:`pandas.DataFrame`, optional - Geometry DataFrame, use if station_coords is None - station_coords : :obj:`numpy.ndarray`, optional - Array of station coordinates [[xr1, yr1, zr1], [xr2, yr2, zr2], ...], use if geometry_df is None - min_amps_p : :obj:`numpy.ndarray`, optional - Minimal measurable displacement amplitudes of on stations [a1, a2, ...] for P waves - min_amps_s : :obj:`numpy.ndarray`, optional - Minimal measurable displacement amplitudes on stations [a1, a2, ...] for S waves - rad_pattern_p : :obj:`float`, optional, default: 0.52 - Radiation pattern factor for P waves - rad_pattern_s : :obj:`float`, optional, default: 0.63 - Radiation pattern factor for S waves - sigma_p : :obj:`float`, optional - Picking uncertainty for P waves, in seconds. - If ``None``, estimated from the peak frequency of the P wave. - sigma_s : :obj:`float`, optional - Picking uncertainty for S waves, in seconds. - If ``None``, estimated from the peak frequency of the S wave. - min_sigma : :obj:`float`, optional, default: 1e-3 (when None) - Minimum picking uncertainty, in seconds. - wave_mode: :obj:`str`, optional, default: 'PS' - Wave mode to use, can be 'P','S' or 'PS' - use_free_surface : :obj:`bool`, optional, default: False - Use free surface correction. - If True, this will consider amplification of the amplitudes when recording on the daily surface - due to the free surface boundary condition. - use_station_directionality : :obj:`bool`, optional, default: False - Use station directionality correction. - This will consider sensitivity of the measurement only in direction parallel to the adjacent stations - output_pdf : :obj:`bool`, optional, default: False - Return probability density function - strict_nan_check: :obj:`bool`, optional, default: False - flag to check if NaN values appear in the result - - Returns - ------- - LU_grid : :obj:`numpy.ndarray` - Array of computed location uncertainties (azimuth, sigma_a, sigma_b, sigma_z) on the grid - pdf : :obj:`numpy.ndarray` - 3D probability density function of the event location, if ``output_pdf = True`` - - Raises - ------ - ValueError : - if both station_coords and geometry_df are not None or None - ValueError : - if there is inconsistency in P/S wave parameters - ValueError : - if number of minimal measurable amplitudes do not match number of stations - ValueError : - if density, wave velocity or attenuation factor is not positive - ValueError : - if radiation pattern factor for P or S wave is less or equal to 0 or higher than 1 - ValueError : - if number of minimal measurable amplitudes does not match number of stations - RuntimeError : - if strict_nan_check is True and NaN appears in the result - - Notes - ----- - If ``station_coords`` is not provided, ``geometry_df`` must be provided and vice versa. - If ``output_pdf`` is True, the function will return a tuple of two elements: LU_grid and pdf. - """ - - if (station_coords is not None and geometry_df is not None) or (station_coords is None and geometry_df is None): - raise ValueError("You must provide either `station_coords` or `geometry_df`, but not both.") - - # Get station coordinates from geometry - if station_coords is None: - station_coords = geometry_df[['X', 'Y', 'Z']].to_numpy() - - # Handle the case if station_coords has only one point (3 elements) - if station_coords.ndim == 1 and station_coords.size == 3: - station_coords = station_coords[np.newaxis, :] # reshape to (1, 3) - - # Handle the case if grid_coords has only one point (3 elements) - if grid_coords.ndim == 1 and grid_coords.size == 3: - grid_coords = grid_coords[np.newaxis, :] # reshape to (1, 3) - - # Check velocity model - if len(velocity_df) != 1: - raise ValueError("Velocity model is not homogeneous!") - - # Retrieve medium parameters - density = velocity_df.at[0,'Rho'] - v_p = velocity_df.at[0,'Vp'] - Q_p = velocity_df.at[0,'Qp'] - v_s = velocity_df.at[0,'Vp']/velocity_df.at[0,'VpVsRatio'] - Q_s = velocity_df.at[0,'Qs'] - - # Validate parameters - validate_parameters(density=density, - v_p=v_p, - v_s=v_s, - Q_s=Q_s, - Q_p=Q_p, - min_amps_p=min_amps_p, - min_amps_s=min_amps_s, - rad_pattern_p=rad_pattern_p, - rad_pattern_s=rad_pattern_s, - wave_mode=wave_mode, - station_coords=station_coords) - - # Get minimal seismic moment - if source_mw is None: - raise ValueError("Source seismic moment can't be None.") - else: - sourceSeismicMoment = calculate_seismic_moment(Mw=source_mw) - - # Check minimal picking uncertainty - if min_sigma is None: - min_sigma = 1e-3 - - # Generate location grid - loc_grid_coords = generate_grid(x=loc_gx, y=loc_gy, z=loc_gz) - - # Precompute distances between all imaging grid points and all stations - # Shape: (n_grid_points, n_stations) - img_distances = np.linalg.norm(grid_coords[:, np.newaxis] - station_coords, axis=2) - - # Precompute distances between all location grid points and all stations - # Shape: (n_loc_grid_points, n_stations) - loc_distances = np.linalg.norm(loc_grid_coords[:, np.newaxis] - station_coords, axis=2) - - # Check for NaNs - if strict_nan_check: - if np.any(img_distances == 0): - raise ValueError("Zero imaging distance: source and receiver must have different coordinates") - if np.any(loc_distances == 0): - raise ValueError("Zero location distance: source and receiver must have different coordinates") - - # Calculate free surface correction coefficient if needed - if use_free_surface: - # Simple approximation - fs_coef_p = fs_coef_s = 2*np.ones_like(img_distances) - else: - # If not using free surface correction, use ones (no effect on calculations) - fs_coef_p = fs_coef_s = np.ones_like(img_distances) - - # Calculate station directionality if needed - if use_station_directionality: - stations_dir_coef_p, stations_dir_coef_s = get_ray_station_directionality(station_coords=station_coords, - grid_coords=grid_coords, - strict_nan_check=strict_nan_check) - else: - # If not using directionality, use ones (no effect on calculations) - stations_dir_coef_p = stations_dir_coef_s = np.ones_like(img_distances) - - # Compute traveltimes and pick uncertainties - if wave_mode in ('P', 'PS'): - tpl = loc_distances / v_p # travel times for P waves from location grid to receivers - tpi = img_distances / v_p # travel times for P waves from imaging grid to receivers - - # Get picking uncertainty for P waves from the peak frequency - if sigma_p is None: - sigma_p = 1.0/calculate_fpeak(v=v_p, Q=Q_p, r=img_distances) - sigma_p = np.maximum(sigma_p, min_sigma) - - # Estimate minimal detectable seismic moment for P wave - minSeismicMomentP = calculate_M0(density=density, - v=v_p, - Q=Q_p, - rad_pattern=rad_pattern_p, - r=img_distances, - amps = min_amps_p * stations_dir_coef_p / fs_coef_p) - - # Remove P wave arrivals if they are under detection threshold - tpi[minSeismicMomentP > sourceSeismicMoment] = np.nan - else: - tpl = None - tpi = None - - if wave_mode in ('S', 'PS'): - tsl = loc_distances / v_s # travel times for S waves from location grid to receivers - tsi = img_distances / v_s # travel times for S waves from imaging grid to receivers - # Get picking uncertainty for S waves from the peak frequency - if sigma_s is None: - sigma_s = 1.0/calculate_fpeak(v=v_s, Q=Q_s, r=img_distances) - sigma_s = np.maximum(sigma_s, min_sigma) - # Estimate minimal detectable seismic moment for P wave - minSeismicMomentS = calculate_M0(density=density, - v=v_s, - Q=Q_s, - rad_pattern=rad_pattern_s, - r=img_distances, - amps = min_amps_s * stations_dir_coef_s / fs_coef_s) - # Remove S wave arrivals if they are under detection threshold - tsi[minSeismicMomentS > sourceSeismicMoment] = np.nan - else: - tsl = None - tsi = None - - # Calculate the 3D probability density function - pdf = calculate_pdf(sigma_p=sigma_p, - sigma_s=sigma_s, - tpl=tpl, - tsl=tsl, - tps=tpi, - tss=tsi) - - # Get location uncertainty - azimuth_a, sigma_a, sigma_b, sigma_z = get_uncertainty_from_pdf( - pdf, - loc_gx=loc_gx, - loc_gy=loc_gy, - loc_gz=loc_gz) - - # Combine uncertainties to a single array - LU_grid = np.array([azimuth_a, sigma_a, sigma_b, sigma_z]) - - # Output depending on the request - if output_pdf: - return LU_grid, pdf - else: - return LU_grid - def calculate_fpeak(v: float, Q: float, r: float, @@ -561,17 +305,18 @@ def calculate_M0(density: float, fcorner: float = None, amps_type: str='displacement'): r""" - Calculates the seismic moment :math:`M_0` given the density :math:`\rho`, wave velocity :math:`v` and distance from the source :math:`r` using: + Calculates the seismic moment :math:`M_0` in a homogeneous medium given the density :math:`\rho`, wave velocity :math:`v`, distance from the source :math:`r` and the displacement amplitude spectrum :math:`\left|U(f)\right|` using: .. math:: - M_0 = \frac{4 \pi \rho v^3 r \Omega_0}{R}, + M_0 = \frac{4 \pi \rho v^3 r}{\left|R\right|}\,\left|U(f)\right|\,e^{\pi f t^*}, - where :math:`R` is the provided radiation pattern factor and the source spectra :math:`\Omega_0` given the frequency :math:`f` and displacement amplitude :math:`A` is calculated as - - .. math:: - \Omega_0 = \frac{A}{2 \pi f e^{-\pi f t^*}} + where :math:`R` is the provided radiation pattern factor and the exponential term corrects + for the intrinsic attenuation accumulated along the ray path. - coming from a zero-frequency limit approximation. + The amplitudes passed in ``amps`` are the displacement amplitude spectrum + :math:`\left|U(f)\right|`. Amplitudes of another type are converted to it with the + coefficients of :func:`~pynetdesign.modelling.utils.get_scaling_displacement`, which + follow from :math:`\left|V(f)\right| = 2 \pi f \left|U(f)\right|`. The term :math:`t^*` is the integral along the ray path of the inverse value of attenuation factor :math:`Q` multiplied by the inverse of the wave velocity :math:`v`: @@ -590,17 +335,18 @@ def calculate_M0(density: float, Parameters ---------- density : :obj:`float` - Fensity + Density (kg/m^3) v : :obj:`float` - Wave velocity + Wave velocity (m/s) Q : :obj:`float` Attenuation factor rad_pattern : :obj:`float` Radiation pattern factor r : :obj:`float` or :obj:`numpy.ndarray` - Distance(s) from the source + Distance(s) from the source (m) amps : :obj:`float` or :obj:`numpy.ndarray` - Displacement amplitudes, single value or array of amplitudes, one for each ray. + Displacement amplitude spectrum :math:`\left|U(f)\right|`, single value or array of + amplitudes, one for each ray. If ``r`` is a single float, ``amps`` must be a single float If ``r`` is an array, ``amps`` must be an array of the same shape. f : :obj:`float`, optional, default: ``None`` @@ -614,7 +360,7 @@ def calculate_M0(density: float, Returns ------- M0 : :obj:`float` or array_like - Seismic moment + Seismic moment (N m) """ # Check validity of 'r' and 'amps' if amps is None: @@ -661,9 +407,8 @@ def calculate_M0(density: float, # Compute exp(pi * f * t*) exp_term = np.exp(pif * t_star) - # Compute the source spectra (Omega_0) and seismic moment (M0) - Omega_0 = scaling * amps * exp_term / (2 * pif) - M0 = constant * r * Omega_0 + # Compute seismic moment: M0 = (4 pi rho v^3 r / |R|) * |U(f)| * exp(pi f t*) + M0 = constant * r * scaling * amps * exp_term return M0 diff --git a/pynetdesign/visualization/sensitivity.py b/pynetdesign/visualization/sensitivity.py index fbec601..c9a9386 100644 --- a/pynetdesign/visualization/sensitivity.py +++ b/pynetdesign/visualization/sensitivity.py @@ -165,16 +165,15 @@ def plot_sensitivity_slice(sens_data: np.ndarray, levels = plt_levels # Plot the sensitivity data as a heatmap in desired style - match plt_style: - case "contourf": - im = ax.contourf(ci1,ci2,pdata, levels, cmap=cmap, extent=extent) - case "imshow": - im = ax.imshow(pdata, - origin='lower', - extent=extent, - cmap=cmap) - case _: - raise ValueError("Unknown plot style: " + plt_style) + if plt_style == "contourf": + im = ax.contourf(ci1,ci2,pdata, levels, cmap=cmap, extent=extent) + elif plt_style == "imshow": + im = ax.imshow(pdata, + origin='lower', + extent=extent, + cmap=cmap) + else: + raise ValueError("Unknown plot style: " + plt_style) # Plot contours if plt_contours: diff --git a/pytests/test_magnitude.py b/pytests/test_magnitude.py new file mode 100644 index 0000000..242346c --- /dev/null +++ b/pytests/test_magnitude.py @@ -0,0 +1,226 @@ +"""Tests for the magnitude / seismic-moment relations. + +The closed-form checks in this module guard against the regression reported in +issue #3, where ``calculate_M0`` applied the amplitude-to-displacement spectral +conversion twice and therefore underestimated the seismic moment by a factor of +``2 * pi * f``. +""" +import numpy as np +import pandas as pd +import pytest + +from pynetdesign.modelling.homo import calculate_M0, calculate_fpeak +from pynetdesign.modelling.utils import ( + calculate_Mw, + calculate_seismic_moment, + get_scaling_displacement, + retrieve_min_amps, +) + +# Medium of test case 1 (Hallo, 2012) +RHO = 2300.0 +VP = 4300.0 +QP = 100.0 +RAD_P = 0.52 + + +def m0_closed_form(density, v, Q, rad_pattern, r, u, f): + """M0 = (4 pi rho v^3 r / |R|) * |U(f)| * exp(pi f t*), with t* = r / (v Q). + + Written out independently of the package so that it cannot drift with it. + """ + t_star = r / (v * Q) + return 4.0 * np.pi * density * v**3 * r / rad_pattern * u * np.exp(np.pi * f * t_star) + + +class TestCalculateM0: + def test_matches_closed_form_explicit_frequency(self): + """The core relation, with the representative frequency supplied.""" + r, u, f = 1000.0, 1e-9, 10.0 + expected = m0_closed_form(RHO, VP, QP, RAD_P, r, u, f) + actual = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=f, amps_type='displacement') + assert actual == pytest.approx(expected, rel=1e-12) + + def test_no_spurious_two_pi_f_factor(self): + """Regression guard for issue #3. + + The buggy implementation returned ``expected / (2 * pi * f)``. Assert the + ratio is one and, explicitly, that it is not ``2 * pi * f``. + """ + r, u, f = 1000.0, 1e-9, 10.0 + expected = m0_closed_form(RHO, VP, QP, RAD_P, r, u, f) + actual = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=f, amps_type='displacement') + ratio = expected / actual + assert ratio == pytest.approx(1.0, rel=1e-12) + assert ratio != pytest.approx(2.0 * np.pi * f, rel=1e-6) + + def test_peak_frequency_path_below_corner(self): + """With f=None and f_peak below the corner, pi*f_peak*t* == 1 exactly, + so the attenuation term collapses to Euler's number.""" + r, u = 5000.0, 1e-9 # f_peak = vQ/(pi r) = 27.4 Hz < 100 Hz + f_peak = VP * QP / (np.pi * r) + assert f_peak < 100.0 + expected = 4.0 * np.pi * RHO * VP**3 * r / RAD_P * u * np.e + actual = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=None, amps_type='displacement') + assert actual == pytest.approx(expected, rel=1e-12) + + def test_peak_frequency_path_clamped_at_corner(self): + """Close to the receiver the peak frequency is clamped to the corner.""" + r, u, fcorner = 500.0, 1e-9, 100.0 + assert VP * QP / (np.pi * r) > fcorner + expected = m0_closed_form(RHO, VP, QP, RAD_P, r, u, fcorner) + actual = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=None, fcorner=fcorner, + amps_type='displacement') + assert actual == pytest.approx(expected, rel=1e-12) + + def test_explicit_and_automatic_paths_agree(self): + """Supplying f explicitly must match the automatic path evaluated at the + same frequency. The two paths scale the amplitudes in different places, + so this pins them together.""" + r = 5000.0 + f_peak = VP * QP / (np.pi * r) + noise_velocity = 5e-9 + # automatic path: raw velocity amplitude, converted inside calculate_M0 + auto = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=noise_velocity, f=None, + amps_type='velocity') + # explicit path: amplitude pre-converted to a displacement spectrum + u = noise_velocity * get_scaling_displacement(amps_type='velocity', f=f_peak) + explicit = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=f_peak, amps_type='displacement') + assert auto == pytest.approx(explicit, rel=1e-12) + + def test_scales_linearly_with_distance_and_amplitude(self): + r, u, f = 1000.0, 1e-9, 10.0 + base = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=f, amps_type='displacement') + doubled = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=2 * u, f=f, amps_type='displacement') + assert doubled == pytest.approx(2.0 * base, rel=1e-12) + + def test_inversely_proportional_to_radiation_pattern(self): + r, u, f = 1000.0, 1e-9, 10.0 + a = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=0.52, + r=r, amps=u, f=f, amps_type='displacement') + b = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=0.63, + r=r, amps=u, f=f, amps_type='displacement') + assert a / b == pytest.approx(0.63 / 0.52, rel=1e-12) + + def test_array_input(self): + r = np.array([1000.0, 2000.0, 4000.0]) + u = np.full_like(r, 1e-9) + f = 10.0 + expected = m0_closed_form(RHO, VP, QP, RAD_P, r, u, f) + actual = calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=r, amps=u, f=f, amps_type='displacement') + assert actual.shape == r.shape + np.testing.assert_allclose(actual, expected, rtol=1e-12) + + def test_rejects_mismatched_shapes(self): + with pytest.raises(ValueError): + calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=np.array([1000.0, 2000.0]), amps=1e-9, f=10.0) + with pytest.raises(ValueError): + calculate_M0(density=RHO, v=VP, Q=QP, rad_pattern=RAD_P, + r=1000.0, amps=None, f=10.0) + + +class TestMomentMagnitude: + def test_kanamori_relation(self): + """Mw = (2/3) (log10 M0 - 9.1).""" + assert calculate_Mw(10**9.1) == pytest.approx(0.0, abs=1e-12) + assert calculate_Mw(10**(9.1 + 1.5)) == pytest.approx(1.0, rel=1e-12) + assert calculate_Mw(10**(9.1 - 1.5)) == pytest.approx(-1.0, rel=1e-12) + + @pytest.mark.parametrize("mw", [-2.0, -1.0, -0.5, 0.0, 1.0, 2.5, 5.0]) + def test_roundtrip(self, mw): + assert calculate_Mw(calculate_seismic_moment(mw)) == pytest.approx(mw, rel=1e-12) + + def test_non_positive_moment_gives_nan(self): + out = calculate_Mw(np.array([-1.0, 0.0, 10**9.1])) + assert np.isnan(out[0]) and np.isnan(out[1]) + assert out[2] == pytest.approx(0.0, abs=1e-12) + + def test_moment_is_monotonic_in_magnitude(self): + mws = np.array([-1.0, 0.0, 1.0, 2.0]) + assert np.all(np.diff(calculate_seismic_moment(mws)) > 0) + + +class TestScalingToDisplacement: + @pytest.mark.parametrize("amps_type,expected", [ + ('displacement', 1.0), + ('velocity', 1.0 / (2 * np.pi * 10.0)), + ('acceleration', 1.0 / (2 * np.pi * 10.0) ** 2), + ('strain', 1.0), + ('strain rate', 1.0 / (2 * np.pi * 10.0)), + ]) + def test_values(self, amps_type, expected): + assert get_scaling_displacement(amps_type=amps_type, f=10.0) == pytest.approx( + expected, rel=1e-12) + + def test_rejects_unknown_type(self): + with pytest.raises(ValueError): + get_scaling_displacement(amps_type='nonsense', f=10.0) + + +class TestPeakFrequency: + def test_formula_below_corner(self): + r = 5000.0 + assert calculate_fpeak(v=VP, Q=QP, r=r) == pytest.approx( + VP * QP / (np.pi * r), rel=1e-12) + + def test_clamped_to_corner(self): + assert calculate_fpeak(v=VP, Q=QP, r=100.0) == pytest.approx(100.0, rel=1e-12) + assert calculate_fpeak(v=VP, Q=QP, r=100.0, fcorner=50.0) == pytest.approx( + 50.0, rel=1e-12) + + def test_decreases_with_distance(self): + r = np.array([2000.0, 4000.0, 8000.0]) + assert np.all(np.diff(calculate_fpeak(v=VP, Q=QP, r=r)) < 0) + + +class TestDasThreshold: + """End-to-end noise -> displacement-spectrum threshold for a DAS channel, + locking in S(f) = L_g / (2 pi f) for strain rate.""" + + @staticmethod + def das_geometry(noise_type, noise_level, gauge_length): + df = pd.DataFrame({'Name': ['C1'], 'X': [0.0], 'Y': [0.0], 'Z': [0.0], + 'NoiseLevel': [noise_level]}) + df.attrs['noise_type'] = noise_type + df.attrs['gauge_length'] = gauge_length + return df + + def test_strain_rate_threshold(self): + a_n, l_g, snr, f = 1.95e-9, 10.0, 2.0, 10.0 + df = self.das_geometry('strain rate', a_n, l_g) + min_p, min_s = retrieve_min_amps(geometry_df=df, SNR_p=snr, SNR_s=snr, + f_p=f, f_s=f, sample_interval=l_g) + expected = snr * a_n * l_g / (2 * np.pi * f) + assert min_p[0] == pytest.approx(expected, rel=1e-12) + assert min_s[0] == pytest.approx(expected, rel=1e-12) + + def test_strain_threshold(self): + a_n, l_g, snr, f = 5e-9, 10.0, 2.0, 10.0 + df = self.das_geometry('strain', a_n, l_g) + min_p, _ = retrieve_min_amps(geometry_df=df, SNR_p=snr, SNR_s=snr, + f_p=f, f_s=f, sample_interval=l_g) + assert min_p[0] == pytest.approx(snr * a_n * l_g, rel=1e-12) + + def test_requires_gauge_length_for_strain(self): + df = self.das_geometry('strain rate', 1.95e-9, 10.0) + with pytest.raises(ValueError): + retrieve_min_amps(geometry_df=df, SNR_p=2.0, SNR_s=2.0, + f_p=10.0, f_s=10.0, sample_interval=None) + + def test_threshold_is_linear_in_snr(self): + df = self.das_geometry('strain rate', 1.95e-9, 10.0) + a, _ = retrieve_min_amps(geometry_df=df, SNR_p=2.0, SNR_s=2.0, + f_p=10.0, f_s=10.0, sample_interval=10.0) + b, _ = retrieve_min_amps(geometry_df=df, SNR_p=4.0, SNR_s=4.0, + f_p=10.0, f_s=10.0, sample_interval=10.0) + assert b[0] == pytest.approx(2.0 * a[0], rel=1e-12) From 00527437a9d397ea82224e3f29ca62874eb7ed44 Mon Sep 17 00:00:00 2001 From: Denis Anikiev Date: Wed, 9 Sep 2026 23:13:47 +0300 Subject: [PATCH 4/4] Update README to reflect DOI changes and citation guidelines for version 1.0.0 --- README.md | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index f7cdfdb..8fafe1e 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # PyNetDesign -[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14945917.svg)](https://doi.org/10.5281/zenodo.14945917) +[![DOI](https://zenodo.org/badge/958300487.svg)](https://zenodo.org/badge/latestdoi/958300487) +[![Pytest](https://github.com/danikiev/PyNetDesign/actions/workflows/pytest.yml/badge.svg)](https://github.com/danikiev/PyNetDesign/actions/workflows/pytest.yml) [![License: MIT](https://img.shields.io/badge/License-MIT-yellowgreen.svg)](https://github.com/danikiev/PyNetDesign/blob/main/LICENSE) PyNetDesign is an open-source Python framework for testing and designing microseismic monitoring networks for homogeneous velocity models. @@ -197,7 +198,12 @@ All notable changes are documented in [CHANGELOG.md](CHANGELOG.md), following ## Citing -If you use PyNetDesign in your research, please cite it as follows: +If you use PyNetDesign in your research, please cite the release you used. The DOI badge +at the top of this page always resolves to the most recently archived release, and the +Zenodo record of each release also carries its own version DOI. + +Version 1.0.0 was deposited manually, before the automated GitHub–Zenodo integration was +enabled, so it has a DOI of its own: **Anikiev, D. (2025).** *PyNetDesign: a Python framework for microseismic monitoring network design* (Version v1.0.0) [Computer software]. Zenodo. @@ -217,6 +223,11 @@ In BibTeX format: } ``` +**Note:** magnitudes computed with version 1.0.0 are affected by an error in the seismic +moment calculation that was corrected in version 1.0.1. If you are citing or reusing +results obtained with 1.0.0, please see [CHANGELOG.md](CHANGELOG.md) and issue +[#3](https://github.com/danikiev/PyNetDesign/issues/3). + ## License PyNetDesign is released under the MIT License. See [LICENSE](LICENSE) for details.