Merge branch 'dev' of https://github.com/gimli-org/gimli into dev

Author: carsten-forty2

.github/workflows/main.yml

@@ -39,7 +39,7 @@ jobs:
       OPENBLAS_CORETYPE: "ARMV8" # current bug in OpenBlas core detection
     steps:
       - name: Run pg.test()
-        run: python3 -c "import pygimli; pygimli.test(show=False)"
+        run: python3 -c "import pygimli; pygimli.test(show=False, abort=True)"
         working-directory: source
       - name: Install as development package
         working-directory: source

INSTALLATION.rst

@@ -21,14 +21,13 @@ Installation
    </p><br>
 
 
-On all platforms, we recommend to install pyGIMLi via the conda package manager
-contained in the Anaconda Python distribution. For details on how to install
-Anaconda, see `this page <https://docs.anaconda.com/anaconda/install/>`_.
-
-Note that Anaconda comes with many (great) packages, many of which you likely
-will not use. If you want to save space, you can install the `light-weight
-version Miniconda
-<https://docs.anaconda.com/free/miniconda/miniconda-install/>`_.
+On all platforms, we recommend install pyGIMLi via the conda package manager
+and the cleanest way is to create a new environment.
+Besides large installers like `Anaconda <https://anaconda.org/>`_
+that come with many (great) packages, you can save space and install a 
+light-weight installer like
+`Miniconda <https://docs.anaconda.com/free/miniconda/miniconda-install/>`_
+or `Miniforge <https://conda-forge.org/download/>`_.
 
 .. note::
 
@@ -91,12 +90,19 @@ follow these steps:
 
 pyGIMLi on Google Colab
 -----------------------
-Even though still experimental, pyGIMLi can be run on Google Colab without any
-installation on your own computer. Just create a new Notebook and install the
-pyGIMLi package via pip:
+pyGIMLi can be run on Google Colab without any installation on your own
+computer. Just create a new Notebook and install pyGIMLi via pip:
+
+.. code:: python
+
+    !pip install pygimli tetgen
+
+It turns out that there are some packages preinstalled that lead to some
+incompatibl numpy version, so you might have to uninstall them first.
 
 .. code:: python
 
+    !pip uninstall -y numba tensorflow pytensor thinc
     !pip install pygimli tetgen
 
 Staying up-to-date

core/src/inversion.cpp

@@ -248,7 +248,6 @@ const RVector & RInversion::run(){ ALLOW_PYTHON_THREADS
 
 bool RInversion::oneStep() {
     iter_++;
-
     deltaModelIter_.resize(model_.size());
     deltaModelIter_ *= 0.0;
     deltaDataIter_ = (tD_->trans(data_) - tD_->trans(response_));
@@ -257,19 +256,19 @@ bool RInversion::oneStep() {
         if (verbose_) std::cout << "sum(abs(deltaDataIter_)) == Zero" << std::endl;
         return false;
     }
-//    Vec deltaModel0(model_.size());
+    //    Vec deltaModel0(model_.size());
     Vec modelNew(   model_.size());
     Vec responseNew( data_.size());
     Vec roughness(this->constraintsCount(), 0.0);
 
     this->checkJacobian((recalcJacobian_ && iter_ > 1) || jacobiNeedRecalc_);
 
     // if ((recalcJacobian_ && iter_ > 1) || jacobiNeedRecalc_ ) {
-    //     Stopwatch swatch(true);
-    //     if (verbose_) std::cout << "Recalculating Jacobian matrix ...";
-    //     forward_->createJacobian(model_);
-    //     if (verbose_) std::cout << swatch.duration(true) << " s" << std::endl;
-    // }
+        //     Stopwatch swatch(true);
+        //     if (verbose_) std::cout << "Recalculating Jacobian matrix ...";
+        //     forward_->createJacobian(model_);
+        //     if (verbose_) std::cout << swatch.duration(true) << " s" << std::endl;
+        // }
 
     if (!localRegularization_) {
         DOSAVE echoMinMax(model_, "model: ");
@@ -304,10 +303,7 @@ bool RInversion::oneStep() {
         DOSAVE save(tM_->deriv(model_), "modelTrans");
         DOSAVE save(tD_->deriv(response_), "responseTrans");
 
-
-
         if (verbose_) std::cout << "solve CGLSCDWWtrans with lambda = " << lambda_ << std::endl;
-
         try{
             solveCGLSCDWWhtrans(*forward_->jacobian(), *forward_->constraints(),
                                 dataWeight_,
@@ -322,19 +318,18 @@ bool RInversion::oneStep() {
             // __MS("Debug halt! 1700")
             // forward_->mesh()->save("S1700.bms");
             // forward_->regionManager().mesh().save("R1700.bms");
-            
+
             // __MS("Debug halt! 17708")
             // forward_->mesh()->save("S17708.bms");
             // forward_->regionManager().mesh().save("R17708.bms");
-            
+
             // std::cout<<forward_->mesh()->boundary(2121).marker() << std::endl;
             // std::cout<<forward_->mesh()->boundary(2122).marker() << std::endl;
 
             // exit(1);
             log(Error, "solveCGLSCDWWhtrans throws unknown exception.");
         }
     } // else no optimization
-
     //restrictMax(deltaModelIter_, 50.0); // only work for log trans
     DOSAVE echoMinMax(deltaModelIter_, "dm");
     modelNew = tM_->update(model_, deltaModelIter_);
@@ -394,7 +389,6 @@ bool RInversion::oneStep() {
         DOSAVE forward_->mesh()->exportVTK("fop-model" + str(iter_));
         DOSAVE forward_->mesh()->clearData();
     }
-
     return true;
 }
 

core/src/meshgenerators.cpp

@@ -246,7 +246,7 @@ Mesh createMesh3D(const Mesh & mesh, const RVector & z, int topMarker, int botto
             }
         }
     }
-
+    mesh3.createNeighborInfos();
     return mesh3;
 }
 

core/src/ttdijkstramodelling.cpp

@@ -608,8 +608,8 @@ void TravelTimeDijkstraModelling::createJacobian(RSparseMapMatrix & jacobian,
     }
     if (this->verbose()){
         std::cout << "/" << swatch.duration(true) << " ";
+        std::cout << std::endl;
     }
-    std::cout << std::endl;
 }
 
 TTModellingWithOffset::TTModellingWithOffset(Mesh & mesh, DataContainer & dataContainer, bool verbose)

