Merge branch 'NIFTy_7' into flexible_FieldZeroPadder

.gitlab-ci.yml

@@ -5,11 +5,20 @@ variables:
   OMP_NUM_THREADS: 1
 
 stages:
+  - static_checks
   - build_docker
   - test
   - release
   - demo_runs
 
+check_no_asserts:
+  image: debian:testing-slim
+  stage: static_checks
+  before_script:
+    - ls
+  script:
+    - if [ `grep -r "^[[:space:]]*assert[ (]" src demos | wc -l` -ne 0 ]; then echo "Have found assert statements. Don't use them! Use \`utilities.myassert\` instead." && exit 1; fi
+
 build_docker_from_scratch:
   only:
     - schedules
@@ -110,6 +119,14 @@ run_getting_started_mf:
     paths:
       - '*.png'
 
+run_getting_density:
+  stage: demo_runs
+  script:
+    - python3 demos/getting_started_density.py
+  artifacts:
+    paths:
+      - '*.png'
+
 run_bernoulli:
   stage: demo_runs
   script:
@@ -126,7 +143,7 @@ run_curve_fitting:
     paths:
       - '*.png'
 
-run_visual_mgvi:
+run_visual_vi:
   stage: demo_runs
   script:
-    - python3 demos/mgvi_visualized.py
+    - python3 demos/variational_inference_visualized.py

ChangeLog.md

 Changes since NIFTy 6
 =====================
 
+New parametric amplitude model
+------------------------------
+
+The `ift.CorrelatedFieldMaker` now features two amplitude models. In addition
+to the non-parametric one, one may choose to use a Matern kernel instead. The
+method is aptly named `add_fluctuations_matern`. The major advantage of the
+parametric model is its more intuitive scaling with the size of the position
+space.
+
 CorrelatedFieldMaker interface change
 -------------------------------------
 
-The interface of `ift.CorrelatedFieldMaker.make` and
-`ift.CorrelatedFieldMaker.add_fluctuations` changed; it now expects the mean
-and the standard deviation of their various parameters not as separate
-arguments but as a tuple.
+The interface of `ift.CorrelatedFieldMaker` changed and instances of it may now
+be instantiated directly without the previously required `make` method. Upon
+initialization, no zero-mode must be specified as the normalization for the
+different axes of the power respectively amplitude spectrum now only happens
+once in the `finalize` method. There is now a new call named
+`set_amplitude_total_offset` to set the zero-mode. The method accepts either an
+instance of `ift.Operator` or a tuple parameterizing a log-normal parameter.
+Methods which require the zero-mode to be set raise a `NotImplementedError` if
+invoked prior to having specified a zero-mode.
+
+Furthermore, the interface of `ift.CorrelatedFieldMaker.add_fluctuations`
+changed; it now expects the mean and the standard deviation of their various
+parameters not as separate arguments but as a tuple. The same applies to all
+new and renamed methods of the `CorrelatedFieldMaker` class.
 
 Furthermore, it is now possible to disable the asperity and the flexibility
 together with the asperity in the correlated field model. Note that disabling
@@ -45,13 +64,32 @@ The implementation tests for nonlinear operators are now available in
 `ift.extra.check_operator()` and for linear operators
 `ift.extra.check_linear_operator()`.
 
-
 MetricGaussianKL interface
 --------------------------
 
