Commit ab18a1b7 authored by Philipp Arras's avatar Philipp Arras
Browse files

Improvements Bernoulli demo

parent 59a4902c
...@@ -15,34 +15,33 @@ ...@@ -15,34 +15,33 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
##################################################################### #####################################################################
# Bernoulli reconstruction # Bernoulli reconstruction
# reconstruct an event probability field with values between 0 and 1 # Reconstruct an event probability field with values between 0 and 1
# from the observed events # from the observed events
# 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2) # 1D (set mode=0), 2D (mode=1), or on the sphere (mode=2)
##################################################################### #####################################################################
import nifty5 as ift
import numpy as np import numpy as np
import nifty5 as ift
if __name__ == '__main__': if __name__ == '__main__':
# FIXME ABOUT THIS CODE
np.random.seed(41) np.random.seed(41)
# Set up the position space of the signal # Set up the position space of the signal
mode = 2 mode = 2
if mode == 0: if mode == 0:
# One dimensional regular grid # One-dimensional regular grid
position_space = ift.RGSpace(1024) position_space = ift.RGSpace(1024)
elif mode == 1: elif mode == 1:
# Two dimensional regular grid # Two-dimensional regular grid
position_space = ift.RGSpace([512, 512]) position_space = ift.RGSpace([512, 512])
else: else:
# Sphere # Sphere
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
# Defining harmonic space and transform # Define harmonic space and transform
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space) HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
...@@ -50,15 +49,13 @@ if __name__ == '__main__': ...@@ -50,15 +49,13 @@ if __name__ == '__main__':
# Define power spectrum and amplitudes # Define power spectrum and amplitudes
def sqrtpspec(k): def sqrtpspec(k):
return 1. / (20. + k**2) return 1./(20. + k**2)
A = ift.create_power_operator(harmonic_space, sqrtpspec) A = ift.create_power_operator(harmonic_space, sqrtpspec)
# Set up a sky model # Set up a sky model and instrumental response
sky = ift.positive_tanh(HT(A)) sky = ift.positive_tanh(HT(A))
GR = ift.GeometryRemover(position_space) GR = ift.GeometryRemover(position_space)
# Set up instrumental response
R = GR R = GR
# Generate mock data # Generate mock data
...@@ -71,8 +68,8 @@ if __name__ == '__main__': ...@@ -71,8 +68,8 @@ if __name__ == '__main__':
# Compute likelihood and Hamiltonian # Compute likelihood and Hamiltonian
position = ift.from_random('normal', harmonic_space) position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(data)(p) likelihood = ift.BernoulliEnergy(data)(p)
ic_newton = ift.DeltaEnergyController(name='Newton', iteration_limit=100, ic_newton = ift.DeltaEnergyController(
tol_rel_deltaE=1e-8) name='Newton', iteration_limit=100, tol_rel_deltaE=1e-8)
minimizer = ift.NewtonCG(ic_newton) minimizer = ift.NewtonCG(ic_newton)
ic_sampling = ift.GradientNormController(iteration_limit=100) ic_sampling = ift.GradientNormController(iteration_limit=100)
...@@ -88,5 +85,4 @@ if __name__ == '__main__': ...@@ -88,5 +85,4 @@ if __name__ == '__main__':
plot.add(reconstruction, title='reconstruction') plot.add(reconstruction, title='reconstruction')
plot.add(GR.adjoint_times(data), title='data') plot.add(GR.adjoint_times(data), title='data')
plot.add(sky(mock_position), title='truth') plot.add(sky(mock_position), title='truth')
plot.output(nx=3, xsize=16, ysize=9, title="results", plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
name="bernoulli.png")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment