Scheduled maintenance on Monday 2019-06-24 between 10:00-11:00 CEST

...
 
Commits (7)
Subproject commit 1933591786679abc846689661db0b2a6e749c0aa
Subproject commit f42d04fe958511f0cd211c5ff5b58f4f7a02c29c
......@@ -89,9 +89,12 @@ def pickle(obj, fname):
pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)
def load_pickle(f):
def load_pickle(fname):
import pickle
try:
return np.load(f)
with open(fname, 'rb') as f:
pos = pickle.load(f)
return pos
except FileNotFoundError:
raise FileNotFoundError('Input file not found.')
except OSError:
......
......@@ -8,8 +8,8 @@ import nifty5 as ift
import resolve as rve
from problem import Problem
from utils import (paper_plot_antenna, paper_plot_calib_solutions,
plot_amplitude_op, paper_plot_pspec, paper_plot_timeline,
paper_plot_uv, plot_integrated, plot_sky)
paper_plot_pspec, paper_plot_timeline, paper_plot_uv,
plot_amplitude_op, plot_integrated, plot_sky, plot_weighted)
def myplots(directory, cglimit, nsamps, position):
......@@ -30,7 +30,8 @@ def myplots(directory, cglimit, nsamps, position):
if nsamps > 1:
print('Draw samples')
ic = ift.GradientNormController(iteration_limit=cglimit)
res_samples = rve.MetricGaussianKL(position, lh, ic, nsamps).samples
ham = ift.StandardHamiltonian(lh, ic)
res_samples = rve.MetricGaussianKL(position, ham, nsamps).samples
with open(fname, 'wb') as f:
pickle.dump(res_samples, f, pickle.HIGHEST_PROTOCOL)
else:
......@@ -77,12 +78,11 @@ def myplots(directory, cglimit, nsamps, position):
res = abs(gt - sc.mean)
plot_sky(res, directory + 'sky_residual', disable_axes, cblabel)
res = gt - sc.mean
plot_weighted(directory + 'weighted', res/sc.var.sqrt())
plot_sky(
res/sc.var.sqrt(),
directory + 'sky_resoversd',
disable_axes,
'',
vmin=0)
abs(gt - sc.mean)/sc.var.sqrt(), directory + 'sky_weighted',
disable_axes, cblabel)
n = res.domain.size
onesig = np.sum((res < sc.var.sqrt()).to_global_data())/n
......@@ -146,8 +146,8 @@ if __name__ == '__main__':
np.random.seed(42)
plot_amplitude_op()
np.random.seed(42)
myplots('mock', 400, 100, '20_imaging')
np.random.seed(17)
myplots('mock', 400, 100, '10_imaging')
np.random.seed(42)
myplots('sn', 400, 100, '32_imaging_adjvar')
......@@ -7,7 +7,7 @@ data column = CORRECTED_DATA
number random rows = 30000
spectral window = 0
channel = -1
snr = 1
snr = 20
[Response]
w planes = 1
......@@ -60,15 +60,15 @@ log-zeromode variance alpha = 2
log-zeromode variance q = 2e1
[Minimization pre-calibration]
sampling rounds = 2
n samples = 2
sampling rounds = 3
n samples = 3
cg sampling = 100
tolerance = 1e-8
maxsteps = 10
[Minimization imaging]
sampling rounds = 10
n samples = 20
cg sampling = 400
sampling rounds = 4
n samples = 5
cg sampling = 200
tolerance = 1e-7
maxsteps = 25
\ No newline at end of file
maxsteps = 12
\ No newline at end of file
# 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) 2019 Max-Planck-Society
from os.path import join
import numpy as np
import nifty5 as ift
import resolve as rve
from problem import Problem
def _print_section(s):
print(80*'-')
print(' {}'.format(s))
print(80*'-')
if __name__ == '__main__':
cfg_file = 'mock.cfg'
np.seterr(all='raise', under='warn')
np.random.seed(17)
ift.fft.enable_fftw()
cfg = rve.ConfigurationParser(cfg_file)
plm = Problem(cfg)
dhs = plm.dhs
cal_ops = plm.cal_ops
rs = plm.rs
diffuse = plm.diffuse
lh, lh_c = plm.lh, plm.lh_c
rve.data_plot(join(cfg.out, 'dataplot'), dhs, 1, 3, 0, cal_ops)
plot_overview = lambda position, name: rve.plot_overview(position, join(cfg.out, name), diffuse, None, diffuse, cal_ops, dhs, rs)
try:
plot_overview(plm.MOCK_POSITION, 'mock')
except AttributeError:
pass
class Saver:
def __init__(self, cfg):
self._counter = 0
self._cfg = cfg
def save(self, position, name):
rve.pickle(
position,
join(self._cfg.out, '{:02}_{}.pickle'.format(
self._counter, name)))
plot_overview(position, '{:02}_{}'.format(self._counter, name))
self._counter += 1
def plot_samples(samples, name):
p = ift.Plot()
for samp in samples:
p.add(diffuse.force(samp))
p.output(name=join(cfg.out, name) + '.png', xsize=20, ysize=20)
###########################################################################
# MINIMIZATION
###########################################################################
saver = Saver(cfg)
position = ift.full(lh.domain, 0.)
plot_samples(
[ift.from_random('normal', diffuse.domain) for _ in range(16)],
'priorsamples')
_print_section('Pre-calibration')
for ii in range(cfg.min_pre_calib['sampling rounds']):
minimizer = cfg.min_pre_calib['minimizer']
c = cfg.min_pre_calib
dct = {
'position': position.extract(lh_c.domain),
'ham': ift.StandardHamiltonian(lh_c, c['ic sampling']),
'nsamp': c['n samples'],
}
e = rve.MetricGaussianKL(**dct)
e, _ = minimizer(e)
position = ift.MultiField.union([position, e.position])
saver.save(position, 'precal')
_print_section('Pre-imaging')
for ii in range(4):
_, lh_only_image = lh.simplify_for_constant_input(
position.extract(lh_c.domain))
minimizer = cfg.min_imaging['minimizer']
pe = rve.key_handler.calib + rve.key_handler.zeromodes
c = cfg.min_imaging
dct = {
'ham': ift.StandardHamiltonian(lh_only_image, c['ic sampling']),
'point_estimates': pe,
'nsamp': 3,
'constants': rve.key_handler.calib
}
e = rve.MetricGaussianKL(position, **dct)
plot_samples([samp + position for samp in e.samples],
'samples_{}'.format(ii))
e, _ = minimizer(e)
position = e.position
saver.save(position, 'preimaging')
# Adjust variances
dct['nsamp'] = 3
e = rve.MetricGaussianKL(position, **dct)
plot_samples([samp + position for samp in e.samples],
'samples_adjvar_{}'.format(ii))
position = diffuse.adjust_variances(position, minimizer, e.samples)
saver.save(position, 'preimaging_adjvar')
# Main imaging
_print_section('Joint calibration and imaging')
for ii in range(cfg.min_imaging['sampling rounds']):
minimizer = cfg.min_imaging['minimizer']
c = cfg.min_imaging
dct = {
'position': position,
'ham': ift.StandardHamiltonian(lh, c['ic sampling']),
'nsamp': c['n samples'],
}
e = rve.MetricGaussianKL(**dct)
plot_samples([samp + position for samp in e.samples],
'samples_{}'.format(ii))
e, _ = minimizer(e)
position = e.position
saver.save(position, 'imaging')
# Adjust variances
e = rve.MetricGaussianKL(**dct)
plot_samples([samp + position for samp in e.samples],
'samples_adjvar_{}'.format(ii))
position = diffuse.adjust_variances(position, minimizer, e.samples)
saver.save(position, 'imaging_adjvar')
......@@ -28,7 +28,8 @@ def _print_section(s):
print(80*'-')
def resolve(cfg_file, preimaging):
if __name__ == '__main__':
cfg_file = 'sn.cfg'
np.seterr(all='raise', under='warn')
np.random.seed(42)
ift.fft.enable_fftw()
......@@ -95,7 +96,7 @@ def resolve(cfg_file, preimaging):
saver.save(position, 'precal')
_print_section('Pre-imaging')
for ii in range(preimaging):
for ii in range(5):
_, lh_only_image = lh.simplify_for_constant_input(
position.extract(lh_c.domain))
minimizer = cfg.min_imaging['minimizer']
......@@ -145,8 +146,3 @@ def resolve(cfg_file, preimaging):
'samples_adjvar_{}'.format(ii))
position = diffuse.adjust_variances(position, minimizer, e.samples)
saver.save(position, 'imaging_adjvar')
if __name__ == '__main__':
resolve('mock.cfg', 1)
resolve('sn.cfg', 5)
......@@ -291,13 +291,13 @@ def paper_plot_pspec(position, name, op, res_samples, ground_truth):
linewidth=0.5)
plt.plot(xcoord, sc.mean.to_global_data(), color=c_mean, linewidth=1)
if ground_truth is not None:
from power_spectrum import pspec
from powerspectrum import pspec
plt.plot(xcoord, pspec(xcoord), color=c_truth, linewidth=1)
plt.xlabel('Spatial frequency [1/rad]')
plt.xlim([xcoord[1], np.max(xcoord)])
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r'Power $[\mathrm{Jy}]$')
plt.ylabel(r'Power $[1]$')
plt.savefig('{}.{}'.format(name, savefig_cfg['format']), **savefig_cfg)
plt.close()
......@@ -373,6 +373,26 @@ def plot_integrated(position, name, diffuse, ground_truth, samples, xmi, xma,
plt.close()
def plot_weighted(name, weighted_diff):
plt.figure(figsize=set_size())
lim = 6
plt.xlim([-lim, lim])
plt.hist(
weighted_diff.to_global_data().flatten(),
bins=30,
density=True)
xs = np.linspace(-lim, lim)
ys = np.exp(-0.5*xs**2)/np.sqrt(2*np.pi)
plt.plot(xs, ys, color=c_truth)
plt.xlabel('Normalized error [1]')
plt.ylabel(r'\# pixels/bin')
plt.savefig('{}.{}'.format(name, savefig_cfg['format']), **savefig_cfg)
plt.close()
def plot_amplitude_op():
a, k0 = 4., 1.
slm, ym = -2.5, 10
......