Merge branch 'demos' into 'NIFTy_5'

new demos

See merge request ift/nifty-dev!11

.gitlab-ci.yml

@@ -63,133 +63,6 @@ pages:
 before_script:
   - export MPLBACKEND="agg"
 
-run_critical_filtering:
-  stage: demo_runs
-  script:
-    - ls
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/critical_filtering.py
-    - python3 demos/critical_filtering.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_nonlinear_critical_filter:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/nonlinear_critical_filter.py
-    - python3 demos/nonlinear_critical_filter.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_nonlinear_wiener_filter:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/nonlinear_wiener_filter.py
-    - python3 demos/nonlinear_wiener_filter.py
-  only:
-    - run_demos
-  artifacts:
-    paths:
-      - '*.png'
-
-# FIXME: disable for now. Fixing it is part of issue #244.
-#run_poisson_demo:
-#  stage: demo_runs
-#  script:
-#    - python setup.py install --user -f
-#    - python3 setup.py install --user -f
-#    - python demos/poisson_demo.py
-#    - python3 demos/poisson_demo.py
-#  artifacts:
-#    paths:
-#      - '*.png'
-
-run_probing:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/probing.py
-    - python3 demos/probing.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_sampling:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/sampling.py
-    - python3 demos/sampling.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_tomography:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/tomography.py
-    - python3 demos/tomography.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_wiener_filter_data_space_noiseless:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/wiener_filter_data_space_noiseless.py
-    - python3 demos/wiener_filter_data_space_noiseless.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_wiener_filter_easy.py:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/wiener_filter_easy.py
-    - python3 demos/wiener_filter_easy.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_wiener_filter_via_curvature.py:
-  stage: demo_runs
-  script:
-    - pip install --user numericalunits
-    - pip3 install --user numericalunits
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/wiener_filter_via_curvature.py
-    - python3 demos/wiener_filter_via_curvature.py
-  artifacts:
-    paths:
-      - '*.png'
-
-run_wiener_filter_via_hamiltonian.py:
-  stage: demo_runs
-  script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
-    - python demos/wiener_filter_via_hamiltonian.py
-    - python3 demos/wiener_filter_via_hamiltonian.py
-  artifacts:
-    paths:
-      - '*.png'
-
 run_ipynb:
   stage: demo_runs
   script:

demos/Wiener_Filter.ipynb

@@ -717,21 +717,21 @@
  "metadata": {
   "celltoolbar": "Slideshow",
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "python2"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.15"
+   "pygments_lexer": "ipython3",
+   "version": "3.6.5"
   }
  },
  "nbformat": 4,

demos/critical_filtering.py

-import nifty5 as ift
-import numpy as np
-from nifty5 import Exponential, Linear, Tanh
-
-np.random.seed(42)
-
-
-def adjust_zero_mode(m0, t0):
-    mtmp = m0.to_global_data().copy()
-    zero_position = len(m0.shape)*(0,)
-    zero_mode = mtmp[zero_position]
-    mtmp[zero_position] = zero_mode / abs(zero_mode)
-    ttmp = t0.to_global_data().copy()
-    ttmp[0] += 2 * np.log(abs(zero_mode))
-    return (ift.Field.from_global_data(m0.domain, mtmp),
-            ift.Field.from_global_data(t0.domain, ttmp))
-
-
-if __name__ == "__main__":
-    noise_level = 1.
-    p_spec = (lambda k: (.5 / (k + 1) ** 3))
-
-    nonlinearity = Linear()
-    # Set up position space
-    s_space = ift.RGSpace((128, 128))
-    h_space = s_space.get_default_codomain()
-
-    # Define harmonic transformation and associated harmonic space
-    HT = ift.HarmonicTransformOperator(h_space, target=s_space)
-
-    # Setting up power space
-    p_space = ift.PowerSpace(h_space,
-                             binbounds=ift.PowerSpace.useful_binbounds(
-                                 h_space, logarithmic=True))
-
-    # Choosing the prior correlation structure and defining
-    # correlation operator
-    p = ift.PS_field(p_space, p_spec)
-    log_p = ift.log(p)
-    S = ift.create_power_operator(h_space, power_spectrum=lambda k: 1e-5)
-
-    # Drawing a sample sh from the prior distribution in harmonic space
-    sh = S.draw_sample()
-
-    # Choosing the measurement instrument
-    # Instrument = SmoothingOperator(s_space, sigma=0.01)
-    mask = np.ones(s_space.shape)
-    # mask[6000:8000] = 0.
-    mask[30:70, 30:70] = 0.
-    mask = ift.Field.from_global_data(s_space, mask)
-
-    MaskOperator = ift.DiagonalOperator(mask)
-    R = ift.GeometryRemover(s_space)
-    R = R*MaskOperator
-    # R = R*HT
-    # R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
-    #                                                response_sigma)
-    MeasurementOperator = R
-
-    d_space = MeasurementOperator.target
-
-    Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
-    power = Distributor(ift.exp(0.5*log_p))
-    # Creating the mock data
-    true_sky = nonlinearity(HT(power*sh))
-    noiseless_data = MeasurementOperator(true_sky)
-    noise_amplitude = noiseless_data.val.std()*noise_level
-    N = ift.ScalingOperator(noise_amplitude**2, d_space)
-    n = N.draw_sample()
-    # Creating the mock data
-    d = noiseless_data + n
-
-    m0 = ift.full(h_space, 1e-7)
-    t0 = ift.full(p_space, -4.)
-    power0 = Distributor.times(ift.exp(0.5 * t0))
-
-    plotdict = {"colormap": "Planck-like"}
-    zmin = true_sky.min()
-    zmax = true_sky.max()
-    ift.plot(true_sky, title="True sky", name="true_sky.png", **plotdict)
-    ift.plot(MeasurementOperator.adjoint_times(d), title="Data",
-             name="data.png", **plotdict)
-
-    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
-                                     tol_abs_gradnorm=1e-3)
-    LS = ift.LineSearchStrongWolfe(c2=0.02)
-    minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
-
-    IC = ift.GradientNormController(iteration_limit=500,
-                                    tol_abs_gradnorm=1e-3)
-
-    for i in range(20):
-        power0 = Distributor(ift.exp(0.5*t0))
-        map0_energy = ift.library.NonlinearWienerFilterEnergy(
-            m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
-
-        # Minimization with chosen minimizer
-        map0_energy, convergence = minimizer(map0_energy)
-        m0 = map0_energy.position
-
-        # Updating parameters for correlation structure reconstruction
-        D0 = map0_energy.curvature
-
-        # Initializing the power energy with updated parameters
-        power0_energy = ift.library.NonlinearPowerEnergy(
-            position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
-            Instrument=MeasurementOperator, nonlinearity=nonlinearity,
-            Distributor=Distributor, sigma=1., samples=2,
-            iteration_controller=IC)
-
-        power0_energy = minimizer(power0_energy)[0]
-
-        # Setting new power spectrum
-        t0 = power0_energy.position
-
-        # break degeneracy between amplitude and excitation by setting
-        # excitation monopole to 1
-        m0, t0 = adjust_zero_mode(m0, t0)
-
-    ift.plot(nonlinearity(HT(power0*m0)), title="Reconstructed sky",
-             name="reconstructed_sky.png", zmin=zmin, zmax=zmax, **plotdict)
-    ymin = np.min(p.to_global_data())
-    ift.plot([ift.exp(t0), p], title="Power spectra",
-             name="ps.png", ymin=ymin, **plotdict)

demos/getting_started_1.py

+import nifty5 as ift
+import numpy as np
+from global_newton.models_other.apply_data import ApplyData
+from global_newton.models_energy.hamiltonian import Hamiltonian
+from nifty5.library.unit_log_gauss import UnitLogGauss
+if __name__ == '__main__':
+    # s_space = ift.RGSpace([1024])
+    s_space = ift.RGSpace([128,128])
+    # s_space = ift.HPSpace(64)
+
+    h_space = s_space.get_default_codomain()
+    total_domain = ift.MultiDomain.make({'xi': h_space})
+    HT = ift.HarmonicTransformOperator(h_space, s_space)
+
+    def sqrtpspec(k):
+        return 16. / (20.+k**2)
+
+    GR = ift.GeometryRemover(s_space)
+
+    d_space = GR.target
+    B = ift.FFTSmoothingOperator(s_space,0.1)
+    mask = np.ones(s_space.shape)
+    mask[64:89,76:100] = 0.
+    mask = ift.Field(s_space,val=mask)
+    Mask = ift.DiagonalOperator(mask)
+    R = GR * Mask * B
+    noise = 1.
+    N = ift.ScalingOperator(noise, d_space)
+
+    p_space = ift.PowerSpace(h_space)
+    pd = ift.PowerDistributor(h_space, p_space)
+    position = ift.from_random('normal', total_domain)
+    xi = ift.Variable(position)['xi']
+    a = ift.Constant(position, ift.PS_field(p_space, sqrtpspec))
+    A = pd(a)
+    s_h = A * xi
+    s = HT(s_h)
+    Rs = R(s)
+
+
+
+    MOCK_POSITION = ift.from_random('normal',total_domain)
+    data = Rs.at(MOCK_POSITION).value + N.draw_sample()
+
+    NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs)
+
+    INITIAL_POSITION = ift.from_random('normal',total_domain)
+    likelihood = UnitLogGauss(INITIAL_POSITION, NWR)
+
+    IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
+    inverter = ift.ConjugateGradient(controller=IC)
+    IC2 = ift.GradientNormController(name='Newton', iteration_limit=15)
+    minimizer = ift.RelaxedNewton(IC2)
+
+
+    H = Hamiltonian(likelihood, inverter)
+    H, convergence = minimizer(H)
+    result = s.at(H.position).value
+
+

