Foundation for inversion module (#339)

CHANGELOG.rst

@@ -9,6 +9,8 @@ Changelog
 latest
 ------
 
+- Created foundation for new module ``inversion``.
+
 Maintenance
 
   - Add notes for ``ipympl`` (interactive plots in modern Jupyter).

docs/api/index.rst

@@ -27,6 +27,7 @@ API reference
    surveys
    time
    utils
+   inversion/index
 
 
 .. grid:: 1
@@ -55,3 +56,7 @@ API reference
     .. grid-item-card::
 
         All sources and receivers are in :mod:`emg3d.electrodes`
+
+    .. grid-item-card::
+
+        Inversion: :mod:`emg3d.inversion`

docs/api/inversion/index.rst

@@ -0,0 +1,26 @@
+.. _inversionapi:
+
+Inversion
+#########
+
+.. automodapi:: emg3d.inversion
+   :no-inheritance-diagram:
+   :no-heading:
+
+.. toctree::
+   :maxdepth: 1
+   :hidden:
+
+   simpeg
+   pygimli
+
+.. grid:: 1
+    :gutter: 2
+
+    .. grid-item-card::
+
+        SimPEG(emg3d): :mod:`emg3d.inversion.simpeg`
+
+    .. grid-item-card::
+
+        pyGIMLi(emg3d): :mod:`emg3d.inversion.pygimli`

docs/api/inversion/pygimli.rst

@@ -0,0 +1,6 @@
+pyGIMLi(emg3d)
+==============
+
+.. automodapi:: emg3d.inversion.pygimli
+   :no-inheritance-diagram:
+   :no-heading:

docs/api/inversion/simpeg.rst

@@ -0,0 +1,6 @@
+SimPEG(emg3d)
+=============
+
+.. automodapi:: emg3d.inversion.simpeg
+   :no-inheritance-diagram:
+   :no-heading:

docs/conf.py

@@ -40,6 +40,8 @@
     "empymod": ("https://empymod.emsig.xyz/en/stable", None),
     "xarray": ("https://docs.xarray.dev/en/stable", None),
     "numba": ("https://numba.readthedocs.io/en/stable", None),
+    "pygimli": ("https://www.pygimli.org", None),
+    "simpeg": ("https://docs.simpeg.xyz/latest", None),
 }
 
 # ==== 2. General Settings ====

docs/manual/installation.rst

@@ -40,6 +40,9 @@ or via ``pip``:
 
    pip install emg3d[full]
 
