WIP: Advi and hmc

demos/ADVI_demo.py

+# 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-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+############################################################
+# Non-linear tomography
+#
+# The signal is a sigmoid-normal distributed field.
+# The data is the field integrated along lines of sight that are
+# randomly (set mode=0) or radially (mode=1) distributed
+#
+# Demo takes a while to compute
+#############################################################
+
+import numpy as np
+import nifty5 as ift
+from matplotlib import pyplot as plt
+
+if __name__ == '__main__':
+    np.random.seed(42)
+
+    mode = 'fullcovariance' # meanfield or fullcovariance
+
+    filename = "getting_started_3_mode_{}_".format(mode) + "{}.png"
+
+    position_space = ift.RGSpace([128])
+    harmonic_space = position_space.get_default_codomain()
+    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
+    power_space = ift.PowerSpace(harmonic_space)
+
+    # Set up an amplitude operator for the field
+    dct = {
+        'target': power_space,
+        'n_pix': 64,  # 64 spectral bins
+
+        # Spectral smoothness (affects Gaussian process part)
+        'a': 3,  # relatively high variance of spectral curbvature
+        'k0': .4,  # quefrency mode below which cepstrum flattens
+
+        # Power-law part of spectrum:
+        'sm': -3,  # preferred power-law slope
+        'sv': .5,  # low variance of power-law slope
+        'im':  0,  # y-intercept mean, in-/decrease for more/less contrast
+        'iv': .3   # y-intercept variance
+    }
+    A = ift.SLAmplitude(**dct)
+
+    # Build the operator for a correlated signal
+    power_distributor = ift.PowerDistributor(harmonic_space, power_space)
+    vol = harmonic_space.scalar_dvol**-0.5
+    xi = ift.ducktape(harmonic_space, None, 'xi')
+    correlated_field = ht(vol*power_distributor(A)*xi)
+    # Alternatively, one can use:
+    # correlated_field = ift.CorrelatedField(position_space, A)
+
+    # Apply a nonlinearity
+    signal = ift.sigmoid(correlated_field)
+
+    # Build the line-of-sight response and define signal response
+    R = ift.GeometryRemover(signal.target)
+    signal_response = R(signal)
+
+    # Specify noise
+    data_space = R.target
+    noise = .1
+    N = ift.ScalingOperator(noise, data_space)
+
+    # Generate mock signal and data
+    mock_position = ift.from_random('normal', signal_response.domain)
+    data = signal_response(mock_position) + N.draw_sample()
+
+    # Minimization parameters
+    ic_sampling = ift.GradientNormController(iteration_limit=100)
+    ic_newton = ift.GradInfNormController(
+        name='Newton', tol=1e-7, iteration_limit=5)
+    minimizer = ift.NewtonCG(ic_newton)
+
+    # Set up likelihood and information Hamiltonian
+    likelihood = ift.GaussianEnergy(mean=data,
+                                    inverse_covariance=N.inverse)(signal_response)
+    H = ift.StandardHamiltonian(likelihood, ic_sampling)
+
+
+    if mode == 'meanfield':
+        (model, entropy, para, var, mean) = ift.build_meanfield(H)
+    if mode == 'fullcovariance':
+        (model, entropy, para, var, mean) = ift.build_fullcovariance(H)
+
+    # Draw new samples to approximate the KL five times
+    N_samp = 10
+    for i in range(50):
+        # Draw new samples and minimize KL
+        KL = ift.ParametrizedGaussianKL(para, H, model,
+                                             entropy, N_samp)
+        KL, convergence = minimizer(KL)
+        para = KL.position
+
+        plt.close('all')
+        plt.figure()
+        plt.subplot(211)
+        plt.plot(signal(mock_position).val, 'r-', alpha=1)
+        plt.plot(data.val, 'kx', alpha=1)
+
+        for sample in KL.samples:
+            plt.plot(signal(model(KL.position + sample)).val, 'k-', alpha=0.2)
+
+        plt.subplot(212)
+        plt.plot(A.force(mock_position).val, 'r-', alpha=1)
+
+        for sample in KL.samples:
+            plt.plot(A.force(model(KL.position + sample)).val, 'k-', alpha=0.2)
+        plt.yscale('log')
+        plt.xscale('log')
+        plt.pause(0.01)
+
+
+
+
+
+
+
+
+
+