Optimize

CHANGELOG.rst

@@ -5,11 +5,26 @@ Changelog
 recent versions
 """""""""""""""
 
-*latest*
---------
+
+*v0.12.0* : Survey & Simulation
+-------------------------------
+
+**2020-07-25**
 
 This is a big release with many new features, and unfortunately not completely
-backwards compatible.
+backwards compatible. The main new features are the new **Survey** and
+**Simulation** classes, as well as some initial work for **optimization**
+(misfit, gradient). Also, a **Model** can now be a resistivity model, a
+conductivity model, or the logarithm (natural or base 10) therefore. Receivers
+can now be arbitrarily rotated, just as the sources. In addition to the
+existing **soft-dependencies** `empymod`, `discretize`, and `h5py` there are
+the new soft-dependencies `xarray` and `tqm`; `discretize` is now much tighter
+integrated. For the new survey and simulation classes `xarray` is a required
+dependency. However, the only hard dependency remain `scipy` and `numba`, if
+you use `emg3d` purely as a solver. Data reading and writing has new a
+JSON-backend, in addition to the existing HDF5 and NumPy-backends.
+
+In more detail:
 
 - Modules:
 
@@ -19,16 +34,17 @@ backwards compatible.
     - Class `surveys.Dipole`, which defines electric or magnetic point dipoles
       and finite length dipoles.
 
-  - `simulations` (**new**; soft-dependencies `discretize` & `tqdm`):
+  - `simulations` (**new**; requires `xarray`; soft-dependency `tqdm`):
 
-    - Class `simulations.Simulation`, which combines a survey with a model.
-      A simulation computes the e-field (and h-field) asynchronously using
-      `concurrent.futures`. To do so it creates the required meshes, source and
-      frequency-dependent, interpolates the model accordingly, and computes the
-      source-fields. If `tqdm` is installed it displays a progress bar for the
-      asynchronous computation.
+    - Class `simulations.Simulation`, which combines a survey with a model. A
+      simulation computes the e-field (and h-field) asynchronously using
+      `concurrent.futures`. This class will include automatic, source- and
+      frequency-dependent gridding in the future. If `tqdm` is installed it
+      displays a progress bar for the asynchronous computation. Note that the
+      simulation class has still some limitations, consult the class
+      documentation.
 
-  - `models` & `maps`:
+  - `models`:
 
     - Model instances take new the parameters `property_{x;y;z}` instead of
       `res_{x;y;z}`. The properties can be either resistivity, conductivity, or
@@ -38,11 +54,17 @@ backwards compatible.
       the moment. The attributes `model.res_{x;y;z}` are still available too,
       but equally **deprecated**. However, it is **no longer possible to
       assign values to these attributes**, which is a **backwards
-      incompatible** change. (The corresponding mappings that make this
-      possible live in `maps`).
+      incompatible** change.
     - A model knows now how to interpolate itself from its grid to another grid
       (`interpolate2grid`).
 
+  - `maps`:
+
+    - **New** mappings for `models.Model` instances: The mappings take care of
+      how to transform the investigation variable to conductivity and back, and
+      how it affects its derivative.
+    - **New** interpolation routine `edges2cellaverages`.
+
   - `fields`:
 
     - Function `get_receiver_response` (**new**), which returns the response
@@ -69,15 +91,25 @@ backwards compatible.
     - `meshes.TensorMesh` **new** inherits from `discretize` if installed.
     - Added `__eq__` to `models.TensorMesh` to compare meshes.
 
+  - `optimize` (**new**)
+
+    - Functionalities related to inversion (data misfit, gradient, data
+      weighting, and depth weighting). This module is in an early stage, and
+      the API will likely change in the future. Current functions are `misfit`,
+      `gradient` (using the adjoint-state method), and `data_weighting`. These
+      functionalities are best accessed through the `Simulation` class.
+
 - Dependencies:
 
-  - `empymod` is new a soft dependency, only required for `utils.Fourier`
-    (time-domain modelling).
-  - New soft dependency `xarray` for the `Survey` class.
+  - `empymod` is now a soft dependency (no longer a hard dependency), only
+    required for `utils.Fourier` (time-domain modelling).
+  - Existing soft dependency `discretize` is now baked straight into `meshes`.
+  - New soft dependency `xarray` for the `Survey` class (and therefore also for
+    the `Simulation` class and the `optimize` module).
   - New soft dependency `tqdm` for nice progress bars in asynchronous
     computation.
 
-- Deprecations and removals:
+- **Deprecations** and removals:
 
   - Removed deprecated functions `data_write` and `data_read`.
   - Removed all deprecated functions from `utils`.

docs/code.rst

@@ -31,6 +31,10 @@ Code
    :no-inheritance-diagram:
    :headings:=^
 
+.. automodapi:: emg3d.optimize
+   :no-inheritance-diagram:
+   :headings:=^
+
 .. automodapi:: emg3d.utils
    :no-inheritance-diagram:
    :headings:=^

docs/references.rst

@@ -50,6 +50,10 @@ References
 .. [PlDM07] Plessix, R.-É., M. Darnet, and W. A. Mulder, 2007, An approach for
    3D multisource, multifrequency CSEM modeling: Geophysics, 72, SM177--SM184;
    DOI: `10.1190/1.2744234 <https://doi.org/10.1190/1.2744234>`_.
+.. [PlMu08] Plessix, R.-É., and W. A. Mulder, 2008, Resistivity imaging with
+   controlled-source electromagnetic data: depth and data weighting: Inverse
+   Problems, 24, no. 3, 034012; DOI: `10.1088/0266-5611/24/3/034012
+   <https://doi.org/10.1088/0266-5611/24/3/034012>`_.
 .. [SlHM10] Slob, E., J. Hunziker, and W. A. Mulder, 2010, Green's tensors for
    the diffusive electric field in a VTI half-space: PIER, 107, 1--20: DOI:
    `10.2528/PIER10052807 <http://doi.org/10.2528/PIER10052807>`_.

emg3d/__init__.py

@@ -23,6 +23,8 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
+import warnings
+
 # Import modules
 from emg3d import io
 from emg3d import maps
@@ -32,6 +34,7 @@
 from emg3d import meshes
 from emg3d import models
 from emg3d import surveys
+from emg3d import optimize
 from emg3d import simulations
 
 # Import most important functions and classes
@@ -44,11 +47,17 @@
 from emg3d.fields import Field  # noqa
 from emg3d.models import Model  # noqa
 from emg3d.utils import Fourier  # noqa
+from emg3d.surveys import Survey  # noqa
 from emg3d.meshes import TensorMesh  # noqa
+from emg3d.simulations import Simulation  # noqa
 from emg3d.fields import get_source_field, get_receiver, get_h_field  # noqa
 
 __all__ = ['solve', 'solver', 'utils', 'io', 'fields', 'maps', 'meshes',
-           'models', 'Report', 'save', 'load', 'surveys', 'simulations']
+           'models', 'Report', 'save', 'load', 'surveys', 'simulations',
+           'optimize']
+
+# Ensure users see Deprecation Warnings, but just once.
+warnings.filterwarnings("once", module='emg3d', category=DeprecationWarning)
 
 # Version defined in utils, so we can easier use it within the package itself.
 __version__ = utils.__version__

emg3d/core.py

@@ -23,7 +23,6 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
-
 import numba as nb
 import numpy as np
 

emg3d/fields.py

@@ -23,9 +23,9 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
+from copy import deepcopy
 
 import numpy as np
-from copy import deepcopy
 from scipy.constants import mu_0
 
 from emg3d import maps, models, utils

emg3d/io.py

@@ -22,21 +22,20 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
-
 import os
 import json
 import warnings
-import numpy as np
 from datetime import datetime
 
-from emg3d import fields, models, utils, meshes, surveys, simulations
+import numpy as np
 
 try:
     import h5py
 except ImportError:
     h5py = ("'.h5'-files require `h5py`. Install it via\n"
             "`pip install h5py` or `conda install -c conda-forge h5py`.")
 
+from emg3d import fields, models, utils, meshes, surveys, simulations
 
 __all__ = ['save', 'load']
 
@@ -358,8 +357,9 @@ def _dict_serialize(inp, out=None, collect_classes=False):
                 try:
                     to_dict = {'hx': value.hx, 'hy': value.hy, 'hz': value.hz,
                                'x0': value.x0, '__class__': name}
-                except AttributeError:  # Gracefully fail.
-                    print(f"* WARNING :: Could not serialize <{key}>")
+                except AttributeError as e:  # Gracefully fail.
+                    print(f"* WARNING :: Could not serialize <{key}>.\n"
+                          f"             {e}")
                     continue
 
             # If we are in the root-directory put them in their own category.
@@ -434,8 +434,10 @@ def _dict_deserialize(inp, first_call=True):
                     inp[key] = inst.from_dict(value)
                     continue
 
-                except (AttributeError, KeyError):  # Gracefully fail.
-                    print(f"* WARNING :: Could not de-serialize <{key}>")
+                except (NotImplementedError, AttributeError, KeyError) as e:
+                    # Gracefully fail.
+                    print(f"* WARNING :: Could not de-serialize <{key}>.\n"
+                          f"             {e}")
 
             # In no __class__-key or de-serialization fails, use recursion.
             _dict_deserialize(value, False)

emg3d/maps.py

@@ -23,21 +23,19 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
-
 import numba as nb
 import numpy as np
 from scipy import interpolate, ndimage
 
-from emg3d import fields
+__all__ = ['grid2grid', 'interp3d', 'MapConductivity', 'MapLgConductivity',
+           'MapLnConductivity', 'MapResistivity', 'MapLgResistivity',
+           'MapLnResistivity', 'edges2cellaverages']
 
 # Numba-settings
 _numba_setting = {'nogil': True, 'fastmath': True, 'cache': True}
 
-__all__ = ['grid2grid', 'interp3d', 'MapConductivity', 'MapLgConductivity',
-           'MapLnConductivity', 'MapResistivity', 'MapLgResistivity',
-           'MapLnResistivity']
-
 
+# INTERPOLATIONS
 def grid2grid(grid, values, new_grid, method='linear', extrapolate=True,
               log=False):
     """Interpolate `values` located on `grid` to `new_grid`.
