Commit 682a34b7 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'python_bestpractise' into 'NIFTy_7'

Wrap all scripts in functions

See merge request !524
parents be9f5044 de77fae1
Pipeline #75988 passed with stages
in 13 minutes and 38 seconds
......@@ -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) 2019 Max-Planck-Society
# Copyright(C) 2019-2020 Max-Planck-Society
# Author: Martin Reinecke
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -23,68 +23,73 @@ import numpy as np
import nifty7 as ift
N0s, a0s, b0s, c0s = [], [], [], []
for ii in range(10, 26):
nu = 1024
nv = 1024
N = int(2**ii)
print('N = {}'.format(N))
rng = ift.random.current_rng()
uv = rng.uniform(-.5, .5, (N,2))
vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
uvspace = ift.RGSpace((nu, nv))
visspace = ift.UnstructuredDomain(N)
img = rng.standard_normal((nu, nv))
img = ift.makeField(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
vis = ift.makeField(visspace, vis)
op = GM.getFull().adjoint
t1 = time()
op(img).val
t2 = time()
op.adjoint(vis).val
t3 = time()
print(t2-t1, t3-t2)
N0s.append(N)
a0s.append(t1 - t0)
b0s.append(t2 - t1)
c0s.append(t3 - t2)
print('Measure rest operator')
sc = ift.StatCalculator()
op = GM.getRest().adjoint
for _ in range(10):
t0 = time()
res = op(img)
sc.add(time() - t0)
t_fft = sc.mean
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.legend()
# no idea why this is necessary, but if it is omitted, the range is wrong
plt.ylim(min(a0s), max(a0s))
plt.ylabel('time [s]')
plt.title('Initialization')
plt.loglog()
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.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean + np.sqrt(sc.var))
plt.axhline(sc.mean - np.sqrt(sc.var))
plt.legend()
plt.ylabel('time [s]')
plt.title('Apply')
plt.loglog()
plt.savefig('bench1.png')
plt.close()
def main():
N0s, a0s, b0s, c0s = [], [], [], []
for ii in range(10, 26):
nu = 1024
nv = 1024
N = int(2**ii)
print('N = {}'.format(N))
rng = ift.random.current_rng()
uv = rng.uniform(-.5, .5, (N, 2))
vis = rng.normal(0., 1., N) + 1j*rng.normal(0., 1., N)
uvspace = ift.RGSpace((nu, nv))
visspace = ift.UnstructuredDomain(N)
img = rng.standard_normal((nu, nv))
img = ift.makeField(uvspace, img)
t0 = time()
GM = ift.GridderMaker(uvspace, eps=1e-7, uv=uv)
vis = ift.makeField(visspace, vis)
op = GM.getFull().adjoint
t1 = time()
op(img).val
t2 = time()
op.adjoint(vis).val
t3 = time()
print(t2-t1, t3-t2)
N0s.append(N)
a0s.append(t1 - t0)
b0s.append(t2 - t1)
c0s.append(t3 - t2)
print('Measure rest operator')
sc = ift.StatCalculator()
op = GM.getRest().adjoint
for _ in range(10):
t0 = time()
res = op(img)
sc.add(time() - t0)
print('FFT shape', res.shape)
plt.scatter(N0s, a0s, label='Gridder mr')
plt.legend()
# no idea why this is necessary, but if it is omitted, the range is wrong
plt.ylim(min(a0s), max(a0s))
plt.ylabel('time [s]')
plt.title('Initialization')
plt.loglog()
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.axhline(sc.mean, label='FFT')
plt.axhline(sc.mean + np.sqrt(sc.var))
plt.axhline(sc.mean - np.sqrt(sc.var))
plt.legend()
plt.ylabel('time [s]')
plt.title('Apply')
plt.loglog()
plt.savefig('bench1.png')
plt.close()
if __name__ == '__main__':
main()
......@@ -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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -26,7 +26,8 @@ import numpy as np
import nifty7 as ift
if __name__ == '__main__':
def main():
# Set up the position space of the signal
mode = 2
if mode == 0:
......@@ -84,3 +85,7 @@ if __name__ == '__main__':
plot.add(GR.adjoint_times(data), title='data')
plot.add(sky(mock_position), title='truth')
plot.output(nx=3, xsize=16, ysize=9, title="results", name="bernoulli.png")
if __name__ == '__main__':
main()
......@@ -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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -38,14 +38,14 @@ def make_checkerboard_mask(position_space):
return mask
def make_random_mask():
def make_random_mask(domain):
# Random mask for spherical mode
mask = ift.from_random(position_space, 'pm1')
mask = ift.from_random(domain, 'pm1')
mask = (mask + 1)/2
return mask.val
if __name__ == '__main__':
def main():
# Choose space on which the signal field is defined
if len(sys.argv) == 2:
mode = int(sys.argv[1])
......@@ -63,7 +63,7 @@ if __name__ == '__main__':
else:
# Sphere with half of its pixels randomly masked
position_space = ift.HPSpace(128)
mask = make_random_mask()
mask = make_random_mask(position_space)
# Specify harmonic space corresponding to signal
harmonic_space = position_space.get_default_codomain()
......@@ -148,3 +148,7 @@ if __name__ == '__main__':
title='Residuals')
plot.output(nx=2, ny=2, xsize=10, ysize=10, name=filename)
print("Saved results as '{}'.".format(filename))
if __name__ == '__main__':
main()
......@@ -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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -28,22 +28,20 @@ import numpy as np
import nifty7 as ift
def exposure_2d():
def exposure_2d(domain):
# Structured exposure for 2D mode
x_shape, y_shape = position_space.shape
exposure = np.ones(position_space.shape)
x_shape, y_shape = domain.shape
exposure = np.ones(domain.shape)
exposure[x_shape//3:x_shape//2, :] *= 2.
exposure[x_shape*4//5:x_shape, :] *= .1
exposure[x_shape//2:x_shape*3//2, :] *= 3.
exposure[:, x_shape//3:x_shape//2] *= 2.
exposure[:, x_shape*4//5:x_shape] *= .1
exposure[:, x_shape//2:x_shape*3//2] *= 3.
return ift.Field.from_raw(position_space, exposure)
return ift.Field.from_raw(domain, exposure)
if __name__ == '__main__':
def main():
# Choose space on which the signal field is defined
if len(sys.argv) == 2:
mode = int(sys.argv[1])
......@@ -57,7 +55,7 @@ if __name__ == '__main__':
elif mode == 1:
# Two-dimensional regular grid with inhomogeneous exposure
position_space = ift.RGSpace([512, 512])
exposure = exposure_2d()
exposure = exposure_2d(position_space)
else:
# Sphere with uniform exposure of 100
position_space = ift.HPSpace(128)
......@@ -118,3 +116,7 @@ if __name__ == '__main__':
plot.add(reconst - signal, title='Residuals')
plot.output(xsize=12, ysize=10, name=filename)
print("Saved results as '{}'.".format(filename))
if __name__ == '__main__':
main()
......@@ -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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -44,7 +44,7 @@ def radial_los(n_los):
return starts, ends
if __name__ == '__main__':
def main():
# Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1)
if len(sys.argv) == 2:
......@@ -60,10 +60,10 @@ if __name__ == '__main__':
# see 'getting_started_4_CorrelatedFields.ipynb'.
cfmaker = ift.CorrelatedFieldMaker.make(
offset_mean = 0.0, # 0.
offset_std_mean = 1e-3, # 1e-3
offset_std_std = 1e-6, # 1e-6
prefix = '')
offset_mean= 0.0,
offset_std_mean= 1e-3,
offset_std_std= 1e-6,
prefix='')
fluctuations_dict = {
# Amplitude of field fluctuations
......@@ -72,7 +72,7 @@ if __name__ == '__main__':
# Exponent of power law power spectrum component
'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
'flexibility_mean': 2.5, # 1.0
......@@ -163,3 +163,7 @@ if __name__ == '__main__':
linewidth=[1.]*len(powers) + [3., 3.])
plot.output(ny=1, nx=3, xsize=24, ysize=6, name=filename_res)
print("Saved results as '{}'.".format(filename_res))
if __name__ == '__main__':
main()
......@@ -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-2019 Max-Planck-Society
# Copyright(C) 2013-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -55,7 +55,7 @@ def radial_los(n_los):
return starts, ends
if __name__ == '__main__':
def main():
# Choose between random line-of-sight response (mode=0) and radial lines
# of sight (mode=1)
if len(sys.argv) == 2:
......@@ -83,7 +83,7 @@ if __name__ == '__main__':
A2 = cfmaker.normalized_amplitudes[1]
DC = SingleDomain(correlated_field.target, position_space)
## Apply a nonlinearity
# Apply a nonlinearity
signal = DC @ ift.sigmoid(correlated_field)
# Build the line-of-sight response and define signal response
......@@ -118,7 +118,7 @@ if __name__ == '__main__':
ic_newton.enable_logging()
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
# Set up likelihood and information Hamiltonian
......@@ -157,7 +157,6 @@ if __name__ == '__main__':
name=filename.format("loop_{:02d}".format(i)))
# Done, draw posterior samples
Nsamples = 20
KL = ift.MetricGaussianKL(mean, H, N_samples)
sc = ift.StatCalculator()
scA1 = ift.StatCalculator()
......@@ -189,3 +188,7 @@ if __name__ == '__main__':
linewidth=[1.]*len(powers2) + [3., 3.])
plot.output(ny=2, nx=2, xsize=15, ysize=15, name=filename_res)
print("Saved results as '{}'.".format(filename_res))
if __name__ == '__main__':
main()
......@@ -34,7 +34,7 @@ from matplotlib.colors import LogNorm
import nifty7 as ift
if __name__ == '__main__':
def main():
dom = ift.UnstructuredDomain(1)
scale = 10
......@@ -119,3 +119,7 @@ if __name__ == '__main__':
ift.logger.info('Finished')
# Uncomment the following line in order to leave the plots open
# plt.show()
if __name__ == '__main__':
main()
......@@ -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) 2019 Max-Planck-Society
# Copyright(C) 2019-2020 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -37,55 +37,52 @@ def convtest(test_signal, delta, func):
# Plot results
plot = ift.Plot()
plot.add(signal, title='Signal')
plot.add(test_signal, title='Signal')
plot.add(conv_signal, title='Signal Convolved')
plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
plot.add(conv_delta, title='Kernel')
plot.output()
# Healpix test
nside = 64
npix = 12 * nside * nside
def main():
# Healpix test
nside = 64
npix = 12 * nside * nside
domain = ift.HPSpace(nside)
domain = ift.HPSpace(nside)
# Define test signal (some point sources)
signal_vals = np.zeros(npix, dtype=np.float64)
for i in range(0, npix, npix//12 + 27):
signal_vals[i] = 500.
# Define test signal (some point sources)
signal_vals = np.zeros(npix, dtype=np.float64)
for i in range(0, npix, npix//12 + 27):
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[0] = 1.0
delta = ift.makeField(domain, delta_vals)
delta_vals = np.zeros(npix, dtype=np.float64)
delta_vals[0] = 1.0
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
def func(theta):
ct = np.cos(theta)
return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
convtest(signal, delta, func)
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))
# 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, lambda d: 1. * np.logical_and(d > 0.1, d <= 0.2))
# 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)
# Define kernel function
def func(dist):
return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
convtest(signal, delta, func)
if __name__ == '__main__':
main()
......@@ -11,7 +11,8 @@
# 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
# Copyright(C) 2013-2020 Max-Planck-Society
# Author: Philipp Arras
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
......@@ -56,7 +57,7 @@ class PolynomialResponse(ift.LinearOperator):
def __init__(self, domain, sampling_points):
if not (isinstance(domain, ift.UnstructuredDomain)
and isinstance(x, np.ndarray)):
and isinstance(sampling_points, np.ndarray)):
raise TypeError
self._domain = ift.DomainTuple.make(domain)
tgt = ift.UnstructuredDomain(sampling_points.shape)
......@@ -80,7 +81,7 @@ class PolynomialResponse(ift.LinearOperator):
return ift.makeField(self._tgt(mode), out)
if __name__ == '__main__':
def main():
# Generate some mock data
N_params = 10
N_samples = 100
......@@ -140,3 +141,7 @@ if __name__ == '__main__':
sigma = np.sqrt(sc.var.val)
for ii in range(len(mean)):
print('Coefficient x**{}: {:.2E} +/- {:.2E}'.format(ii, mean[ii], sigma[ii]))
if __name__ == '__main__':
main()
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