refactor: use np.einsum for spherical harmonic summations#173
Open
tsutterley wants to merge 4 commits into
Open
refactor: use np.einsum for spherical harmonic summations#173tsutterley wants to merge 4 commits into
np.einsum for spherical harmonic summations#173tsutterley wants to merge 4 commits into
Conversation
refactor: use `np.radians` to convert from degrees to radians refactor: use `np.hypot` to calculate the sum of two squares feat: add dunder (magic) methods for mathematical operations
There was a problem hiding this comment.
Pull request overview
This PR refactors multiple spherical-harmonic and sea level equation code paths to use complex/Euler forms and np.einsum-based contractions, while also standardizing many angle and norm computations (np.radians/np.degrees, np.hypot) and adding arithmetic dunder methods on core data containers.
Changes:
- Refactor spherical harmonic summations/transforms (and related geocenter / sea-level computations) toward
np.einsum+ complex representations. - Standardize degree↔radian conversions and “sum of squares” computations across scripts, library code, and docs/tests.
- Add operator overloads (dunder methods) on
spatialandharmonicsobjects and introduce a sea level equation regression test.
Reviewed changes
Copilot reviewed 74 out of 76 changed files in this pull request and generated 13 comments.
Show a summary per file
| File | Description |
|---|---|
| utilities/quick_mascon_regress.py | Use np.hypot/np.degrees for amplitude/phase error calculations. |
| test/test_sea_level.py | New regression test for sea_level_equation output against reference grid. |
| test/test_masks.py | Replace manual radian conversions with np.radians. |
| test/test_harmonics.py | Update theta conversion and avoid mutating shared harmonics in transform test. |
| test/test_download_and_read.py | Add pytest import; skip deprecated GFZ FTP test; update ITSG host path. |
| scripts/sea_level_stokes.py | Use np.radians for colatitude conversion. |
| scripts/sea_level_regress.py | Use np.hypot/np.degrees; refactor phase-wrapping logic. |
| scripts/sea_level_error.py | Use np.radians for colatitude conversion. |
| scripts/scale_grace_maps.py | Use np.radians and minor formatting tweaks. |
| scripts/run_sea_level_equation.py | Use np.radians for colatitude conversion. |
| scripts/regress_grace_maps.py | Use np.where phase wrap; use np.hypot/np.degrees. |
| scripts/regional_spherical_caps.py | Use np.radians/np.degrees for angle conversions. |
| scripts/plot_SLR_azimuthal.py | Make LaTeX labels raw strings. |
| scripts/piecewise_grace_maps.py | Use np.flatnonzero; np.where phase wrap; np.hypot/np.degrees. |
| scripts/grace_spatial_maps.py | Use np.radians for colatitude conversion. |
| scripts/grace_spatial_error.py | Use np.radians and minor formatting tweaks. |
| scripts/convert_harmonics.py | Use np.radians for colatitude conversion. |
| scripts/combine_sea_level_data.py | Use np.radians for colatitude conversion. |
| scripts/combine_HEX_spherical_caps.py | Use np.hypot for combined error. |
| scripts/combine_harmonics.py | Use np.radians for colatitude conversion. |
| scripts/calc_sensitivity_kernel.py | Use np.radians for colatitude conversion. |
| mapping/plot_QML_grid_3maps.py | Use np.radians for cell-area computations and cosine latitude term. |
| mapping/plot_GrIS_grid_movie.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_GrIS_grid_maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_GrIS_grid_5maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_GrIS_grid_3maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_movie.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_9maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_5maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_4maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_global_grid_3maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_regional_movie.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_regional_maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_GrIS_maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_grid_movie.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_grid_maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| mapping/plot_AIS_grid_4maps.py | Same np.radians area/cos-lat refactor; minor formatting cleanup. |
| mapping/plot_AIS_grid_3maps.py | Same np.radians area/cos-lat refactor as other mapping scripts. |
| gravity_toolkit/tools.py | Use np.radians/np.degrees in wrap_longitudes. |
| gravity_toolkit/time_series/smooth.py | Use np.hypot/np.degrees for amplitude/phase. |
| gravity_toolkit/time_series/piecewise.py | Use np.hypot when combining error terms. |
| gravity_toolkit/time_series/amplitude.py | Use np.degrees(np.arctan2(...)) for phase. |
| gravity_toolkit/spatial.py | Add to_ascii(float_format=...) and add arithmetic dunder methods. |
| gravity_toolkit/SLR/C50.py | Unpack scalar from convert_calendar_decimal output array. |
| gravity_toolkit/SLR/C30.py | Unpack scalar from convert_calendar_decimal output array. |
| gravity_toolkit/SLR/C20.py | Unpack scalar from convert_calendar_decimal output array. |
| gravity_toolkit/sea_level_equation.py | Major refactor: complex coefficients, np.einsum summations, _clenshaw helper. |
| gravity_toolkit/mascons.py | Use np.where, np.flatnonzero, and np.radians for refactors/robustness. |
| gravity_toolkit/harmonics.py | Add arithmetic dunder methods. |
| gravity_toolkit/harmonic_summation.py | Refactor summations to np.einsum + Euler complex form. |
| gravity_toolkit/harmonic_gradients.py | Refactor gradients/currents computations to np.einsum + Euler complex form. |
| gravity_toolkit/grace_input_months.py | Use np.flatnonzero for month indexing. |
| gravity_toolkit/gen_stokes.py | Refactor integration to complex/Euler form and np.einsum. |
| gravity_toolkit/gen_spherical_cap.py | Refactor summation to complex/Euler form and np.einsum; use np.radians. |
| gravity_toolkit/gen_point_load.py | Refactor to _complex_harmonics and np.einsum; use np.radians. |
| gravity_toolkit/gen_harmonics.py | Refactor integration/fourier paths to complex/Euler form and np.einsum. |
| gravity_toolkit/gen_disc_load.py | Refactor summation to complex/Euler form and np.einsum; use np.radians. |
| gravity_toolkit/clenshaw_summation.py | Refactor clenshaw summation to complex/Euler form and np.einsum; add _clenshaw. |
| geocenter/monte_carlo_degree_one.py | Large np.einsum refactor; add removal/redistribution options; add global attrs. |
| geocenter/model_degree_one.py | Refactor to np.einsum + Euler form; standardize radians conversions. |
| geocenter/kernel_degree_one.py | Refactor IMAT/G to np.einsum; standardize radians conversions. |
| geocenter/geocenter_spatial_maps.py | Use np.radians and np.where for phase wrapping. |
| geocenter/geocenter_monte_carlo.py | Use raw string formatting for LaTeX sigma label. |
| geocenter/delta_degree_one.py | Refactor to np.einsum + Euler form; standardize radians conversions. |
| geocenter/calc_degree_one.py | Refactor to np.einsum + Euler form; add remove/redistribute options. |
| doc/source/notebooks/GRACE-Spatial-Maps.ipynb | Update example code to np.radians. |
| doc/source/notebooks/GRACE-Spatial-Error.ipynb | Update example code to np.radians. |
| doc/source/notebooks/GRACE-Geostrophic-Maps.ipynb | Update example code to np.radians/np.flatnonzero; clear execution count. |
| doc/source/api_reference/sea_level_equation.rst | Add _clenshaw to API reference. |
| doc/source/api_reference/clenshaw_summation.rst | Add _clenshaw to API reference. |
| dealiasing/dealiasing_global_uplift.py | Use np.radians for colatitude conversion. |
| access/itsg_graz_grace_sync.py | Update ITSG server path from outgoing to pub. |
| .github/workflows/python-request.yml | Bump setup-pixi action version. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+95
to
+101
| # real (cosine) and imaginary (sine) components | ||
| Ylm = np.zeros((LMAX+1, MMAX+1), dtype=np.complex128) | ||
| # Truncating harmonics to degree and order LMAX | ||
| # removing coefficients below LMIN and above MMAX | ||
| mm = np.arange(0,MMAX+1) | ||
| clm = np.zeros((LMAX+1, MMAX+1)) | ||
| slm = np.zeros((LMAX+1, MMAX+1)) | ||
| clm[LMIN:LMAX+1,mm] = clm1[LMIN:LMAX+1,mm] | ||
| slm[LMIN:LMAX+1,mm] = slm1[LMIN:LMAX+1,mm] | ||
| # spherical harmonic degree and order | ||
| ll = np.arange(0,LMAX+1)[np.newaxis, :]# lmax+1 to include lmax | ||
| mm = np.arange(0,MMAX+1)[:, np.newaxis]# mmax+1 to include mmax | ||
| Ylm.imag[LMIN:LMAX+1,mm] = -clm1[LMIN:LMAX+1,mm].copy() | ||
| Ylm.real[LMIN:LMAX+1,mm] = -slm1[LMIN:LMAX+1,mm].copy() | ||
| dlm = np.einsum("l...,lm...->lm", np.sqrt((ll+1.0)*ll), Ylm) |
Comment on lines
+1730
to
+1732
| def __idiv__(self, other): | ||
| """In-place divide values from a ``spatial`` object""" | ||
| return self.scale(1.0 / other) |
Comment on lines
+1914
to
+1920
| def __idiv__(self, other): | ||
| """In-place divide values from a ``harmonics`` object""" | ||
| if isinstance(other, (int, float, np.ndarray)): | ||
| return self.scale(1.0 / other) | ||
| else: | ||
| return self.divide(other) | ||
|
|
Comment on lines
+176
to
+178
| # calculate summation for order LMAX | ||
| s_m = s_m_c[:,2*LMAX]*cos_m_phi[:,LMAX] + s_m_c[:,2*LMAX+1]*sin_m_phi[:,LMAX] | ||
| g = np.einsum("h...,p...->ph...", cs_m[:,LMAX], m_phi[:,LMAX]) | ||
| s_m = g.real |
Comment on lines
+182
to
+186
| m_phi[:,m] = cos_phi_2*m_phi[:,m+1] - m_phi[:,m+2] | ||
| a_m = np.sqrt((2.0*m + 3.0)/(2.0*m + 2.0)) | ||
| g = np.einsum("h...,p...->ph...", cs_m[:,m], m_phi[:,m]) | ||
| # update summation and discard imaginary component | ||
| s_m = a_m*u*s_m + g.real |
Comment on lines
+1721
to
+1724
| def __div__(self, other): | ||
| """Divide values from a ``spatial`` object""" | ||
| temp = self.copy() | ||
| return temp.scale(1.0 / other) |
Comment on lines
+1902
to
+1909
| def __div__(self, other): | ||
| """Divide values from a ``harmonics`` object""" | ||
| temp = self.copy() | ||
| if isinstance(other, (int, float, np.ndarray)): | ||
| return temp.scale(1.0 / other) | ||
| else: | ||
| return temp.divide(other) | ||
|
|
feat: add dunder methods to geocenter class
Comment on lines
+99
to
+101
| Ylm.imag[LMIN:LMAX+1,mm] = -clm1[LMIN:LMAX+1,mm].copy() | ||
| Ylm.real[LMIN:LMAX+1,mm] = -slm1[LMIN:LMAX+1,mm].copy() | ||
| dlm = np.einsum("l...,lm...->lm", np.sqrt((ll+1.0)*ll), Ylm) |
Comment on lines
+5
to
+9
| import pytest | ||
| import inspect | ||
| import pathlib | ||
| import numpy as np | ||
| import gravity_toolkit as gravtk |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
refactor: use
np.radiansto convert from degrees to radiansrefactor: use
np.hypotto calculate the sum of two squaresfeat: add dunder (magic) methods for mathematical operations