Add data types 'amplitude' and 'phase'

CHANGELOG.rst

@@ -6,6 +6,16 @@ Changelog
 """"""""""
 
 
+latest
+------
+
+- ``electrodes``:
+
+  - Receivers have a new parameter ``data_type``; by default this is
+    ``complex``, but in addition there are new the (experimental) options
+    ``amplitude`` and ``phase``.
+
+
 v1.3.2: Bugfix CLI-select
 -------------------------
 
@@ -60,7 +70,6 @@ v1.3.0: File-based computations
     ``interpolate_to_grid`` is the same as the current one.
 
 
-
 v1.2.1: Remove optimize & bug fix
 ---------------------------------
 

emg3d/electrodes.py

@@ -589,8 +589,17 @@ class Receiver(Wire):
         sources, such as is the case in a :class:`emg3d.surveys.Survey`.
 
     data_type : str
-        Data type of the measured responses. Currently implemented is only
-        ``'complex'``.
+        Data type of the measured responses. The data are always stored as
+        complex values, but the meaning of the real and imaginary part differs
+        depending on the data type. Currently implemented are:
+
+        - ``'complex'``: Complex values:        Real + j Imag
+        - ``'amplitude'``: Amplitude and phase:  Amp + j 0
+        - ``'phase'``: Amplitude and phase:      Pha + j 0
+
+        The preferred choice is complex, but the latter two are implemented for
+        the case where phase or amplitude information is not available or very
+        poor.
 
     """
 
@@ -601,7 +610,7 @@ def __init__(self, relative, data_type, **kwargs):
         """Initiate a receiver."""
 
         # Check data type is a known type.
-        if data_type.lower() != 'complex':
+        if data_type.lower() not in ['complex', 'amplitude', 'phase']:
             raise ValueError(f"Unknown data type '{data_type}'.")
 
         # Store relative, add a repr-addition.
@@ -623,6 +632,27 @@ def data_type(self):
         """Data type of the measured responses."""
         return self._data_type
 
+    def from_complex(self, complex_data):
+        """Return data in `data_type`-format from complex values."""
+
+        if self.data_type == 'amplitude':
+            return complex(abs(complex_data))
+
+        elif self.data_type == 'phase':
+            return complex(np.angle(complex_data))
+
+        else:
+            return complex_data
+
+    def derivative_chain(self, data, complex_data):
+        """Chain rule for data types other than complex."""
+
+        if self.data_type == 'amplitude':  # Amp + 0j
+            data *= np.real(complex_data.conj()/abs(complex_data))
+
+        elif self.data_type == 'phase':    # Pha + 0j
+            data *= np.real(-1j*complex_data.conj()/abs(complex_data)**2)
+
     def center_abs(self, source):
         """Returns points as absolute positions."""
         if self.relative:
@@ -656,8 +686,17 @@ class RxElectricPoint(Receiver, Point):
         sources, such as is the case in a :class:`emg3d.surveys.Survey`.
 
     data_type : str, default: 'complex'
-        Data type of the measured responses. Currently implemented is only the
-        default value.
+        Data type of the measured responses. The data are always stored as
+        complex values, but the meaning of the real and imaginary part differs
+        depending on the data type. Currently implemented are:
+
+        - ``'complex'``: Complex values:        Real + j Imag
+        - ``'amplitude'``: Amplitude and phase:  Amp + j 0
+        - ``'phase'``: Amplitude and phase:        0 + j Pha
+
+        The preferred choice is complex, but the latter two are implemented for
+        the case where phase or amplitude information is not available or very
+        poor.
 
     """
     _adjoint_source = TxElectricPoint
@@ -688,8 +727,17 @@ class RxMagneticPoint(Receiver, Point):
         sources, such as is the case in a :class:`emg3d.surveys.Survey`.
 
     data_type : str, default: 'complex'
-        Data type of the measured responses. Currently implemented is only the
-        default value.
+        Data type of the measured responses. The data are always stored as
+        complex values, but the meaning of the real and imaginary part differs
+        depending on the data type. Currently implemented are:
+
+        - ``'complex'``: Complex values:        Real + j Imag
+        - ``'amplitude'``: Amplitude and phase:  Amp + j 0
+        - ``'phase'``: Amplitude and phase:        0 + j Pha
+
+        The preferred choice is complex, but the latter two are implemented for
+        the case where phase or amplitude information is not available or very
+        poor.
 
     """
     _adjoint_source = TxMagneticPoint

emg3d/simulations.py

@@ -234,10 +234,15 @@ def __init__(self, survey, model, max_workers=4, gridding='single',
         # Sets self.model and self.gridding_opts.
         self._set_model(model, kwargs)
 
-        # Initiate synthetic data with NaN's if they don't exist.
+        # Initiate synthetic and complex 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(
                     data=np.full(self.survey.shape, np.nan+1j*np.nan))
+        if self.survey._anyrec_non_complex:
+            # Only required if there are non-complex data.
+            if 'complex' not in self.survey.data.keys():
+                self.survey._data['complex'] = self.data.observed.copy(
+                        data=np.full(self.survey.shape, np.nan+1j*np.nan))
 
         # `tqdm`-options; undocumented.
         # Can be used to, e.g., disable tqdm completely via:
@@ -343,6 +348,9 @@ def clean(self, what='computed'):
                     del self.data[key]
             self.data['synthetic'] = self.data.observed.copy(
                     data=np.full(self.survey.shape, np.nan+1j*np.nan))
+            if self.survey._anyrec_non_complex:
+                self.data['complex'] = self.data.observed.copy(
+                        data=np.full(self.survey.shape, np.nan+1j*np.nan))
             for name in ['_gradient', '_misfit']:
                 delattr(self, name)
                 setattr(self, name, None)
@@ -390,7 +398,7 @@ def to_dict(self, what='computed', copy=False):
 
         # Clean unwanted data if plain.
         if what == 'plain':
-            for key in ['synthetic', 'residual', 'weights']:
+            for key in ['synthetic', 'complex', 'residual', 'weights']:
                 if key in out['survey']['data'].keys():
                     del out['survey']['data'][key]
 
@@ -782,7 +790,14 @@ def collect_efield_inputs(inp):
 
             # Store responses at receiver locations.
             resp = self._get_responses(src, freq)
-            self.data['synthetic'].loc[src, :, freq] = resp
+            if self.survey._anyrec_non_complex:
+                syn = np.zeros_like(self.data.synthetic.loc[src, :, freq])
+                for i, rec in enumerate(self.survey.receivers.values()):
+                    syn[i] = rec.from_complex(resp[i])
+                self.data['synthetic'].loc[src, :, freq] = syn
+                self.data['complex'].loc[src, :, freq] = resp
+            else:
+                self.data['synthetic'].loc[src, :, freq] = resp
 
         # Print solver info.
         self.print_solver_info('efield', verb=self.verb)
@@ -794,7 +809,13 @@ def collect_efield_inputs(inp):
             self.data['observed'] = self.data['synthetic'].copy()
 
             # Add noise.
-            if kwargs.pop('add_noise', True):
+            if kwargs.get('add_noise', True):
+
+                # Ensure there are no non-complex data.
+                if self.survey._anyrec_non_complex:
+                    msg = "`add_noise` is only implemented for complex data."
+                    raise NotImplementedError(msg)
+
                 self.survey.add_noise(**kwargs)
 
     # OPTIMIZATION
@@ -1075,6 +1096,8 @@ def _get_rfield(self, source, frequency):
 
         # Get values for this source and frequency.
         grid = self.get_grid(source, frequency)
+        if self.survey._anyrec_non_complex:
+            complx = self.data.complex.loc[source, :, frequency].data
         residual = self.data.residual.loc[source, :, frequency].data
         weight = self.data.weights.loc[source, :, frequency].data
 
@@ -1091,6 +1114,11 @@ def _get_rfield(self, source, frequency):
             if np.isnan(residual[i]):
                 continue
 
+            # Apply chain rule to strength if data_type != complex.
+            if self.survey._anyrec_non_complex:
+                strength[i] = rec.from_complex(strength[i])
+                rec.derivative_chain(strength[i], complx[i])
+
             # Get absolute coordinates as fct of source.
             # (Only relevant in case of "relative" receivers.)
             coords = rec.coordinates_abs(self.survey.sources[source])
@@ -1178,6 +1206,14 @@ def collect_gfield_inputs(inp, vector=vector):
         for i, (src, freq) in enumerate(self._srcfreq):
             gfield = self._load(out[i][0], 'efield')
             resp = self._get_responses(src, freq, gfield)
+
+            # Apply chain rule to responses if data_type != complex.
+            if self.survey._anyrec_non_complex:
+                complx = self.data.complex.loc[src, :, freq].data
+                for ii, rec in enumerate(self.survey.receivers.values()):
+                    resp[ii] = rec.from_complex(resp[ii])
+                    rec.derivative_chain(resp[ii], complx[ii])
+
             self.data['jvec'].loc[src, :, freq] = resp
 
         return self.data['jvec'].data

emg3d/surveys.py

@@ -93,6 +93,10 @@ class Survey:
 
         If None, it will be initiated with NaN's.
 
+        Data can be provided in different formats, which are defined when
+        creating receiver instances. Consult their documentation for more info;
+        default is always ``'complex'``.
+
     noise_floor, relative_error : {float, ndarray}, default: None
         Noise floor and relative error of the data. They can be arrays of a
         shape which can be broadcasted to the data shape, e.g., (nsrc, 1, 1) or
@@ -703,6 +707,14 @@ def _rec_types_coord(self, source):
         indices = self._irec_types
         return [tuple(self._rec_coord[source][ind].T) for ind in indices]
 
+    @property
+    def _anyrec_non_complex(self):
+        """Boolean if any receiver has another data_type than complex."""
+        if getattr(self, '_anyrec_non_complex_bool', None) is None:
+            cmpl = [r.data_type != 'complex' for r in self.receivers.values()]
+            self._anyrec_non_complex_bool = any(cmpl)
+        return self._anyrec_non_complex_bool
+
 
 def random_noise(standard_deviation, mean_noise=0.0, ntype='white_noise'):
     r"""Return random noise for given inputs.