Update CLI for receiver_interpolation (#224)

CHANGELOG.rst

@@ -6,61 +6,63 @@ Changelog
 """"""""""
 
 
-latest
-------
+v1.1.0: Adjoint-fix for electric receivers
+------------------------------------------
 
-**Fields**
-
-- ``get_receiver`` has a new keyword ``method``, which can be ``'cubic'`` or
-  ``'linear'``; default is the former, which is the same behaviour us before.
-  However, if you want to compute the gradient, you should set it to
-  ``'linear'`` in your Simulation parameters. Otherwise the adjoint-state
-  gradient will not exactly be the adjoint state. Note: This will change once
-  the adjoint of a cubic interpolation is implemented as source function.
-
-- ``get_source_field``: If ``frequency=None``, it returns new the real-valued,
-  frequency-independent, source vector.
-
-
-**Electrodes**
+**2021-06-30**
 
-- Re-introduced the point source as ``TxElectricPoint``.
+This release contains, besides the usual small bugfixes, typos, and small
+improvements, an important fix for ``optimize.gradient``. Keep in mind that
+while the forward modelling is regarded as stable, the ``optimize`` module is
+still work in progress.
 
+The fixes with regard to ``optimize.gradient`` ensure that the gradient is
+indeed using the proper adjoint to back-propagate the field. This is currently
+*only* given for electric receivers, not yet for magnetic receivers. These
+improvement happened mainly thanks to the help of Seogi (@sgkang).
 
-**Simulations**
-
-- New keyword ``receiver_interpolation``, which corresponds to the ``method``
-  in ``get_receiver`` (see above). Cubic is more precise. However, if you are
-  interested in the gradient, you need to choose 'linear' at the moment, as
-  there are only linearly interpolated source functions. To be the proper
-  adjoint for the gradient the receiver has to be interpolated linearly too.
-
-- If ``gridding`` is ``'same'`` or ``'input'``, it new checks if the provided
-  grid is a sensible grid for emg3d; if not, it throws a warning.
-
-
-**Meshes**
+The changes in more detail:
 
-- New function ``check_grid`` to verify if a given grid is good for emg3d.
+- ``fields``:
 
+  - ``get_receiver`` has a new keyword ``method``, which can be ``'cubic'`` or
+    ``'linear'``; default is the former, which is the same behaviour as before.
+    However, if you want to compute the gradient, you should set it to
+    ``'linear'`` in your Simulation parameters. Otherwise the adjoint-state
+    gradient will not exactly be the adjoint state.
+  - ``get_source_field`` returns new the real-valued, frequency-independent
+    source vector if ``frequency=None``.
+  - ``get_source_field`` uses the adjoint of trilinear interpolation for point
+    sources (new). For dipoles and wires it the source is distributed onto the
+    cells as fraction of the source length (as before).
 
-**Optimize**
+- ``electrodes``: Re-introduced the point source as ``TxElectricPoint``.
 
-This release contains various improvements to the adjoint-state gradient.
-Electric receivers work fine, but there are still some remaining issues with
-respect to magnetic receivers.
+- ``simulations.Simulation``:
 
-- ``gradient``:
+  - New keyword ``receiver_interpolation``, which corresponds to the ``method``
+    in ``get_receiver`` (see above). Cubic is more precise. However, if you are
+    interested in the gradient, you need to choose linear interpolation at the
+    moment, as the point source is the adjoint of linear interpolation. To be
+    the proper adjoint for the gradient the receiver has to be interpolated
+    linearly too.
+  - If ``gridding`` is ``'same'`` or ``'input'``, it checks now if the provided
+    grid is a sensible grid for emg3d; if not, it throws a warning.
 
-  - Changed order when going from computational grid to inversion grid.
-    Changing the grids at the field stage (cubic interpolation) seems to be
-    better than changing at the cell-averaged stage:
+- ``meshes``: New function ``check_grid`` to verify if a given grid is good for
+  emg3d.
 
-    New: field_comp -> field_inv -> cells_inv
-    Old: field_comp -> cells_comp -> cells_inv
+- ``optimize.gradient``: Changed order when going from computational grid to
+  inversion grid. Changing the grids at the field stage (cubic interpolation)
+  seems to be better than changing at the cell-averaged stage::
 
+      New: field_comp -> field_inv -> cells_inv
+      Old: field_comp -> cells_comp -> cells_inv
 
-- Various small fixes to docs etc.
+- ``cli``: Uses now by default linear receiver interpolation if the
+  ``gradient`` is wanted (new), otherwise it uses cubic interpolation (as
+  before). The new keyword ``receiver_interpolation`` of the simulation can be
+  set in the parameter file, which overwrites the described default behaviour.
 
 
 v1.0.0: Stable API
@@ -1097,7 +1099,7 @@ were removed, however.
   - ``utils.get_cell_numbers`` to get good values of number of cells for given
     primes.
 
-- Speed-up ``njitted.volume_average`` significantly thanks to @jcapriot.
+- Speed-up ``njitted.volume_average`` significantly thanks to Joe (@jcapriot).
 - Bugfixes and other minor things:
 
   - Abort if l2-norm is NaN (only works for MG).

emg3d/cli/parser.py

@@ -183,6 +183,14 @@ def parse_config_file(args_dict):
         _ = all_sim.pop(key)
         simulation[key] = cfg.getfloat('simulation', key)
 
+    key = 'receiver_interpolation'
+    if cfg.has_option('simulation', key):
+        _ = all_sim.pop(key)
+        simulation[key] = cfg.get('simulation', key)
+    elif term['function'] == 'gradient':
+        # Default is 'cubic' - gradient needs 'linear'
+        simulation[key] = 'linear'
+
     # Ensure no keys are left.
     if all_sim:
         raise TypeError(

emg3d/fields.py

@@ -413,9 +413,9 @@ def get_source_field(grid, source, frequency, **kwargs):
           :math:`s = \mathrm{i}\omega = 2\mathrm{i}\pi f` (complex);
         - ``frequency < 0``: Laplace domain, hence
           :math:`s = f` (real).
-        - ``frequency == None``: Returns the real-valued,
-          frequency-independent, source vector. This basically excludes the
-          multiplication by :math:`-\mathrm{i}\omega\mu_0`.
+        - ``frequency == None``: Returns the real-valued, frequency-independent
+          source vector. This basically excludes the multiplication by
+          :math:`-\mathrm{i}\omega\mu_0`.
 
     strength : {float, complex}, default: 1.0
         Source strength (A), put through to

emg3d/simulations.py

@@ -162,7 +162,7 @@ class Simulation:
         interested in the gradient, you need to choose 'linear' at the moment,
         as there are only linearly interpolated source functions. To be the
         proper adjoint for the gradient the receiver has to be interpolated
-        linearly too. (This will change in the future.)
+        linearly too.
 
     """
 
@@ -188,8 +188,8 @@ def __init__(self, survey, model, max_workers=4, gridding='single',
         self.verb = kwargs.pop('verb', 0)
         self.name = kwargs.pop('name', None)
         self.info = kwargs.pop('info', None)
-        self.receiver_interpolation = kwargs.pop(   # Remove once we have
-                'receiver_interpolation', 'cubic')  # cubic source fct.
+        self.receiver_interpolation = kwargs.pop(
+                'receiver_interpolation', 'cubic')
 
         # Assemble solver_opts.
         self.solver_opts = {

tests/test_cli.py

@@ -204,10 +204,13 @@ def test_simulation(self, tmpdir):
         args_dict = self.args_dict.copy()
         args_dict['config'] = config
         cfg, term = cli.parser.parse_config_file(args_dict)
-        assert cfg['simulation_options']['max_workers'] == 5
-        assert cfg['simulation_options']['gridding'] == 'fancything'
-        assert cfg['simulation_options']['name'] == "PyTest simulation"
-        assert cfg['simulation_options']['min_offset'] == 1320.0
+        sim_opts = cfg['simulation_options']
+        assert sim_opts['max_workers'] == 5
+        assert sim_opts['gridding'] == 'fancything'
+        assert sim_opts['name'] == "PyTest simulation"
+        assert sim_opts['min_offset'] == 1320.0
+        with pytest.raises(KeyError, match="receiver_interpolation"):
+            assert sim_opts['receiver_interpolation'] == 'linear'
 
         with pytest.raises(TypeError, match="Unexpected parameter in"):
             with open(config, 'a') as f:
@@ -216,6 +219,26 @@ def test_simulation(self, tmpdir):
             args_dict['config'] = config
             _ = cli.parser.parse_config_file(args_dict)
 
+        # Ensure it sets interpolation to linear for gradient
+        args_dict = self.args_dict.copy()
+        args_dict['gradient'] = True
+        cfg, term = cli.parser.parse_config_file(args_dict)
+        sim_opts = cfg['simulation_options']
+        assert sim_opts['receiver_interpolation'] == 'linear'
+
+        # Ensure config file overrides that
+        config = os.path.join(tmpdir, 'emg3d.cfg')
+        with open(config, 'w') as f:
+            f.write("[simulation]\n")
+            f.write("receiver_interpolation=cubic")
+
+        args_dict = self.args_dict.copy()
+        args_dict['config'] = config
+        args_dict['gradient'] = True
+        cfg, term = cli.parser.parse_config_file(args_dict)
+        sim_opts = cfg['simulation_options']
+        assert sim_opts['receiver_interpolation'] == 'cubic'
+
     def test_solver(self, tmpdir):
 
         # Write a config file.