-Users do not instantiate `MetricGaussianKL` by its constructor anymore. Rather
-`MetricGaussianKL.make()` shall be used. Additionally, `mirror_samples` is not
-set by default anymore.
+`mirror_samples` is not set by default anymore.
+
+
+GeoMetricKL
+-----------
+
+A new posterior approximation scheme, called geometric Variational Inference
+(geoVI) was introduced. `GeoMetricKL` extends `MetricGaussianKL` in the sense
+that it uses (non-linear) geoVI samples instead of (linear) MGVI samples.
+`GeoMetricKL` can be configured such that it reduces to `MetricGaussianKL`.
+`GeoMetricKL` is now used in `demos/getting_started_3.py` and a visual
+comparison to MGVI can be found in `demos/variational_inference_visualized.py`.
+For further details see (<https://arxiv.org/abs/2105.10470>).
+
+
+LikelihoodOperator
+------------------
+
+A new subclass of `EnergyOperator` was introduced and all `EnergyOperator`s
+that are likelihoods are now `LikelihoodOperator`s. A `LikelihoodOperator`
+has to implement the function `get_transformation`, which returns a
+coordinate transformation in which the Fisher metric of the likelihood becomes
+the identity matrix. This is needed for the `GeoMetricKL` algorithm.
 
 
 Changes since NIFTy 5

README.md

@@ -134,17 +134,32 @@ The NIFTy package is licensed under the terms of the
 *without any warranty*.
 
 
+Contributing
+------------
+
+Please note our convention not to use pure Python `assert` statements in
+production code. They are not guaranteed to by executed by Python and can be
+turned off by the user (`python -O` in cPython). As an alternative use
+`ift.myassert`.
+
+
 Contributors
 ------------
 
 ### NIFTy7
 
+- Andrija Kostic
 - Gordian Edenhofer
+- Jakob Knollmüller
+- Jakob Roth
 - Lukas Platz
+- Matteo Guardiani
 - Martin Reinecke
-- [Philipp Arras](https://wwwmpa.mpa-garching.mpg.de/~parras/)
+- [Philipp Arras](https://philipp-arras.de)
 - Philipp Frank
 - [Reimar Heinrich Leike](https://wwwmpa.mpa-garching.mpg.de/~reimar/)
+- Simon Ding
+- Vincent Eberle
 
 ### NIFTy6
 

demos/getting_started_3.py

@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -31,6 +31,7 @@ import numpy as np
 
 import nifty7 as ift
 
+ift.random.push_sseq_from_seed(27)
 
 def random_los(n_los):
     starts = list(ift.random.current_rng().random((n_los, 2)).T)
@@ -63,16 +64,16 @@ def main():
         'offset_std': (1e-3, 1e-6),
 
         # Amplitude of field fluctuations
-        'fluctuations': (2., 1.),  # 1.0, 1e-2
+        'fluctuations': (1., 0.8),  # 1.0, 1e-2
 
         # Exponent of power law power spectrum component
-        'loglogavgslope': (-4., 1),  # -6.0, 1
+        'loglogavgslope': (-3., 1),  # -6.0, 1
 
         # Amplitude of integrated Wiener process power spectrum component
-        'flexibility': (5, 2.),  # 2.0, 1.0
+        'flexibility': (2, 1.),  # 1.0, 0.5
 
         # How ragged the integrated Wiener process component is
-        'asperity': (0.5, 0.5)  # 0.1, 0.5
+        'asperity': (0.5, 0.4)  # 0.1, 0.5
     }
 
     correlated_field = ift.SimpleCorrelatedField(position_space, **args)
@@ -89,7 +90,6 @@ def main():
     # Specify noise
     data_space = R.target
     noise = .001
-    sig = ift.ScalingOperator(data_space, np.sqrt(noise))
     N = ift.ScalingOperator(data_space, noise)
 
     # Generate mock signal and data
@@ -97,11 +97,14 @@ def main():
     data = signal_response(mock_position) + N.draw_sample_with_dtype(dtype=np.float64)
 
     # Minimization parameters
-    ic_sampling = ift.AbsDeltaEnergyController(
-        deltaE=0.05, iteration_limit=100)
-    ic_newton = ift.AbsDeltaEnergyController(
-        name='Newton', deltaE=0.5, iteration_limit=35)
+    ic_sampling = ift.AbsDeltaEnergyController(name="Sampling (linear)",
+            deltaE=0.05, iteration_limit=100)
+    ic_newton = ift.AbsDeltaEnergyController(name='Newton', deltaE=0.5,
+            convergence_level=2, iteration_limit=35)
     minimizer = ift.NewtonCG(ic_newton)
+    ic_sampling_nl = ift.AbsDeltaEnergyController(name='Sampling (nonlin)',
+            deltaE=0.5, iteration_limit=15, convergence_level=2)
+    minimizer_sampling = ift.NewtonCG(ic_sampling_nl)
 
     # Set up likelihood and information Hamiltonian
     likelihood = (ift.GaussianEnergy(mean=data, inverse_covariance=N.inverse) @
@@ -112,18 +115,21 @@ def main():
     mean = initial_mean
 
     plot = ift.Plot()
-    plot.add(signal(mock_position), title='Ground Truth')
+    plot.add(signal(mock_position), title='Ground Truth', zmin = 0, zmax = 1)
     plot.add(R.adjoint_times(data), title='Data')
     plot.add([pspec.force(mock_position)], title='Power Spectrum')
     plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename.format("setup"))
 
     # number of samples used to estimate the KL
-    N_samples = 20
+    N_samples = 10
 
-    # Draw new samples to approximate the KL five times
-    for i in range(5):
+    # Draw new samples to approximate the KL six times
+    for i in range(6):
+        if i==5:
+            # Double the number of samples in the last step for better statistics
+            N_samples = 2*N_samples
         # Draw new samples and minimize KL
-        KL = ift.MetricGaussianKL.make(mean, H, N_samples, True)
+        KL = ift.GeoMetricKL(mean, H, N_samples, minimizer_sampling, True)
         KL, convergence = minimizer(KL)
         mean = KL.position
         ift.extra.minisanity(data, lambda x: N.inverse, signal_response,
@@ -131,7 +137,7 @@ def main():
 
         # Plot current reconstruction
         plot = ift.Plot()
-        plot.add(signal(KL.position), title="Latent mean")
+        plot.add(signal(KL.position), title="Latent mean", zmin = 0, zmax = 1)
         plot.add([pspec.force(KL.position + ss) for ss in KL.samples],
                  title="Samples power spectrum")
         plot.output(ny=1, ysize=6, xsize=16,
@@ -144,16 +150,16 @@ def main():
     # Plotting
     filename_res = filename.format("results")
     plot = ift.Plot()
-    plot.add(sc.mean, title="Posterior Mean")
+    plot.add(sc.mean, title="Posterior Mean", zmin = 0, zmax = 1)
     plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
 
     powers = [pspec.force(s + KL.position) for s in KL.samples]
     sc = ift.StatCalculator()
     for pp in powers:
-        sc.add(pp)
+        sc.add(pp.log())
     plot.add(
         powers + [pspec.force(mock_position),
-                  pspec.force(KL.position), sc.mean],
+                  pspec.force(KL.position), sc.mean.exp()],
         title="Sampled Posterior Power Spectrum",
         linewidth=[1.]*len(powers) + [3., 3., 3.],
         label=[None]*len(powers) + ['Ground truth', 'Posterior latent mean', 'Posterior mean'])

demos/getting_started_4_CorrelatedFields.ipynb

@@ -87,7 +87,8 @@
     "main_sample = None\n",
     "def init_model(m_pars, fl_pars):\n",
     "    global main_sample\n",
-    "    cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
+    "    cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
+    "    cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
     "    cf.add_fluctuations(**fl_pars)\n",
     "    field = cf.finalize(prior_info=0)\n",
     "    main_sample = ift.from_random(field.domain)\n",
@@ -95,7 +96,8 @@
     "    \n",
     "# Helper: field and spectrum from parameter dictionaries + plotting\n",
     "def eval_model(m_pars, fl_pars, title=None, samples=None):\n",
-    "    cf = ift.CorrelatedFieldMaker.make(**m_pars)\n",
+    "    cf = ift.CorrelatedFieldMaker(m_pars[\"prefix\"])\n",
+    "    cf.set_amplitude_total_offset(m_pars[\"offset_mean\"], m_pars[\"offset_std\"])\n",
     "    cf.add_fluctuations(**fl_pars)\n",
     "    field = cf.finalize(prior_info=0)\n",
     "    spectrum = cf.amplitude\n",
@@ -473,7 +475,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.9.2"
   }
  },
  "nbformat": 4,

demos/getting_started_5_mf.py

@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2020 Max-Planck-Society
+# Copyright(C) 2013-2021 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
@@ -73,14 +73,17 @@ def main():
     sp2 = ift.RGSpace(npix2)
 
     # Set up signal model
-    cfmaker = ift.CorrelatedFieldMaker.make(0., (1e-2, 1e-6), '')
-    cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (2, .2), (.01, .5), (-4, 2.), 'amp1')
-    cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (2, .2), (.01, .5),
-                             (-3, 1), 'amp2')
+    cfmaker = ift.CorrelatedFieldMaker('')
+    cfmaker.add_fluctuations(sp1, (0.1, 1e-2), (2, .2), (.01, .5), (-4, 2.),
+                             'amp1')
+    cfmaker.add_fluctuations(sp2, (0.1, 1e-2), (2, .2), (.01, .5), (-3, 1),
+                             'amp2')
+    cfmaker.set_amplitude_total_offset(0., (1e-2, 1e-6))
     correlated_field = cfmaker.finalize()
 
