diff --git a/docs/api/index.rst b/docs/api/index.rst index 4c03f995..5fbd6284 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -27,6 +27,7 @@ API reference surveys time utils + inversion .. grid:: 1 diff --git a/docs/api/inversion.rst b/docs/api/inversion.rst new file mode 100644 index 00000000..7b6271e5 --- /dev/null +++ b/docs/api/inversion.rst @@ -0,0 +1,6 @@ +Inversion +========= + +.. automodapi:: emg3d.inversion.pygimli + :no-inheritance-diagram: + :no-heading: diff --git a/docs/conf.py b/docs/conf.py index caeb9329..ec784ea9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -37,8 +37,10 @@ "scipy": ("https://docs.scipy.org/doc/scipy", None), "discretize": ("https://discretize.simpeg.xyz/en/main", None), "empymod": ("https://empymod.emsig.xyz/en/stable", None), - "xarray": ("https://xarray.pydata.org/en/stable", None), + "xarray": ("https://docs.xarray.dev/en/stable", None), "numba": ("https://numba.readthedocs.io/en/stable", None), + "numba": ("https://numba.readthedocs.io/en/stable", None), + "pygimli": ("https://www.pygimli.org", None), } # ==== 2. General Settings ==== diff --git a/emg3d/inversion/__init__.py b/emg3d/inversion/__init__.py new file mode 100644 index 00000000..309ab528 --- /dev/null +++ b/emg3d/inversion/__init__.py @@ -0,0 +1,17 @@ +# 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. + +from emg3d.inversion import pygimli diff --git a/emg3d/inversion/pygimli.py b/emg3d/inversion/pygimli.py new file mode 100644 index 00000000..c77b2479 --- /dev/null +++ b/emg3d/inversion/pygimli.py @@ -0,0 +1,286 @@ +# 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 + +try: + import pygimli +except ImportError: + pygimli = None + +from emg3d import models, utils + +__all__ = ['Kernel', 'Inversion', 'Jacobian'] + +# Add pygimli and pgcore to the emg3d.Report(). +utils.OPTIONAL.extend(['pygimli', 'pgcore']) + + +def __dir__(): + return __all__ + +# TODO: Create functions +# - convert_model(from/to) +# - convert_data(from/to) + + +@utils._requires('pygimli') +class Jacobian(pygimli.Matrix): + """Create a Jacobian operator for emg3d which is understood by pyGIMLi. + + This never builds the actual Jacobian, but provides functions to compute + the + + - Jacobian times a model vector ``Jm`` + (``jvec`` in emg3d, ``mult`` in pyGIMLi), and the + - Jacobian transposed times a data vector ``Jᵀd`` + (``jtvec`` in emg3d, ``transMult`` in pyGIMLi). + + + Parameters + ---------- + simulation : Simulation + The simulation; a :class:`emg3d.simulations.Simulation` instance. + + mesh : Mesh + The mesh; :func:`pygimli.meshtools.grid.createGrid` instance. + + """ + + def __init__(self, simulation, mesh): + """Initiate a new Jacobian instance.""" + super().__init__() + + # Store pointers to the emg3d-simulation and the pyGIMLi-mesh. + self.simulation = simulation + self.mesh = mesh + + # Store n-cols and n-rows. + self._cols = simulation.model.size + self._rows = simulation.survey.count * 2 + + def cols(self): + """The number of columns corresponds to the model size.""" + return self._cols + + def rows(self): + """The number of rows corresponds to 2x(data size), for [Re; Im].""" + return self._rows + + def mult(self, x): + """Multiply the Jacobian with a vector, Jm.""" + + self.simulation._count_jvec += 1 + + # Resort x to represent the model. + model = x[self.mesh().cellMarkers()] + + # Compute jvec. + jvec = self.simulation.jvec( + vector=np.reshape(model, self.simulation.model.shape, order='F') + ) + + # Get non-NaN data. + data = jvec[self.simulation.survey.isfinite] + + # Return real and imaginary parts stacked. + return np.hstack((data.real, data.imag)) + + def transMult(self, x): + """Multiply Jacobian transposed with a vector, Jᵀd = (dJᵀ)ᵀ.""" + self.simulation._count_jtvec += 1 + + # Cast finite [Re, Im] data from pyGIMLi into the emg3d format. + data = np.ones( + self.simulation.survey.shape, + dtype=self.simulation.data.observed.dtype + )*np.nan + x = np.asarray(x) + xl = x.size//2 + data[self.simulation.survey.isfinite] = x[:xl] + 1j*x[xl:] + + # Compute jtvec. + jtvec = self.simulation.jtvec(data).ravel('F') + + # Resort jtvec according to regions. + out = np.empty(jtvec.size) + out[self.mesh().cellMarkers()] = jtvec + return out + + def save(self, *args): + """There is no save for this pseudo-Jacobian.""" + pass + + +@utils._requires('pygimli') +class Kernel(pygimli.Modelling): + """Create a forward operator to use emg3d within a pyGIMLi inversion. + + + Parameters + ---------- + simulation : Simulation + The emg3d simulation a :class:`emg3d.simulations.Simulation` instance. + + """ + + def __init__(self, simulation, markers=None, pgthreads=2): + """Initialize the pyGIMLi(emg3d)-wrapper.""" + + pygimli.setThreadCount(pgthreads) + + # TODO move this to the Simulation class! + simulation._count_forward = 0 + simulation._count_jvec = 0 + simulation._count_jtvec = 0 + + # Check current limitation 1: isotropic. + mcase = simulation.model.case + if mcase != 'isotropic': + raise NotImplementedError( + f"pyGIMLi(emg3d) not implemented for {mcase} case." + ) + + # Check current limitation 2: conductivity. + mname = simulation.model.map.name + if mname != 'Conductivity': + raise NotImplementedError( + f"pyGIMLi(emg3d) not implemented for {mname} mapping." + ) + + # Initiate first pygimli.Modelling, which will do its magic. + super().__init__() + + # Store the simulation. + self.simulation = simulation + + # Translate discretize TensorMesh to pygimli-Grid. + mesh = pygimli.createGrid( + x=simulation.model.grid.nodes_x, + y=simulation.model.grid.nodes_y, + z=simulation.model.grid.nodes_z, + ) + + if markers is not None: + mesh.setCellMarkers(markers) + + self.setMesh(mesh) + + # Store J and set it + self.J = Jacobian(self.simulation, self.mesh) + self.setJacobian(self.J) + + def response(self, model): + """Create synthetic data for provided model.""" + + # Clean emg3d-simulation, so things are recomputed + self.simulation.clean('computed') + + # Replace model + self.simulation.model = models.Model( + grid=self.simulation.model.grid, + # Resort inversion model array to represent the actual model. + property_x=model[self.mesh().cellMarkers()], + mapping='Conductivity' + ) + + # Compute forward model and set initial residuals. + self.simulation._count_forward += 1 + _ = self.simulation.misfit + + # Return the responses + data = self.simulation.data.synthetic.data[ + self.simulation.survey.isfinite + ] + return np.hstack((data.real, data.imag)) + + def createStartModel(self, dataVals=None): + """Returns the model from the provided simulation.""" + return self.simulation.model.property_x.ravel('F') + + def createJacobian(self, model): + """Dummy to prevent pygimli.Modelling from doing it the hard way.""" + pass # do nothing + + +def _post_step(n, inv): + """TODO""" + + # TODO: save data, model, and everything to re-start inversion. + + # inv.chi2History + # inv.modelHistory + + sim = inv.fop.simulation + + kc = sim._count_forward + sim._count_jvec + sim._count_jtvec + sim.survey.data[f"it{n}"] = sim.survey.data.synthetic + cglsit = max(0, sim._count_jvec-1) + phi = inv.inv.getPhi() + if not hasattr(inv, 'lastphi'): + lastphi = "" + else: + lastphi = f"; Δϕ = {(1-phi/inv.lastphi)*100:.2f}%" + inv.lastphi = phi + pygimli.info( + f"{n}: " + f"χ² = {inv.inv.chi2():7.2f}; " + f"λ = {inv.inv.getLambda()}; " + f"#CGLS {cglsit:2d} ({kc:2d} solves); " + f"ϕ = {inv.inv.getPhiD():.2f} + {inv.inv.getPhiM():.2f}·λ = " + f"{phi:.2f}{lastphi}" + ) + + # Reset counters + sim._count_forward = 0 + sim._count_jvec = 0 + sim._count_jtvec = 0 + + +@utils._requires('pygimli') +class Inversion(pygimli.Inversion): + """TODO""" + def __init__(self, fop=None, inv=None, **kwargs): + super().__init__(fop=fop, inv=inv, **kwargs) + self._postStep = _post_step + + def run(self, dataVals=None, errorVals=None, **kwargs): + timer = utils.Timer() + pygimli.info(":: pyGIMLi(emg3d) START ::") + + # Take data from the survey if not provided. + if dataVals is None: + finite_data = self.fop.simulation.survey.finite_data() + dataVals = np.hstack([finite_data.real, finite_data.imag]) + + # Take the error from the survey if not provided. + if errorVals is None: + # TODO - IS THIS CORRECT? + std_dev_full = self.fop.simulation.survey.standard_deviation + std_dev = std_dev_full.data[self.fop.simulation.survey.isfinite] + errorVals = np.hstack([std_dev, std_dev]) / abs(dataVals) + + # TODO does it make any difference, is it needed? + # To completely ignore big errors + # => Test if it is actually necessary or not + errorVals[errorVals > 0.5] = 1e8 + + # Run the inversion + out = super().run(dataVals=dataVals, errorVals=errorVals, **kwargs) + + # Print passed time and exit + pygimli.info(f":: pyGIMLi(emg3d) END :: runtime = {timer.runtime}") + + return out diff --git a/emg3d/meshes.py b/emg3d/meshes.py index 66ed12e0..3ac7e0e6 100644 --- a/emg3d/meshes.py +++ b/emg3d/meshes.py @@ -180,6 +180,29 @@ def __eq__(self, mesh): return bool(equal) + def contains(self, mesh): + """Check if input mesh is a subset of this mesh. + + The provided ``mesh`` can be either an emg3d or a discretize + TensorMesh instance. + + """ + + # Check if mesh is of the same instance. + equal = mesh.__class__.__name__ == self.__class__.__name__ + + # Check dimensions. + if equal: + equal *= len(mesh.shape_cells) == len(self.shape_cells) + + # Check distances and origin. + if equal: + equal *= np.isin(mesh.nodes_x, self.nodes_x).all() + equal *= np.isin(mesh.nodes_y, self.nodes_y).all() + equal *= np.isin(mesh.nodes_z, self.nodes_z).all() + + return bool(equal) + def copy(self): """Return a copy of the TensorMesh.""" return self.from_dict(self.to_dict(True)) diff --git a/emg3d/surveys.py b/emg3d/surveys.py index 1ffcf437..d8a5b59d 100644 --- a/emg3d/surveys.py +++ b/emg3d/surveys.py @@ -713,6 +713,17 @@ def _rec_types_coord(self, source): indices = self._irec_types return [tuple(self._rec_coord[source][ind].T) for ind in indices] + @property + def isfinite(self): + """TODO: document!""" + if not hasattr(self, '_isfinite'): + self._isfinite = np.isfinite(self.data.observed.data) + return self._isfinite + + def finite_data(self, data='observed'): + """TODO: document!""" + return self.data[data].data[self.isfinite] + def random_noise(standard_deviation, mean_noise=0.0, ntype='white_noise'): r"""Return random noise for given inputs. diff --git a/emg3d/utils.py b/emg3d/utils.py index 356707dc..59b734eb 100644 --- a/emg3d/utils.py +++ b/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) @@ -171,6 +172,7 @@ class Timer: def __init__(self): """Initiate timer with a performance counter.""" self._t0 = perf_counter() + self._previous = 0 def __repr__(self): """Simple representation.""" @@ -191,6 +193,13 @@ def runtime(self): """Return elapsed time as hh:mm:ss string.""" return str(timedelta(seconds=np.round(self.elapsed))) + @property + def laptime(self): + """Return time of this lap as hh:mm:ss string.""" + previous = self._previous + self._previous = self.elapsed + return str(timedelta(seconds=np.round(self._previous-previous))) + @property def elapsed(self): """Return elapsed time in seconds.""" diff --git a/setup.py b/setup.py index d7a0d3ca..504fd6c0 100644 --- a/setup.py +++ b/setup.py @@ -15,7 +15,7 @@ author_email="info@emsig.xyz", 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", diff --git a/tests/test_meshes.py b/tests/test_meshes.py index bbaa2fe1..701d7c8e 100644 --- a/tests/test_meshes.py +++ b/tests/test_meshes.py @@ -101,6 +101,14 @@ def test_TensorMesh(self): assert mesh.cell_volumes.sum() > 69046392 + # Check contains. + tgrid = meshes.TensorMesh( + [cgrid.h[0][:3], cgrid.h[1][:5], cgrid.h[2]], origin=cgrid.origin) + assert emg3dgrid.contains(tgrid) + newgrid = meshes.TensorMesh( + [np.ones(3), np.ones(3), np.ones(3)], np.zeros(3)) + assert not emg3dgrid.contains(newgrid) + def test_TensorMesh_repr(self): # Create some dummy data grid = meshes.TensorMesh(