demos/nonlinear_critical_filter.py

-import nifty5 as ift
-from nifty5 import Exponential
-import numpy as np
-np.random.seed(42)
-
-
-def adjust_zero_mode(m0, t0):
-    mtmp = m0.to_global_data().copy()
-    zero_position = len(m0.shape)*(0,)
-    zero_mode = mtmp[zero_position]
-    mtmp[zero_position] = zero_mode / abs(zero_mode)
-    ttmp = t0.to_global_data().copy()
-    ttmp[0] += 2 * np.log(abs(zero_mode))
-    return (ift.Field.from_global_data(m0.domain, mtmp),
-            ift.Field.from_global_data(t0.domain, ttmp))
-
-
-if __name__ == "__main__":
-    noise_level = 1.
-    p_spec = (lambda k: (1. / (k + 1) ** 2))
-
-    # nonlinearity = Linear()
-    nonlinearity = Exponential()
-    # Set up position space
-    # s_space = ift.RGSpace([1024])
-    s_space = ift.HPSpace(32)
-    h_space = s_space.get_default_codomain()
-
-    # Define harmonic transformation and associated harmonic space
-    HT = ift.HarmonicTransformOperator(h_space, target=s_space)
-
-    # Setting up power space
-    p_space = ift.PowerSpace(h_space,
-                             binbounds=ift.PowerSpace.useful_binbounds(
-                                 h_space, logarithmic=True))
-    # Choosing the prior correlation structure and defining
-    # correlation operator
-    p = ift.PS_field(p_space, p_spec)
-    log_p = ift.log(p)
-    S = ift.create_power_operator(h_space, lambda k: 1.)
-
-    # Drawing a sample sh from the prior distribution in harmonic space
-    sh = S.draw_sample()
-
-    # Choosing the measurement instrument
-    # Instrument = SmoothingOperator(s_space, sigma=0.01)
-    mask = np.ones(s_space.shape)
-    mask[6000:8000] = 0.
-    mask = ift.Field.from_global_data(s_space, mask)
-
-    MaskOperator = ift.DiagonalOperator(mask)
-    R = ift.GeometryRemover(s_space)
-    R = R*MaskOperator
-    # R = R*HT
-    # R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
-    #                                                response_sigma)
-    MeasurementOperator = R
-
-    d_space = MeasurementOperator.target
-
-    Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
-    power = Distributor(ift.exp(0.5*log_p))
-    # Creating the mock data
-    true_sky = nonlinearity(HT(power*sh))
-    noiseless_data = MeasurementOperator(true_sky)
-    noise_amplitude = noiseless_data.val.std()*noise_level
-    N = ift.ScalingOperator(noise_amplitude**2, d_space)
-    n = N.draw_sample()
-    # Creating the mock data
-    d = noiseless_data + n
-
-    m0 = ift.full(h_space, 1e-7)
-    t0 = ift.full(p_space, -4.)
-    power0 = Distributor.times(ift.exp(0.5 * t0))
-
-    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
-                                     tol_abs_gradnorm=1e-3)
-    LS = ift.LineSearchStrongWolfe(c2=0.02)
-    minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
-
-    IC = ift.GradientNormController(iteration_limit=500,
-                                    tol_abs_gradnorm=1e-3)
-
-    for i in range(20):
-        power0 = Distributor(ift.exp(0.5*t0))
-        map0_energy = ift.library.NonlinearWienerFilterEnergy(
-            m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S, IC)
-
-        # Minimization with chosen minimizer
-        map0_energy, convergence = minimizer(map0_energy)
-        m0 = map0_energy.position
-
-        # Updating parameters for correlation structure reconstruction
-        D0 = map0_energy.curvature
-
-        # Initializing the power energy with updated parameters
-        power0_energy = ift.library.NonlinearPowerEnergy(
-            position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
-            Instrument=MeasurementOperator, nonlinearity=nonlinearity,
-            Distributor=Distributor, sigma=1., samples=2, iteration_controller=IC)
-
-        power0_energy = minimizer(power0_energy)[0]
-
-        # Setting new power spectrum
-        t0 = power0_energy.position
-
-        # break degeneracy between amplitude and excitation by setting
-        # excitation monopole to 1
-        m0, t0 = adjust_zero_mode(m0, t0)
-
-    plotdict = {"colormap": "Planck-like"}
-    ift.plot(true_sky, name="true_sky.png", **plotdict)
-    ift.plot(nonlinearity(HT(power0*m0)),
-             name="reconstructed_sky.png", **plotdict)
-    ift.plot(MeasurementOperator.adjoint_times(d), name="data.png", **plotdict)
-    ift.plot([ift.exp(t0), p], name="ps.png")