-    pspec1 = cfmaker.normalized_amplitudes[0]**2
-    pspec2 = cfmaker.normalized_amplitudes[1]**2
+    normalized_amp = cfmaker.get_normalized_amplitudes()
+    pspec1 = normalized_amp[0]**2
+    pspec2 = normalized_amp[1]**2
     DC = SingleDomain(correlated_field.target, position_space)
 
     # Apply a nonlinearity
@@ -131,7 +134,7 @@ def main():
 
     for i in range(5):
         # Draw new samples and minimize KL
-        KL = ift.MetricGaussianKL.make(mean, H, N_samples, True)
+        KL = ift.MetricGaussianKL(mean, H, N_samples, True)
         KL, convergence = minimizer(KL)
         mean = KL.position
 

demos/getting_started_density.py

+#!/usr/bin/env python3
+
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2021 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+############################################################
+# Density estimation
+#
+# Compute a density estimate for a log-normal process measured by a
+# Poissonian likelihood.
+#
+# Demo takes a while to compute
+#############################################################
+
+import numpy as np
+
+import nifty7 as ift
+
+if __name__ == "__main__":
+    filename = "getting_started_density_{}.png"
+    ift.random.push_sseq_from_seed(42)
+
+    # Set up signal domain
+    npix1 = 128
+    npix2 = 128
+    sp1 = ift.RGSpace(npix1)
+    sp2 = ift.RGSpace(npix2)
+    position_space = ift.DomainTuple.make((sp1, sp2))
+
+    signal, ops = ift.density_estimator(position_space)
+    correlated_field = ops["correlated_field"]
+
+    data_space = signal.target
+
+    # Generate mock signal and data
+    rng = ift.random.current_rng()
+    mock_position = ift.from_random(signal.domain, "normal")
+    data = ift.Field.from_raw(data_space, rng.poisson(signal(mock_position).val))
+
+    # Rejoin domains for plotting
+    plotting_domain = ift.DomainTuple.make(ift.RGSpace((npix1, npix2)))
+    plotting_domain_expanded = ift.DomainTuple.make(ift.RGSpace((2 * npix1, 2 * npix2)))
+
+    plot = ift.Plot()
+    plot.add(
+        ift.Field.from_raw(
+            plotting_domain_expanded, ift.exp(correlated_field(mock_position)).val
+        ),
+        title="Pre-Slicing Truth",
+    )
+    plot.add(
+        ift.Field.from_raw(plotting_domain, signal(mock_position).val),
+        title="Ground Truth",
+    )
+    plot.add(ift.Field.from_raw(plotting_domain, data.val), title="Data")
+    plot.output(ny=1, nx=3, xsize=10, ysize=3, name=filename.format("setup"))
+    print("Setup saved as", filename.format("setup"))
+
+    # Minimization parameters
+    ic_sampling = ift.AbsDeltaEnergyController(
+        name="Sampling", deltaE=0.01, iteration_limit=100
+    )
+    ic_newton = ift.AbsDeltaEnergyController(
+        name="Newton", deltaE=0.01, iteration_limit=35
+    )
+    ic_sampling.enable_logging()
+    ic_newton.enable_logging()
+    minimizer = ift.NewtonCG(ic_newton, enable_logging=True)
+
+    # Number of samples used to estimate the KL
+    n_samples = 5
+
+    # Set up likelihood and information Hamiltonian
+    likelihood = ift.PoissonianEnergy(data) @ signal
+    ham = ift.StandardHamiltonian(likelihood, ic_sampling)
+
+    # Start minimization
+    initial_mean = ift.MultiField.full(ham.domain, 0.)
+    mean = initial_mean
+
+    for i in range(5):
+        # Draw new samples and minimize KL
+        kl = ift.MetricGaussianKL(mean, ham, n_samples, True)
+        kl, convergence = minimizer(kl)
+        mean = kl.position
+
+        # Plot current reconstruction
+        plot = ift.Plot()
+        plot.add(
+            ift.Field.from_raw(
+                plotting_domain_expanded, ift.exp(correlated_field(mock_position)).val
+            ),
+            title="Ground truth",
+        )
+        plot.add(
+            ift.Field.from_raw(plotting_domain, signal(mock_position).val),
+            title="Ground truth",
+        )
+        plot.add(
+            ift.Field.from_raw(plotting_domain, signal(kl.position).val),
+            title="Reconstruction",
+        )
+        plot.add(
+            (ic_newton.history, ic_sampling.history, minimizer.inversion_history),
+            label=["kl", "Sampling", "Newton inversion"],
+            title="Cumulative energies",
+            s=[None, None, 1],
+            alpha=[None, 0.2, None],
+        )
+        plot.output(
+            nx=3, ny=2, ysize=10, xsize=15, name=filename.format(f"loop_{i:02d}")
+        )
+
+    # Done, draw posterior samples
+    sc = ift.StatCalculator()
+    sc_unsliced = ift.StatCalculator()
+    for sample in kl.samples:
+        sc.add(signal(sample + kl.position))
+        sc_unsliced.add(ift.exp(correlated_field(sample + kl.position)))
+
+    # Plotting
+    plot = ift.Plot()
+    plot.add(ift.Field.from_raw(plotting_domain, sc.mean.val), title="Posterior Mean")
+    plot.add(
+        ift.Field.from_raw(plotting_domain, ift.sqrt(sc.var).val),
+        title="Posterior Standard Deviation",
+    )
+    plot.add(
+        ift.Field.from_raw(plotting_domain_expanded, sc_unsliced.mean.val),
+        title="Posterior Unsliced Mean",
+    )
+    plot.add(
+        ift.Field.from_raw(plotting_domain_expanded, ift.sqrt(sc_unsliced.var).val),
+        title="Posterior Unsliced Standard Deviation",
+    )
+    plot.output(xsize=15, ysize=15, name=filename.format("results"))
+    print("Saved results as", filename.format("results"))

demos/mgvi_visualized.py → demos/variational_inference_visualized.py

@@ -11,21 +11,21 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2020 Max-Planck-Society
-# Authors: Reimar Leike, Philipp Arras
+# Copyright(C) 2013-2021 Max-Planck-Society
+# Authors: Reimar Leike, Philipp Arras, Philipp Frank
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
 ###############################################################################
-# Metric Gaussian Variational Inference (MGVI)
+# Variational Inference (VI)
 #
-# This script demonstrates how MGVI works for an inference problem with only
-# two real quantities of interest. This enables us to plot the posterior
-# probability density as two-dimensional plot. The posterior samples generated
-# by MGVI are contrasted with the maximum-a-posterior (MAP) solution together
-# with samples drawn with the Laplace method. This method uses the local
-# curvature at the MAP solution as inverse covariance of a Gaussian probability
-# density.
+# This script demonstrates how MGVI and GeoVI work for an inference problem
+# with only two real quantities of interest. This enables us to plot the
+# posterior probability density as two-dimensional plot. The approximate
+# posterior samples are contrasted with the maximum-a-posterior (MAP) solution
+# together with samples drawn with the Laplace method. This method uses the
+# local curvature at the MAP solution as inverse covariance of a Gaussian
+# probability density.
 ###############################################################################
 
 import numpy as np