@@ -115,7 +113,7 @@ def grid2grid(grid, values, new_grid, method='linear', extrapolate=True,
                        extrapolate, log)
 
         # Return a field instance.
-        return fields.Field(fx, fy, fz)
+        return values.__class__(fx, fy, fz)
 
     # If values is a particular field, ensure method is not 'volume'.
     if not np.all(grid.vnC == values.shape) and method == 'volume':
@@ -275,7 +273,7 @@ def interp3d(points, values, new_points, method, fill_value, mode):
     return new_values
 
 
-# Maps
+# MAPS
 class _Map:
     """Maps variable `x` to computational variable `σ` (conductivity)."""
 
@@ -428,7 +426,7 @@ def derivative(self, gradient, conductivity):
         gradient /= -conductivity
 
 
-# Volume averaging
+# VOLUME AVERAGING
 @nb.njit(**_numba_setting)
 def volume_average(edges_x, edges_y, edges_z, values,
                    new_edges_x, new_edges_y, new_edges_z, new_values, new_vol):
@@ -463,9 +461,9 @@ def volume_average(edges_x, edges_y, edges_z, values,
     """
 
     # Get the weights and indices for each direction.
-    wx, ix_in, ix_out = _volume_avg_weights(edges_x, new_edges_x)
-    wy, iy_in, iy_out = _volume_avg_weights(edges_y, new_edges_y)
-    wz, iz_in, iz_out = _volume_avg_weights(edges_z, new_edges_z)
+    wx, ix_in, ix_out = _volume_average_weights(edges_x, new_edges_x)
+    wy, iy_in, iy_out = _volume_average_weights(edges_y, new_edges_y)
+    wz, iz_in, iz_out = _volume_average_weights(edges_z, new_edges_z)
 
     # Loop over the elements and sum up the contributions.
     for iz, w_z in enumerate(wz):
@@ -485,8 +483,8 @@ def volume_average(edges_x, edges_y, edges_z, values,
 
 
 @nb.njit(**_numba_setting)
-def _volume_avg_weights(x1, x2):
-    """Returns the weights for the volume averaging technique.
+def _volume_average_weights(x1, x2):
+    """Returaveragethe weights for the volume averaging technique.
 
 
     Parameters
@@ -552,3 +550,59 @@ def _volume_avg_weights(x1, x2):
             ii += 1
 
     return hs[:ii], ix1[:ii], ix2[:ii]
+
+
+# EDGES => CENTERS
+@nb.njit(**_numba_setting)
+def edges2cellaverages(ex, ey, ez, vol, out_x, out_y, out_z):
+    r"""Interpolate fields defined on edges to volume-averaged cell values.
+
+
+    Parameters
+    ----------
+    ex, ey, ez : ndarray
+        Electric fields in x-, y-, and z-directions, as obtained from
+        :class:`emg3d.fields.Field`.
+
+    vol : ndarray
+        Volumes of the grid, as obtained from :class:`emg3d.meshes.TensorMesh`.
+
+    out_x, out_y, out_z : ndarray
+        Arrays where the results are placed (per direction).
+
+    """
+
+    # Get dimensions
+    nx, ny, nz = vol.shape
+
+    # Loop over dimensions.
+    for iz in range(nz+1):
+        izm = max(0, iz-1)
+        izp = min(nz-1, iz)
+
+        for iy in range(ny+1):
+            iym = max(0, iy-1)
+            iyp = min(ny-1, iy)
+
+            for ix in range(nx+1):
+                ixm = max(0, ix-1)
+                ixp = min(nx-1, ix)
+
+                # Multiply field by volume/4.
+                if ix < nx:
+                    out_x[ix, iym, izm] += vol[ix, iym, izm]*ex[ix, iy, iz]/4
+                    out_x[ix, iyp, izm] += vol[ix, iyp, izm]*ex[ix, iy, iz]/4
+                    out_x[ix, iym, izp] += vol[ix, iym, izp]*ex[ix, iy, iz]/4
+                    out_x[ix, iyp, izp] += vol[ix, iyp, izp]*ex[ix, iy, iz]/4
+
+                if iy < ny:
+                    out_y[ixm, iy, izm] += vol[ixm, iy, izm]*ey[ix, iy, iz]/4
+                    out_y[ixp, iy, izm] += vol[ixp, iy, izm]*ey[ix, iy, iz]/4
+                    out_y[ixm, iy, izp] += vol[ixm, iy, izp]*ey[ix, iy, iz]/4
+                    out_y[ixp, iy, izp] += vol[ixp, iy, izp]*ey[ix, iy, iz]/4
+
+                if iz < nz:
+                    out_z[ixm, iym, iz] += vol[ixm, iym, iz]*ez[ix, iy, iz]/4
+                    out_z[ixp, iym, iz] += vol[ixp, iym, iz]*ez[ix, iy, iz]/4
+                    out_z[ixm, iyp, iz] += vol[ixm, iyp, iz]*ez[ix, iy, iz]/4
+                    out_z[ixp, iyp, iz] += vol[ixp, iyp, iz]*ez[ix, iy, iz]/4

emg3d/meshes.py

@@ -22,12 +22,11 @@
 # License for the specific language governing permissions and limitations under
 # the License.
 
+from copy import deepcopy
 
 import numpy as np
-from copy import deepcopy
 from scipy import optimize
 
-# Import soft dependencies.
 try:
     import discretize.TensorMesh as dTensorMesh
 except ImportError:
@@ -57,9 +56,9 @@ def __init__(self, h, x0):
         self.x0 = x0
 
         # Width of cells.
-        self.hx = h[0]
-        self.hy = h[1]
-        self.hz = h[2]
+        self.hx = np.array(h[0])
+        self.hy = np.array(h[1])
+        self.hz = np.array(h[2])
 
         # Cell related properties.
         self.nCx = int(self.hx.size)