+The wrappers provided in :mod:`emg3d.inversion` may require additional
+packages, please consult the :ref:`inversionapi`-API for more information.
+
 If you are new to Python we recommend using a Python distribution, which will
 ensure that all dependencies are met, specifically properly compiled versions
 of ``NumPy`` and ``SciPy``; we recommend using `Anaconda

emg3d/__init__.py

@@ -22,6 +22,7 @@
     RxElectricPoint, RxMagneticPoint,
 )
 from emg3d.fields import Field, get_source_field, get_magnetic_field
+from emg3d import inversion
 from emg3d.io import save, load, convert
 from emg3d.meshes import TensorMesh, construct_mesh
 from emg3d.models import Model

emg3d/_multiprocessing.py

@@ -43,6 +43,7 @@ def process_map(fn, *iterables, max_workers, **kwargs):
     execution.
 
     """
+    process_map.count += 1
 
     # Parallel
     if max_workers > 1 and tqdm is None:
@@ -64,6 +65,10 @@ def process_map(fn, *iterables, max_workers, **kwargs):
             iterable=map(fn, *iterables), total=len(iterables[0]), **kwargs))
 
 
+# Counter for processing map (used, e.g., for inversions).
+process_map.count = 0
+
+
 def solve(inp):
     """Thin wrapper of `solve` or `solve_source` for a `process_map`.
 

emg3d/inversion/__init__.py

@@ -0,0 +1,52 @@
+"""
+The inversion submodule of emg3d provides wrapper functionalities to use emg3d
+as a forward modelling kernel within third-party inversion frameworks. These
+third-party libraries are not included in emg3d, you have to install them
+separately.
+
+.. note::
+
+    This development is work in progress. Until an official *«inversion»*
+    release, everything can change or disappear without warning.
+
+Currently planned wrappers and their corresponding requirements:
+
+- SimPEG(emg3d): Requires `SimPEG <https://simpeg.xyz>`_ (*Simulation and
+  Parameter Estimation in Geophysics*).
+- pyGIMLi(emg3d): Requires `pyGIMLi <https://pygimli.org>`_ (*Geophysical
+  Inversion & Modelling Library*).
+
+"""
+# Copyright 2024 The emsig community.
+#
+# This file is part of emg3d.
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License.  You may obtain a copy
+# of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+# License for the specific language governing permissions and limitations under
+# the License.
+import importlib as _importlib
+
+
+submodules = [
+    'pygimli',
+    'simpeg',
+]
+
+__all__ = submodules
+
+
+def __dir__():
+    return __all__
+
+
+def __getattr__(name):
+    if name in submodules:
+        return _importlib.import_module(f"emg3d.inversion.{name}")

emg3d/inversion/pygimli.py

@@ -0,0 +1,43 @@
+"""
+Thin wrappers to use emg3d as a forward modelling kernel within the
+*Geophysical Inversion & Modelling Library* `pyGIMLi <https://pygimli.org>`_.
+
+It deals mainly with converting the data and model from the emg3d format to the
+pyGIMLi format and back, and creating the correct classes and functions as
+expected by a pyGIMLi inversion.
+"""
+# Copyright 2024 The emsig community.
+#
+# This file is part of emg3d.
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License.  You may obtain a copy
+# of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+# License for the specific language governing permissions and limitations under
+# the License.
+# import numpy as np
+
+from emg3d import utils  # io, _multiprocessing
+
+try:
+    import pygimli
+    # Add pygimli and pgcore to the emg3d.Report().
+    utils.OPTIONAL.extend(['pygimli', 'pgcore'])
+except ImportError:
+    pygimli = None
+
+__all__ = []
+
+
+def __dir__():
+    return __all__
+
+
+if pygimli is not None:
+    print("NOTE: pyGIMLi(emg3d) is in development.")

emg3d/inversion/simpeg.py

@@ -0,0 +1,47 @@
+"""
+Thin wrappers to use emg3d as a forward modelling kernel within the package
+*Simulation and Parameter Estimation in Geophysics* `SimPEG
+<https://simpeg.xyz>`_.
+
+It deals mainly with converting the data and model from the emg3d format to the
+SimPEG format and back, and creating the correct classes and functions as
+expected by a SimPEG inversion.
+"""
+# Copyright 2024 The emsig community.
+#
+# This file is part of emg3d.
+#
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License.  You may obtain a copy
+# of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+# License for the specific language governing permissions and limitations under
+# the License.
+# import numpy as np
+
+from emg3d import utils  # electrodes, meshes, models, simulations, surveys
+
+try:
+    import simpeg
+    # import discretize
+    # from simpeg.electromagnetics import frequency_domain as simpeg_fd
+    # Add simpeg to the emg3d.Report().
+    utils.OPTIONAL.extend(['simpeg',])
+except ImportError:
+    simpeg = None
+
+
+__all__ = []
+
+
+def __dir__():
+    return __all__
+
+
+if simpeg is not None:
+    print("NOTE: SimPEG(emg3d) is in development.")

emg3d/simulations.py

@@ -301,7 +301,7 @@ def __init__(self, survey, model, max_workers=4, gridding='single',
 
         # Initiate synthetic data with NaN's if they don't exist.
         if 'synthetic' not in self.survey.data.keys():
-            self.survey._data['synthetic'] = self.data.observed.copy(
+            self.survey.data['synthetic'] = self.data.observed.copy(
                     data=np.full(self.survey.shape, np.nan+1j*np.nan))
 
         # `tqdm`-options.

emg3d/utils.py

@@ -42,6 +42,8 @@
 
 __all__ = ['Report', 'EMArray', 'Timer']
 
+OPTIONAL = ['xarray', 'discretize', 'h5py', 'matplotlib', 'tqdm', 'IPython']
+
 
 def __dir__():
     return __all__
@@ -158,8 +160,7 @@ def __init__(self, add_pckg=None, ncol=3, text_width=80, sort=False):
         core = ['numpy', 'scipy', 'numba', 'emg3d', 'empymod']
 
         # Optional packages.
-        optional = ['xarray', 'discretize', 'h5py', 'matplotlib',
-                    'tqdm', 'IPython']
+        optional = OPTIONAL
 
         super().__init__(additional=add_pckg, core=core, optional=optional,
                          ncol=ncol, text_width=text_width, sort=sort)

requirements-dev.txt

@@ -8,6 +8,8 @@ xarray
 discretize
 matplotlib
 ipympl
+pygimli>=1.5.2
+simpeg>=0.22.1
 
 # SETUP RELATED
 setuptools_scm

setup.py

@@ -15,7 +15,7 @@
     author_email="[email protected]",
     url="https://emsig.xyz",
     license="Apache-2.0",
-    packages=["emg3d", "emg3d.cli"],
+    packages=["emg3d", "emg3d.inversion", "emg3d.cli"],
     classifiers=[
         "Development Status :: 5 - Production/Stable",
         "License :: OSI Approved :: Apache Software License",

tests/test_pygimli.py

@@ -0,0 +1,7 @@
+from emg3d import inversion
+from emg3d.inversion import pygimli as ipygimli
+
+
+def test_all_dir():
+    assert set(inversion.__all__) == set(dir(inversion))
+    assert set(ipygimli.__all__) == set(dir(ipygimli))

tests/test_simpeg.py

@@ -0,0 +1,7 @@
+from emg3d import inversion
+from emg3d.inversion import simpeg as isimpeg
+
+
+def test_all_dir():
+    assert set(inversion.__all__) == set(dir(inversion))
+    assert set(isimpeg.__all__) == set(dir(isimpeg))

tests/test_utils.py

@@ -5,6 +5,8 @@
 
 import scooby
 from emg3d import utils
+from emg3d.inversion import pygimli as ipygimli
+from emg3d.inversion import simpeg as isimpeg
 
 
 def test_known_class():
@@ -30,13 +32,20 @@ def dummy():
 def test_Report(capsys):
     out, _ = capsys.readouterr()  # Empty capsys
 
+    add = []
+
+    if ipygimli.pygimli:
+        add.extend(['pygimli', 'pgcore'])
+    if isimpeg.simpeg:
+        add.append('simpeg')
+
     # Reporting is now done by the external package scooby.
     # We just ensure the shown packages do not change (core and optional).
     out1 = utils.Report()
     out2 = scooby.Report(
-            core=['numpy', 'scipy', 'numba', 'emg3d'],
+            core=['numpy', 'scipy', 'numba', 'emg3d', 'empymod'],
             optional=['empymod', 'xarray', 'discretize', 'h5py',
-                      'matplotlib', 'tqdm', 'IPython'],
+                      'matplotlib', 'tqdm', 'IPython'] + add,
             ncol=4)
 
     # Ensure they're the same; exclude time to avoid errors.