@@ -87,36 +87,52 @@ def main():
 
     minimizer = ift.NewtonCG(
         ift.GradientNormController(iteration_limit=2, name='Mini'))
-    pos = ift.from_random(ham.domain, 'normal')
-    plt.figure(figsize=[12, 8])
+    pos = pos1 = ift.from_random(ham.domain, 'normal')
+    fig, axs = plt.subplots(2, 1, figsize=[12, 8])
     for ii in range(15):
         if ii % 3 == 0:
-            mgkl = ift.MetricGaussianKL.make(pos, ham, 40, False)
-
-        plt.cla()
-        plt.imshow(z.T, origin='lower', norm=LogNorm(), vmin=1e-3,
-                   vmax=np.max(z), cmap='gist_earth_r',
-                   extent=x_limits_scaled + y_limits)
-        if ii == 0:
-            cbar = plt.colorbar()
-        cbar.ax.set_ylabel('pdf')
-        xs, ys = [], []
-        for samp in mgkl.samples:
-            samp = (samp + pos).val
-            xs.append(samp['a'])
-            ys.append(samp['b'])
-        plt.scatter(np.array(xs)*scale, np.array(ys), label='MGVI samples')
-        plt.scatter(pos.val['a']*scale, pos.val['b'], label='MGVI latent mean')
-        plt.scatter(np.array(map_xs)*scale, np.array(map_ys),
-                    label='Laplace samples')
-        plt.scatter(MAP.position.val['a']*scale, MAP.position.val['b'],
-                    label='Maximum a posterior solution')
-        plt.legend()
+            # Resample
+            mgkl = ift.MetricGaussianKL(pos, ham, 100, False)
+            mini_samp = ift.NewtonCG(ift.GradientNormController(iteration_limit=5))
+            geokl = ift.GeoMetricKL(pos1, ham, 100, mini_samp, False)
+
+        for axx in axs:
+            axx.clear()
+            im = axx.imshow(z.T, origin='lower', norm=LogNorm(vmin=1e-3, vmax=np.max(z)),
+                               cmap='gist_earth_r', extent=x_limits_scaled + y_limits)
+            if ii == 0:
+                cbar = plt.colorbar(im, ax=axx)
+                cbar.ax.set_ylabel('pdf')
+
+        for jj, nn, kl, pp in ((0, "MGVI", mgkl, pos), (1, "GeoVI", geokl, pos1)):
+            xs, ys = [], []
+            for samp in kl.samples:
+                samp = (samp + pp).val
+                xs.append(samp['a'])
+                ys.append(samp['b'])
+            axs[jj].scatter(np.array(xs)*scale, np.array(ys), label=f'{nn} samples')
+            axs[jj].scatter(pp.val['a']*scale, pp.val['b'], label=f'{nn} latent mean')
+            axs[jj].set_title(nn)
+
+        for axx in axs:
+            axx.scatter(np.array(map_xs)*scale, np.array(map_ys),
+                        label='Laplace samples')
+            axx.scatter(MAP.position.val['a']*scale, MAP.position.val['b'],
+                        label='Maximum a posterior solution')
+            axx.set_xlim(x_limits_scaled)
+            axx.set_ylim(y_limits)
+            axx.set_ylabel('y')
+            axx.legend(loc='lower right')
+        axs[0].xaxis.set_visible(False)
+        axs[1].set_xlabel('x')
+        plt.tight_layout()
         plt.draw()
         plt.pause(1.0)
 
         mgkl, _ = minimizer(mgkl)
+        geokl, _ = minimizer(geokl)
         pos = mgkl.position
+        pos1 = geokl.position
     ift.logger.info('Finished')
     # Uncomment the following line in order to leave the plots open
     # plt.show()