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

Wrap all scripts in functions

This is best practise in python scripts in order no to clutter the global
namespace.
parent 1200b36a
Pipeline #75983 passed with stages
in 13 minutes and 50 seconds
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2019 Max-Planck-Society # Copyright(C) 2019-2020 Max-Planck-Society
# Author: Martin Reinecke # Author: Martin Reinecke
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -23,68 +23,73 @@ import numpy as np ...@@ -23,68 +23,73 @@ import numpy as np
import nifty7 as ift import nifty7 as ift
N0s, a0s, b0s, c0s = [], [], [], []
def main():
for ii in range(10, 26): N0s, a0s, b0s, c0s = [], [], [], []
nu = 1024
nv = 1024 for ii in range(10, 26):
N = int(2**ii) nu = 1024
print('N = {}'.format(N)) nv = 1024
N = int(2**ii)
rng = ift.random.current_rng() print('N = {}'.format(N))
uv = rng.uniform(-.5, .5, (N,2))
vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N) rng = ift.random.current_rng()
uv = rng.uniform(-.5, .5, (N, 2))
uvspace = ift.RGSpace((nu, nv)) vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
visspace = ift.UnstructuredDomain(N) uvspace = ift.RGSpace((nu, nv))
img = rng.standard_normal((nu, nv)) visspace = ift.UnstructuredDomain(N)
img = ift.makeField(uvspace, img)
img = rng.standard_normal((nu, nv))
t0 = time() img = ift.makeField(uvspace, img)
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
vis = ift.makeField(visspace, vis) t0 = time()
op = GM.getFull().adjoint GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
t1 = time() vis = ift.makeField(visspace, vis)
op(img).val op = GM.getFull().adjoint
t2 = time() t1 = time()
op.adjoint(vis).val op(img).val
t3 = time() t2 = time()
print(t2-t1, t3-t2) op.adjoint(vis).val
N0s.append(N) t3 = time()
a0s.append(t1 - t0) print(t2-t1, t3-t2)
b0s.append(t2 - t1) N0s.append(N)
c0s.append(t3 - t2) a0s.append(t1 - t0)
b0s.append(t2 - t1)
print('Measure rest operator') c0s.append(t3 - t2)
sc = ift.StatCalculator()
op = GM.getRest().adjoint print('Measure rest operator')
for _ in range(10): sc = ift.StatCalculator()
t0 = time() op = GM.getRest().adjoint
res = op(img) for _ in range(10):
sc.add(time() - t0) t0 = time()
t_fft = sc.mean res = op(img)
print('FFT shape', res.shape) sc.add(time() - t0)
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.legend() plt.scatter(N0s, a0s, label='Gridder mr')
# no idea why this is necessary, but if it is omitted, the range is wrong plt.legend()
plt.ylim(min(a0s), max(a0s)) # no idea why this is necessary, but if it is omitted, the range is wrong
plt.ylabel('time [s]') plt.ylim(min(a0s), max(a0s))
plt.title('Initialization') plt.ylabel('time [s]')
plt.loglog() plt.title('Initialization')
plt.savefig('bench0.png') plt.loglog()
plt.close() plt.savefig('bench0.png')
plt.close()
plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times') plt.scatter(N0s, b0s, color='k', marker='^', label='Gridder mr times')
plt.axhline(sc.mean, label='FFT') plt.scatter(N0s, c0s, color='k', label='Gridder mr adjoint times')
plt.axhline(sc.mean + np.sqrt(sc.var)) plt.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean - np.sqrt(sc.var)) plt.axhline(sc.mean + np.sqrt(sc.var))
plt.legend() plt.axhline(sc.mean - np.sqrt(sc.var))
plt.ylabel('time [s]') plt.legend()
plt.title('Apply') plt.ylabel('time [s]')
plt.loglog() plt.title('Apply')
plt.savefig('bench1.png') plt.loglog()
plt.close() plt.savefig('bench1.png')
plt.close()
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -26,7 +26,8 @@ import numpy as np ...@@ -26,7 +26,8 @@ import numpy as np
import nifty7 as ift import nifty7 as ift
if __name__ == '__main__':
def main():
# 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:
...@@ -84,3 +85,7 @@ if __name__ == '__main__': ...@@ -84,3 +85,7 @@ if __name__ == '__main__':
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", name="bernoulli.png") plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -38,14 +38,14 @@ def make_checkerboard_mask(position_space): ...@@ -38,14 +38,14 @@ def make_checkerboard_mask(position_space):
return mask return mask
def make_random_mask(): def make_random_mask(domain):
# Random mask for spherical mode # Random mask for spherical mode
mask = ift.from_random(position_space, 'pm1') mask = ift.from_random(domain, 'pm1')
mask = (mask + 1)/2 mask = (mask + 1)/2
return mask.val return mask.val
if __name__ == '__main__': def main():
# Choose space on which the signal field is defined # Choose space on which the signal field is defined
if len(sys.argv) == 2: if len(sys.argv) == 2:
mode = int(sys.argv[1]) mode = int(sys.argv[1])
...@@ -63,7 +63,7 @@ if __name__ == '__main__': ...@@ -63,7 +63,7 @@ if __name__ == '__main__':
else: else:
# Sphere with half of its pixels randomly masked # Sphere with half of its pixels randomly masked
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
mask = make_random_mask() mask = make_random_mask(position_space)
# Specify harmonic space corresponding to signal # Specify harmonic space corresponding to signal
harmonic_space = position_space.get_default_codomain() harmonic_space = position_space.get_default_codomain()
...@@ -148,3 +148,7 @@ if __name__ == '__main__': ...@@ -148,3 +148,7 @@ if __name__ == '__main__':
title='Residuals') title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename) plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename)
print("Saved results as '{}'.".format(filename)) print("Saved results as '{}'.".format(filename))
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -28,22 +28,20 @@ import numpy as np ...@@ -28,22 +28,20 @@ import numpy as np
import nifty7 as ift import nifty7 as ift
def exposure_2d(): def exposure_2d(domain):
# Structured exposure for 2D mode # Structured exposure for 2D mode
x_shape, y_shape = position_space.shape x_shape, y_shape = domain.shape
exposure = np.ones(domain.shape)
exposure = np.ones(position_space.shape)
exposure[x_shape//3:x_shape//2, :] *= 2. exposure[x_shape//3:x_shape//2, :] *= 2.
exposure[x_shape*4//5:x_shape, :] *= .1 exposure[x_shape*4//5:x_shape, :] *= .1
exposure[x_shape//2:x_shape*3//2, :] *= 3. exposure[x_shape//2:x_shape*3//2, :] *= 3.
exposure[:, x_shape//3:x_shape//2] *= 2. exposure[:, x_shape//3:x_shape//2] *= 2.
exposure[:, x_shape*4//5:x_shape] *= .1 exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3. exposure[:, x_shape//2:x_shape*3//2] *= 3.
return ift.Field.from_raw(domain, exposure)
return ift.Field.from_raw(position_space, exposure)
if __name__ == '__main__': def main():
# Choose space on which the signal field is defined # Choose space on which the signal field is defined
if len(sys.argv) == 2: if len(sys.argv) == 2:
mode = int(sys.argv[1]) mode = int(sys.argv[1])
...@@ -57,7 +55,7 @@ if __name__ == '__main__': ...@@ -57,7 +55,7 @@ if __name__ == '__main__':
elif mode == 1: elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure # Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512]) position_space = ift.RGSpace([512, 512])
exposure = exposure_2d() exposure = exposure_2d(position_space)
else: else:
# Sphere with uniform exposure of 100 # Sphere with uniform exposure of 100
position_space = ift.HPSpace(128) position_space = ift.HPSpace(128)
...@@ -118,3 +116,7 @@ if __name__ == '__main__': ...@@ -118,3 +116,7 @@ if __name__ == '__main__':
plot.add(reconst - signal, title='Residuals') plot.add(reconst - signal, title='Residuals')
plot.output(xsize=12, ysize=10, name=filename) plot.output(xsize=12, ysize=10, name=filename)
print("Saved results as '{}'.".format(filename)) print("Saved results as '{}'.".format(filename))
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -44,7 +44,7 @@ def radial_los(n_los): ...@@ -44,7 +44,7 @@ def radial_los(n_los):
return starts, ends return starts, ends
if __name__ == '__main__': def main():
# Choose between random line-of-sight response (mode=0) and radial lines # Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1) # of sight (mode=1)
if len(sys.argv) == 2: if len(sys.argv) == 2:
...@@ -60,10 +60,10 @@ if __name__ == '__main__': ...@@ -60,10 +60,10 @@ if __name__ == '__main__':
# see 'getting_started_4_CorrelatedFields.ipynb'. # see 'getting_started_4_CorrelatedFields.ipynb'.
cfmaker = ift.CorrelatedFieldMaker.make( cfmaker = ift.CorrelatedFieldMaker.make(
offset_mean = 0.0, # 0. offset_mean= 0.0,
offset_std_mean = 1e-3, # 1e-3 offset_std_mean= 1e-3,
offset_std_std = 1e-6, # 1e-6 offset_std_std= 1e-6,
prefix = '') prefix='')
fluctuations_dict = { fluctuations_dict = {
# Amplitude of field fluctuations # Amplitude of field fluctuations
...@@ -72,7 +72,7 @@ if __name__ == '__main__': ...@@ -72,7 +72,7 @@ if __name__ == '__main__':
# Exponent of power law power spectrum component # Exponent of power law power spectrum component
'loglogavgslope_mean': -2.0, # -3.0 'loglogavgslope_mean': -2.0, # -3.0
'loglogavgslope_stddev': 0.5, # 0.5 'loglogavgslope_stddev': 0.5, # 0.5
# Amplitude of integrated Wiener process power spectrum component # Amplitude of integrated Wiener process power spectrum component
'flexibility_mean': 2.5, # 1.0 'flexibility_mean': 2.5, # 1.0
...@@ -163,3 +163,7 @@ if __name__ == '__main__': ...@@ -163,3 +163,7 @@ if __name__ == '__main__':
linewidth=[1.]*len(powers) + [3., 3.]) linewidth=[1.]*len(powers) + [3., 3.])
plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res) plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
print("Saved results as '{}'.".format(filename_res)) print("Saved results as '{}'.".format(filename_res))
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -55,7 +55,7 @@ def radial_los(n_los): ...@@ -55,7 +55,7 @@ def radial_los(n_los):
return starts, ends return starts, ends
if __name__ == '__main__': def main():
# Choose between random line-of-sight response (mode=0) and radial lines # Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1) # of sight (mode=1)
if len(sys.argv) == 2: if len(sys.argv) == 2:
...@@ -83,7 +83,7 @@ if __name__ == '__main__': ...@@ -83,7 +83,7 @@ if __name__ == '__main__':
A2 = cfmaker.normalized_amplitudes[1] A2 = cfmaker.normalized_amplitudes[1]
DC = SingleDomain(correlated_field.target, position_space) DC = SingleDomain(correlated_field.target, position_space)
## Apply a nonlinearity # Apply a nonlinearity
signal = DC @ ift.sigmoid(correlated_field) signal = DC @ ift.sigmoid(correlated_field)
# Build the line-of-sight response and define signal response # Build the line-of-sight response and define signal response
...@@ -118,7 +118,7 @@ if __name__ == '__main__': ...@@ -118,7 +118,7 @@ if __name__ == '__main__':
ic_newton.enable_logging() ic_newton.enable_logging()
minimizer = ift.NewtonCG(ic_newton, enable_logging=True) minimizer = ift.NewtonCG(ic_newton, enable_logging=True)
## number of samples used to estimate the KL # number of samples used to estimate the KL
N_samples = 20 N_samples = 20
# Set up likelihood and information Hamiltonian # Set up likelihood and information Hamiltonian
...@@ -157,7 +157,6 @@ if __name__ == '__main__': ...@@ -157,7 +157,6 @@ if __name__ == '__main__':
name=filename.format("loop_{:02d}".format(i))) name=filename.format("loop_{:02d}".format(i)))
# Done, draw posterior samples # Done, draw posterior samples
Nsamples = 20
KL = ift.MetricGaussianKL(mean, H, N_samples) KL = ift.MetricGaussianKL(mean, H, N_samples)
sc = ift.StatCalculator() sc = ift.StatCalculator()
scA1 = ift.StatCalculator() scA1 = ift.StatCalculator()
...@@ -189,3 +188,7 @@ if __name__ == '__main__': ...@@ -189,3 +188,7 @@ if __name__ == '__main__':
linewidth=[1.]*len(powers2) + [3., 3.]) linewidth=[1.]*len(powers2) + [3., 3.])
plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res) plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res)
print("Saved results as '{}'.".format(filename_res)) print("Saved results as '{}'.".format(filename_res))
if __name__ == '__main__':
main()
...@@ -34,7 +34,7 @@ from matplotlib.colors import LogNorm ...@@ -34,7 +34,7 @@ from matplotlib.colors import LogNorm
import nifty7 as ift import nifty7 as ift
if __name__ == '__main__': def main():
dom = ift.UnstructuredDomain(1) dom = ift.UnstructuredDomain(1)
scale = 10 scale = 10
...@@ -119,3 +119,7 @@ if __name__ == '__main__': ...@@ -119,3 +119,7 @@ if __name__ == '__main__':
ift.logger.info('Finished') ift.logger.info('Finished')
# Uncomment the following line in order to leave the plots open # Uncomment the following line in order to leave the plots open
# plt.show() # plt.show()
if __name__ == '__main__':
main()
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2019 Max-Planck-Society # Copyright(C) 2019-2020 Max-Planck-Society
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
...@@ -37,55 +37,52 @@ def convtest(test_signal, delta, func): ...@@ -37,55 +37,52 @@ def convtest(test_signal, delta, func):
# Plot results # Plot results
plot = ift.Plot() plot = ift.Plot()
plot.add(signal, title='Signal') plot.add(test_signal, title='Signal')
plot.add(conv_signal, title='Signal Convolved') plot.add(conv_signal, title='Signal Convolved')
plot.add(cac_signal, title='Signal, Conv, Adj-Conv') plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
plot.add(conv_delta, title='Kernel') plot.add(conv_delta, title='Kernel')
plot.output() plot.output()
# Healpix test def main():
nside = 64 # Healpix test
npix = 12 * nside * nside nside = 64
npix = 12 * nside * nside
domain = ift.HPSpace(nside) domain = ift.HPSpace(nside)
# Define test signal (some point sources) # Define test signal (some point sources)
signal_vals = np.zeros(npix, dtype=np.float64) signal_vals = np.zeros(npix, dtype=np.float64)
for i in range(0, npix, npix//12 + 27): for i in range(0, npix, npix//12 + 27):
signal_vals[i] = 500. signal_vals[i] = 500.
signal = ift.makeField(domain, signal_vals) signal = ift.makeField(domain, signal_vals)
delta_vals = np.zeros(npix, dtype=np.float64) delta_vals = np.zeros(npix, dtype=np.float64)
delta_vals[0] = 1.0 delta_vals[0] = 1.0
delta = ift.makeField(domain, delta_vals) delta = ift.makeField(domain, delta_vals)
# Define kernel function
def func(theta):
ct = np.cos(theta)
return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
# Define kernel function convtest(signal, delta, func)
def func(theta):
ct = np.cos(theta)
return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
domain = ift.RGSpace((100, 100))
# Define test signal (some point sources)
signal_vals = np.zeros(domain.shape, dtype=np.float64)
signal_vals[35, 70] = 5000.
signal_vals[45, 8] = 5000.
signal = ift.makeField(domain, signal_vals)
convtest(signal, delta, func) # Define delta signal, generate kernel image
delta_vals = np.zeros(domain.shape, dtype=np.float64)
delta_vals[0, 0] = 1.0
delta = ift.makeField(domain, delta_vals)
domain = ift.RGSpace((100, 100)) convtest(signal, delta, lambda d: 1. * np.logical_and(d > 0.1, d <= 0.2))
# Define test signal (some point sources)
signal_vals = np.zeros(domain.shape, dtype=np.float64)
signal_vals[35, 70] = 5000.
signal_vals[45, 8] = 5000.
signal = ift.makeField(domain, signal_vals)
# Define delta signal, generate kernel image
delta_vals = np.zeros(domain.shape, dtype=np.float64)
delta_vals[0, 0] = 1.0
delta = ift.makeField(domain, delta_vals)
if __name__ == '__main__':
# Define kernel function main()
def func(dist):
return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
convtest(signal, delta, func)
...@@ -11,7 +11,8 @@ ...@@ -11,7 +11,8 @@
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2020 Max-Planck-Society