demos/nonlinear_wiener_filter.py

-import nifty5 as ift
-from nifty5.library.nonlinearities import Linear, Exponential, Tanh
-import numpy as np
-np.random.seed(20)
-
-if __name__ == "__main__":
-
-    noise_level = 0.3
-    correlation_length = 0.1
-    p_spec = lambda k: (1. / (k*correlation_length + 1) ** 4)
-
-    nonlinearity = Tanh()
-    # nonlinearity = Linear()
-    # nonlinearity = Exponential()
-
-    # Set up position space
-    s_space = ift.RGSpace(1024)
-    h_space = s_space.get_default_codomain()
-
-    # Define harmonic transformation and associated harmonic space
-    HT = ift.HarmonicTransformOperator(h_space, target=s_space)
-
-    S = ift.ScalingOperator(1., h_space)
-
-    # Drawing a sample sh from the prior distribution in harmonic space
-    sh = S.draw_sample()
-
-    # Choosing the measurement instrument
-    # Instrument = SmoothingOperator(s_space, sigma=0.01)
-    mask = np.ones(s_space.shape)
-    mask[600:800] = 0.
-    mask = ift.Field.from_global_data(s_space, mask)
-
-    R = ift.GeometryRemover(s_space) * ift.DiagonalOperator(mask)
-
-    d_space = R.target
-
-    p_op = ift.create_power_operator(h_space, p_spec)
-    power = ift.sqrt(p_op(ift.full(h_space, 1.)))
-
-    # Creating the mock data
-    true_sky = nonlinearity(HT(power*sh))
-    noiseless_data = R(true_sky)
-    noise_amplitude = noiseless_data.val.std()*noise_level
-    N = ift.ScalingOperator(noise_amplitude**2, d_space)
-    n = N.draw_sample()
-    # Creating the mock data
-    d = noiseless_data + n
-
-    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
-                                     tol_abs_gradnorm=1e-4)
-    LS = ift.LineSearchStrongWolfe(c2=0.02)
-    minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
-
-    IC = ift.GradientNormController(iteration_limit=2000,
-                                    tol_abs_gradnorm=1e-3)
-
-    # initial guess
-    m = ift.full(h_space, 1e-7)
-    map_energy = ift.library.NonlinearWienerFilterEnergy(
-        m, d, R, nonlinearity, HT, power, N, S, IC)
-
-    # Minimization with chosen minimizer
-    map_energy, convergence = minimizer(map_energy)
-    m = map_energy.position
-
-    recsky = nonlinearity(HT(power*m))
-    data = R.adjoint_times(d)
-    lo = np.min([true_sky.min(), recsky.min(), data.min()])
-    hi = np.max([true_sky.max(), recsky.max(), data.max()])
-    plotdict = {"colormap": "Planck-like", "ymin": lo, "ymax": hi}
-    ift.plot(true_sky, name="true_sky.png", **plotdict)
-    ift.plot(recsky, name="reconstructed_sky.png", **plotdict)
-    ift.plot(data, name="data.png", **plotdict)