doc/examples/3_ert/plot_02_ert_field_data.py

@@ -9,8 +9,11 @@
 - 2D inversion with topography
 - geometric factor generation
 - topography effect
+
+The data is the profile 11 already shown by Günther et al. (2006, Fig. 11).
 """
-# sphinx_gallery_thumbnail_number = 5
+# sphinx_gallery_thumbnail_number = 7
+import matplotlib.pyplot as plt
 import pygimli as pg
 from pygimli.physics import ert
 
@@ -21,36 +24,46 @@
 data = pg.getExampleData('ert/slagdump.ohm', verbose=True)
 print(data)
 
+###############################################################################
+# Let us first have a look at the topography contained in the data
+plt.plot(pg.x(data), pg.z(data), 'x-')
+
 ###############################################################################
 # The data file does not contain geometric factors (token field 'k'),
 # so we create them based on the given topography.
+k0 = ert.createGeometricFactors(data)  # the analytical one
 data['k'] = ert.createGeometricFactors(data, numerical=True)
 
-###############################################################################
-# We initialize the ERTManager for further steps and eventually inversion.
-mgr = ert.ERTManager(sr=False)
-
 ###############################################################################
 # It might be interesting to see the topography effect, i.e the ratio between
-# the numerically computed geometry factor and the analytical formula
-k0 = ert.createGeometricFactors(data)
-_ = ert.showData(data, vals=k0/data['k'], label='Topography effect')
+# the numerically computed geometry factor and the analytical formula after
+# Rücker et al. (2006). We display it using a colormap with neutral white.
+_ = ert.showData(data, vals=k0/ data['k'], label='Topography effect',
+                 cMin=2/3, cMax=3/2, logScale=True, cMap="bwr")
 
 ###############################################################################
-# The data container has no apparent resistivities (token field 'rhoa') yet.
-# We can let the Manager fix this later for us (as we now have the 'k' field),
-# or we do it manually.
-mgr.checkData(data)
-print(data)
+# We can now compute the apparent resistivity and display it, once with the
+# wrong analytical formula and once with the numerical values in data['k']
+data['rhoa'] = data['r'] * data['k']
+kw = dict(cMin=6, cMax=33)
+fig, ax = plt.subplots(ncols=2)
+data.show(data['r']*k0, ax=ax[0], **kw);
+data.show(ax=ax[1], **kw)
+ax[0].set_title('Uncorrected')
+ax[1].set_title('Corrected');
 
 ###############################################################################
 # The data container does not necessarily contain data errors data errors
 # (token field 'err'), requiring us to enter data errors. We can let the
 # manager guess some defaults for us automaticly or set them manually
-data['err'] = ert.estimateError(data, relativeError=0.03, absoluteUError=5e-5)
-# or manually:
-# data['err'] = data_errors  # somehow
-_ = ert.show(data, data['err']*100)
+data.estimateError(relativeError=0.03, absoluteUError=5e-5)
+# which internally calls
+# data['err'] = ert.estimateError(data, ...)  # can also set manually
+_ = data.show(data['err']*100, label='error estimate (%)')
+
+###############################################################################
+# We initialize the ERTManager for further steps and eventually inversion.
+mgr = ert.ERTManager(data)
 
 ###############################################################################
 # Now the data have all necessary fields ('rhoa', 'err' and 'k') so we can run
@@ -59,8 +72,7 @@
 #
 mod = mgr.invert(data, lam=10, verbose=True,
                  paraDX=0.3, paraMaxCellSize=10, paraDepth=20, quality=33.6)
-
-_ = mgr.showResult()
+ax, cb = mgr.showResult()
 
 ###############################################################################
 # We can view the resulting model in the usual way.
@@ -69,4 +81,4 @@
 
 ###############################################################################
 # Or just plot the model only using your own options.
-_ = mgr.showResult(mod, cMin=5, cMax=30, cMap="Spectral_r", logScale=True)
+ax, cb = mgr.showResult(mod, cMin=5, cMax=30, cMap="Spectral_r", logScale=True)

doc/tutorials/2_modelling/plot_3-mod-fem-BC.py

@@ -72,9 +72,9 @@ def uDirichlet(boundary):
 # Alternatively we can define the gradients of the solution on the boundary,
 # i.e., Neumann type BC. This is done with another dictionary {marker: value} and passed by the bc dictionary.
 neumannBC = {1: -0.5,  # left
-             4: 2.5}  # bottom
+             3: 2.5}  # bottom
 
-dirichletBC = {3: 1.0}  # top
+dirichletBC = {4: 1.0}  # top
 
 u = solve(grid, f=0., bc={'Dirichlet': dirichletBC, 'Neumann': neumannBC})
 

doc/tutorials/2_modelling/plot_5-mod-fem-heat-2d.py

@@ -38,7 +38,7 @@
 ###############################################################################
 # Call :py:func:`pygimli.solver.solveFiniteElements` to solve the heat
 # diffusion equation :math:`\nabla\cdot(a\nabla T)=0` with :math:`T(bottom)=0`
-# (boundary marker 8) and :math:`T(top)=1` (boundary marker 4), where :math:`a`
+# (boundary marker 4) and :math:`T(top)=1` (boundary marker 8), where :math:`a`
 # is the thermal diffusivity and :math:`T` is the temperature distribution.
 # We assign thermal diffusivities to the four # regions using their marker
 # numbers in a dictionary (a) and the fixed temperatures at the boundaries

pygimli/CMakeLists.txt

@@ -49,7 +49,7 @@ else()
     add_custom_target(whl${TARGET_NAME}PackageBuild DEPENDS whlpgcore
         COMMAND "${Python_EXECUTABLE}" -m pip wheel . --wheel-dir=${WHEELHOUSE} --find-links ${WHEELHOUSE}
         WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}"
-        COMMENT "Building python wheel package for ${TARGET_NAME}"
+        COMMENT "Building python wheel package (${WHEELFILE}) for ${TARGET_NAME}"
     )
 
     add_custom_target(whl${TARGET_NAME} DEPENDS whl${TARGET_NAME}PackageBuild)
@@ -59,12 +59,12 @@ else()
     add_custom_target(whl${TARGET_NAME}TestInstall DEPENDS whl${TARGET_NAME} DEPENDS whlpgcoreTestRun
         COMMAND "${Python3_EXECUTABLE}" -m pip uninstall -y ${TARGET_NAME}
         COMMAND "${Python3_EXECUTABLE}" -m pip install ${WHEELHOUSE}/${PY_WHEELFILE} --quiet --find-links ${WHEELHOUSE}
-        COMMENT "Installing ${TARGET_NAME} in virtual test environment"
+        COMMENT "Installing ${TARGET_NAME} ((${WHEELFILE}) in virtual test environment"
     )
     add_custom_target(whl${TARGET_NAME}TestRun DEPENDS whl${TARGET_NAME}TestInstall
         COMMAND ${Python3_EXECUTABLE} -c "import pygimli as pg; pg.version()"
         VERBATIM
-        COMMENT "Testing ${TARGET_NAME} installation in virtual test environment. ${Python3_EXECUTABLE} -c \'import pygimli as pg; pg.version()\'"
+        COMMENT "Testing ${TARGET_NAME} installation (${WHEELFILE}) in virtual test environment. ${Python3_EXECUTABLE} -c \'import pygimli as pg; pg.version()\'"
     )
     add_custom_target(whl${TARGET_NAME}Test DEPENDS whl${TARGET_NAME}TestRun)
 

pygimli/core/datacontainer.py

@@ -220,3 +220,28 @@ def __DataContainer_subset(self, **kwargs):
     return new
 
 DataContainer.subset = __DataContainer_subset
+
+def __DataContainer_markSensorInvalid(self, idx):
+    """Mark all measurements using a specific sensor invalid."""
+    for tok in self.dataMap().keys():
+        if self.isSensorIndex(tok):
+            self.markInvalid(self[tok] == idx)
+
+DataContainer.markSensorInvalid = __DataContainer_markSensorInvalid
+
+def __DataContainer_removeSensorData(self, idx):
+    """Remove all measurements using a specific sensor."""
+    self.markSensorInvalid(idx)
+    self.removeInvalid()  # TODO4Nico: do it correctly
+    ## It will also remove any previously invalid measurements
+    ## So you better use only remove or markValid
+
+DataContainer.removeSensorData = __DataContainer_removeSensorData
+
+def __DataContainer_removeSensor(self, idx):
+    """Remove a specific sensor."""
+    self.removeSensorData(idx)
+    self.removeUnusedSensors()
+    ## Same argument as above, there might be other unused ones
+
+DataContainer.removeSensor = __DataContainer_removeSensor