demos/paper_demos/cartesian_wiener_filter.py

-import numpy as np
-import nifty5 as ift
-
-if __name__ == "__main__":
-    signal_to_noise = 0.5  # The signal to noise ratio
-
-    # Setting up parameters
-    L_1 = 2.                   # Total side-length of the domain
-    N_pixels_1 = 512           # Grid resolution (pixels per axis)
-    L_2 = 2.                   # Total side-length of the domain
-    N_pixels_2 = 512           # Grid resolution (pixels per axis)
-    correlation_length_1 = 1.
-    field_variance_1 = 2.      # Variance of field in position space
-    response_sigma_1 = 0.05    # Smoothing length of response
-    correlation_length_2 = 1.
-    field_variance_2 = 2.      # Variance of field in position space
-    response_sigma_2 = 0.01    # Smoothing length of response
-
-    def power_spectrum_1(k):   # note: field_variance**2 = a*k_0/4.
-        a = 4 * correlation_length_1 * field_variance_1**2
-        return a / (1 + k * correlation_length_1) ** 4.
-
-    def power_spectrum_2(k):  # note: field_variance**2 = a*k_0/4.
-        a = 4 * correlation_length_2 * field_variance_2**2
-        return a / (1 + k * correlation_length_2) ** 2.5
-
-    signal_space_1 = ift.RGSpace([N_pixels_1], distances=L_1/N_pixels_1)
-    harmonic_space_1 = signal_space_1.get_default_codomain()
-    signal_space_2 = ift.RGSpace([N_pixels_2], distances=L_2/N_pixels_2)
-    harmonic_space_2 = signal_space_2.get_default_codomain()
-
-    signal_domain = ift.DomainTuple.make((signal_space_1, signal_space_2))
-    harmonic_domain = ift.DomainTuple.make((harmonic_space_1,
-                                            harmonic_space_2))
-
-    ht_1 = ift.HarmonicTransformOperator(harmonic_domain, space=0)
-    ht_2 = ift.HarmonicTransformOperator(ht_1.target, space=1)
-    ht = ht_2*ht_1
-
-    S = (ift.create_power_operator(harmonic_domain, power_spectrum_1, 0) *
-         ift.create_power_operator(harmonic_domain, power_spectrum_2, 1))
-
-    np.random.seed(10)
-    mock_signal = S.draw_sample()
-
-    # Setting up a exemplary response
-    N1_10 = int(N_pixels_1/10)
-    mask_1 = np.ones(signal_space_1.shape)
-    mask_1[N1_10*7:N1_10*9] = 0.
-    mask_1 = ift.Field.from_global_data(signal_space_1, mask_1)
-
-    N2_10 = int(N_pixels_2/10)
-    mask_2 = np.ones(signal_space_2.shape)
-    mask_2[N2_10*7:N2_10*9] = 0.
-    mask_2 = ift.Field.from_global_data(signal_space_2, mask_2)
-
-    R = ift.GeometryRemover(signal_domain)
-    R = R*ift.DiagonalOperator(mask_1, signal_domain, spaces=0)
-    R = R*ift.DiagonalOperator(mask_2, signal_domain, spaces=1)
-    R = R*ht
-    R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 0,
-                                                   response_sigma_1)
-    R = R * ift.create_harmonic_smoothing_operator(harmonic_domain, 1,
-                                                   response_sigma_2)
-    data_domain = R.target
-
-    noiseless_data = R(mock_signal)
-    noise_amplitude = noiseless_data.val.std()/signal_to_noise
-    # Setting up the noise covariance and drawing a random noise realization
-    N = ift.ScalingOperator(noise_amplitude**2, data_domain)
-    noise = N.draw_sample()
-    data = noiseless_data + noise
-
-    # Wiener filter
-    j = R.adjoint_times(N.inverse_times(data))
-    ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
-    sampling_ctrl = ift.GradientNormController(name="sampling",
-                                               tol_abs_gradnorm=1e2)
-    wiener_curvature = ift.library.WienerFilterCurvature(
-        S=S, N=N, R=R, iteration_controller=ctrl,
-        iteration_controller_sampling=sampling_ctrl)
-
-    m_k = wiener_curvature.inverse_times(j)
-    m = ht(m_k)
-
-    plotdict = {"colormap": "Planck-like"}
-    plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))