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

Simplifications

parent 6fdbf96c
......@@ -15,40 +15,23 @@
# Author: Philipp Arras, Philipp Frank, Philipp Haim, Reimar Leike,
# Jakob Knollmueller
import os
import nifty6 as ift
import util as vres
oversampling_factor = 2
nthreads = int(os.getenv('OMP_NUM_THREADS', 1))
ift.fft.set_nthreads(nthreads)
ift.fft.set_nthreads(1)
eps = 1e-7
npix = 128
dt = 6 # hours
fov = 256*vres.MUAS2RAD # RAD
fov = 256*vres.MUAS2RAD
doms = ift.RGSpace(2*(npix,), 2*(fov/npix,))
npixt = 7*24//dt
startt = 0
if dt is not None:
dt = int(dt)
_fm, _fs, _flm, _fls = 1.5, 1., 0.3, 0.001
_am, _as, _lm, _ls = 0.01, 0.001, -1.5, 0.5
if npix is not None and dt is not None:
npix = int(npix)
doms = ift.RGSpace(2*(npix,), 2*(fov/npix,))
npixt = 7*24//dt
startt = 0
domt = ift.RGSpace(npixt, dt)
dom = ift.makeDomain((domt, doms))
print('Sky domain:')
print(dom)
domt_zeropadded = ift.RGSpace(int(oversampling_factor*npixt), dt)
cfm = ift.CorrelatedFieldMaker.make(0.2, 1e-1, '')
cfm.add_fluctuations(domt_zeropadded,
cfm = ift.CorrelatedFieldMaker.make(0.2, 1e-1, '')
cfm.add_fluctuations(ift.RGSpace(int(2*npixt), dt),
fluctuations_mean=.2,
fluctuations_stddev=1.,
flexibility_mean=0.1,
......@@ -58,32 +41,30 @@ if npix is not None and dt is not None:
loglogavgslope_mean=-4,
loglogavgslope_stddev=0.5,
prefix='time')
cfm.add_fluctuations(doms,
fluctuations_mean=cfm.moment_slice_to_average(_fm),
fluctuations_stddev=_fs,
flexibility_mean=_flm,
flexibility_stddev=_fls,
asperity_mean=_am,
asperity_stddev=_as,
loglogavgslope_mean=_lm,
loglogavgslope_stddev=_ls,
cfm.add_fluctuations(doms,
fluctuations_mean=cfm.moment_slice_to_average(1.5),
fluctuations_stddev=1.,
flexibility_mean=0.3,
flexibility_stddev=0.001,
asperity_mean=0.01,
asperity_stddev=0.001,
loglogavgslope_mean=-1.5,
loglogavgslope_stddev=0.5,
prefix='space')
logsky = cfm.finalize(prior_info=0)
logsky = cfm.finalize(prior_info=0)
FA_lo = ift.FieldAdapter(logsky.target, 'lo')
FA_hi = ift.FieldAdapter(logsky.target, 'hi')
xi_hi = ift.FieldAdapter(logsky.domain['xi'], 'xi_hi')
id_hi = ift.FieldAdapter(logsky.domain['xi'], 'xi').adjoint @ xi_hi
xi_lo = ift.FieldAdapter(logsky.domain['xi'], 'xi_lo')
id_lo = ift.FieldAdapter(logsky.domain['xi'], 'xi').adjoint @ xi_lo
FA_lo = ift.FieldAdapter(logsky.target, 'lo')
FA_hi = ift.FieldAdapter(logsky.target, 'hi')
xi_hi = ift.FieldAdapter(logsky.domain['xi'], 'xi_hi')
id_hi = ift.FieldAdapter(logsky.domain['xi'], 'xi').adjoint @ xi_hi
xi_lo = ift.FieldAdapter(logsky.domain['xi'], 'xi_lo')
id_lo = ift.FieldAdapter(logsky.domain['xi'], 'xi').adjoint @ xi_lo
expected_difference = 1e-2
logsky_1 = (FA_hi.adjoint @ logsky).partial_insert(id_hi)
logsky_2 = (FA_lo.adjoint
@ logsky).partial_insert(id_lo).scale(expected_difference)
logsky_lo = FA_lo.adjoint @ (FA_hi + FA_lo)
logsky_hi = FA_hi.adjoint @ (FA_hi - FA_lo)
logsky_mf = (logsky_hi + logsky_lo) @ (logsky_1 + logsky_2)
logsky_1 = (FA_hi.adjoint @ logsky).partial_insert(id_hi)
logsky_2 = (FA_lo.adjoint @ logsky).partial_insert(id_lo).scale(0.01)
logsky_lo = FA_lo.adjoint @ (FA_hi + FA_lo)
logsky_hi = FA_hi.adjoint @ (FA_hi - FA_lo)
logsky_mf = (logsky_hi + logsky_lo) @ (logsky_1 + logsky_2)
sky_movie_mf = vres.normalize(logsky_mf.target) @ logsky_mf.exp()
pspec_movie = cfm.normalized_amplitudes
sky_movie_mf = vres.normalize(logsky_mf.target) @ logsky_mf.exp()
pspec_movie = cfm.normalized_amplitudes
......@@ -21,13 +21,12 @@ os.environ["OMP_NUM_THREADS"] = "1"
import sys
from functools import reduce
from operator import add
from time import time
import numpy as np
import nifty6 as ift
import util as vres
from config import doms, dt, eps, npixt, nthreads
from config import doms, dt, eps, npixt
from config import sky_movie_mf as sky
from config import startt
......@@ -40,42 +39,35 @@ if __name__ == '__main__':
raise RuntimeError
path = f'data/{pre}'
lh = []
active_inds = []
# Build data model
lh, active_inds = [], []
for freq in vres.data.FREQS:
rawd = vres.combined_data(path, [freq], identify_short_baselines=['APAA', 'SMJC'])
args = {'tmin': startt, 'tmax': npixt*dt, 'delta_t': dt}
for ii, dd in enumerate(vres.time_binning(rawd, **args)):
for ii, dd in enumerate(vres.time_binning(rawd, dt, startt, npixt*dt)):
if len(dd) == 0:
continue
ind = str(ii) + "_" + freq
active_inds.append(ind)
v2ph, icovph = vres.Visibilities2ClosurePhases(dd)
v2ampl, icovampl = vres.Visibilities2ClosureAmplitudes(dd)
R = vres.RadioResponse(doms, dd['uv'], eps).ducktape(ind)
vis = ift.makeField(R.target, dd['vis'])
ll_ph = ift.GaussianEnergy(v2ph(vis), icovph) @ v2ph
ll_ampl = ift.GaussianEnergy(v2ampl(vis), icovampl) @ v2ampl
lh.append((ll_ph + ll_ampl) @ R)
lh = reduce(add, lh) @ vres.DomainTuple2MultiField(sky.target, active_inds)
pb = vres.gaussian_profile(sky.target['hi'][1], 50*vres.MUAS2RAD)
pb = ift.ContractionOperator(sky.target['hi'], 0).adjoint(pb)
pb = ift.makeOp(ift.MultiField.from_dict({'hi': pb, 'lo': pb}))
vis2closph, invcov_ph = vres.Visibilities2ClosurePhases(dd)
vis2closampl, invcov_ampl = vres.Visibilities2ClosureAmplitudes(dd)
nfft = vres.RadioResponse(doms, dd['uv'], eps, nthreads)
vis = ift.makeField(nfft.target, dd['vis'])
ll = ift.GaussianEnergy(mean=vis2closph(vis), inverse_covariance=invcov_ph) @ vis2closph
ll = ll + ift.GaussianEnergy(mean=vis2closampl(vis), inverse_covariance=invcov_ampl) @ vis2closampl
lh.append(ll @ nfft.ducktape(ind))
conv = vres.DomainTuple2MultiField(sky.target, active_inds)
lh = reduce(add, lh) @ conv
# Plot prior samples
for ii in range(4):
pp = ift.from_random('normal', sky.domain)
vres.save_state(sky, pp, pre, f'prior{ii}', samples=[0*pp, 0*pp])
vres.save_state(sky, pp, pre, f'prior{ii}', [])
# Minimization
pos = 0.1*ift.from_random('normal', sky.domain)
vres.save_state(sky, pos, pre, 'initial', samples=[0*pos, 0*pos])
ic = ift.GradientNormController(iteration_limit=40, name='Sampling')
pb = vres.gaussian_profile(sky.target['hi'][1], 50*vres.MUAS2RAD)
pb = ift.ContractionOperator(sky.target['hi'], 0).adjoint(pb)
pb = ift.MultiField.from_dict({'hi': pb, 'lo': pb})
t0 = time()
(lh @ ift.makeOp(pb) @ sky)(pos)
print(f'Likelihood call: {1000*(time()-t0):.0f} ms')
vres.save_state(sky, pos, pre, 'initial', [])
for ii in range(30):
if ii < 10:
N_samples = 1
......@@ -93,23 +85,16 @@ if __name__ == '__main__':
N_samples = 20
N_iterations = 30
print(f'Iter: {ii}, N_samples: {N_samples}, N_iter: {N_iterations}')
if ii < 15:
print('Primary beam hack active')
ll = lh @ ift.makeOp(pb) @ sky
else:
print('No primary beam anymore')
ll = lh @ sky
cstkeys = ('spaceasperity', 'spaceflexibility', 'spacefluctuations', 'spaceloglogavgslope', 'spacespectrum', 'timeasperity', 'timeflexibility', 'timefluctuations', 'timeloglogavgslope', 'timespectrum', 'zeromode')
pe = cstkeys if ii < 20 else []
cstkeys = set(pos.domain.keys()) - set(['xi_lo', 'xi_hi'])
cst = cstkeys if ii < 20 else []
H = ift.StandardHamiltonian(ll, ic)
H = ift.StandardHamiltonian(lh @ pb @ sky if ii < 15 else lh @ sky,
ift.GradientNormController(iteration_limit=40))
KL = ift.MetricGaussianKL(pos, H, N_samples, mirror_samples=True,
lh_sampling_dtype=np.complex128, point_estimates=pe, constants=cst)
ic_newton = ift.AbsDeltaEnergyController(0.1,
iteration_limit=N_iterations,
name=f'it_{ii}',
convergence_level=2)
minimizer = ift.NewtonCG(ic_newton)
lh_sampling_dtype=np.complex128,
point_estimates=cst, constants=cst)
dct = {'deltaE': 0.1, 'iteration_limit': N_iterations,
'name': f'it_{ii}', 'convergence_level': 2}
minimizer = ift.NewtonCG(ift.AbsDeltaEnergyController(**dct))
KL, _ = minimizer(KL)
pos = KL.position
vres.save_state(sky, pos, pre, ii, samples=KL.samples)
vres.save_state(sky, pos, pre, ii, KL.samples)
......@@ -23,10 +23,8 @@ import nifty6 as ift
import util as eht
pmp = pytest.mark.parametrize
bools = [False, True]
# cb = [[], ['APAA'], ['SMJC'], ['APAA', 'SMJC']]
cb = [['APAA', 'SMJC']]
d3 = {
cb = [[], ['APAA', 'SMJC']]
ds = [{
'uv': np.random.randn(11, 2)*1e7,
'vis': np.random.randn(11) + 1j*np.random.rand(11),
'sigma': np.random.randn(11)**2,
......@@ -34,12 +32,9 @@ d3 = {
'time': np.array(5*[7.24027777] + 6*[7.4236114]),
'ant1': np.array([0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 2]),
'ant2': np.array([1, 2, 3, 2, 3, 1, 2, 3, 2, 3, 3])
}
ds = [
eht.read_data(dd, ff, o2)
for dd, ff, o2 in product(eht.DAYS, eht.FREQS, cb)
] + [eht.combined_data(eht.FREQS, identify_short_baselines=o2)
for o2 in cb] + [d3]
}]
ds += [eht.read_data('data/disk', dd, ff, o2) for dd, ff, o2 in product(eht.DAYS, eht.FREQS, cb)]
ds += [eht.combined_data('data/disk', eht.FREQS, identify_short_baselines=o2) for o2 in cb]
seeds = [189, 4938]
dtypes = [np.complex128, np.float64]
......@@ -82,8 +77,9 @@ def test_visibility_closure_matrices(n):
Psi = eht.closure.closure_phase_design_matrix(n)
np.testing.assert_allclose(Psi @ Phi, 0)
@pmp('seed', seeds)
@pmp('shape', [(3,4), (10,20)])
@pmp('shape', [(3, 4), (10, 20)])
def test_normalize(seed, shape):
np.random.seed(seed)
sp = ift.RGSpace(shape, distances=np.exp(np.random.randn(len(shape))))
......@@ -92,7 +88,6 @@ def test_normalize(seed, shape):
np.testing.assert_allclose(s_normalized.integrate(), 1.)
@pmp('seed', seeds)
@pmp('d', ds)
@pmp('mode', ['ph', 'ampl'])
......@@ -105,16 +100,19 @@ def test_closure_consistency(seed, d, mode):
pos = ift.from_random('normal', op.domain, dtype=np.complex128)
ift.extra.check_jacobian_consistency(op, pos)
@pmp('seed', seeds)
@pmp('shape', [(3,4), (10,20)])
@pmp('shape', [(3, 4), (10, 20)])
def test_sparse_cov(seed, shape):
np.random.seed(seed)
matrix = np.random.randn(*shape)
sigma_sq = np.exp(np.random.randn(shape[1]))
d, r, c = eht.closure.sparse_cov(matrix, sigma_sq)
mat = np.zeros((shape[0],)*2)
mat[(r,c)] = d
np.testing.assert_allclose(mat.T @ mat @ matrix @ np.diag(sigma_sq) @ matrix.T, np.eye(shape[0]), rtol = 1e-7, atol=1e-7)
mat[(r, c)] = d
res0 = mat.T @ mat @ matrix @ np.diag(sigma_sq) @ matrix.T
res1 = np.eye(shape[0])
np.testing.assert_allclose(res0, res1, rtol=1e-7, atol=1e-7)
@pmp('seed', seeds)
......@@ -179,24 +177,18 @@ def test_nfft(npix, epsilon, dd):
only_r_linear=True)
def test_closure():
np.random.seed(42)
def vis2closure(d, ind):
def _vis2closure(d, ind):
ind = np.sort(ind)
i, j, k = ind
vis_ind = tuple(
np.where(np.logical_and(d['ant1'] == a, d['ant2'] == b) == 1)[0]
for a, b in ((i, j), (j, k), (i, k)))
phase = d['vis']/abs(d['vis'])
return phase[vis_ind[0]]*phase[vis_ind[1]]*phase[
vis_ind[2]].conjugate()
# Antenna setup
# 0 -- 3 -- 4
# | \/ |
# | /\ |
# 1 -- 2
return phase[vis_ind[0]]*phase[vis_ind[1]]*phase[vis_ind[2]].conjugate()
def test_closure():
np.random.seed(42)
ant1 = []
ant2 = []
for ii in range(4):
......@@ -217,9 +209,9 @@ def test_closure():
d['sigma'] = abs(d['vis'])/np.sqrt(ww)
closure = np.empty(3, dtype=np.complex128)
closure[0] = vis2closure(d, [1, 2, 3])[0]
closure[1] = vis2closure(d, [0, 2, 3])[0]
closure[2] = vis2closure(d, [0, 1, 3])[0]
closure[0] = _vis2closure(d, [1, 2, 3])[0]
closure[1] = _vis2closure(d, [0, 2, 3])[0]
closure[2] = _vis2closure(d, [0, 1, 3])[0]
closure_op = eht.Visibilities2ClosurePhases(d)[0]
vis = ift.makeField(ift.UnstructuredDomain(len(ant1)), d['vis'])
......
......@@ -15,13 +15,8 @@
import numpy as np
DEG2RAD = np.pi/180.
SPEEDOFLIGHT = 299792458.
ARCMIN2RAD = 1/60*DEG2RAD
AS2RAD = 1/3600*DEG2RAD
MUAS2RAD = 1/1e6/3600*DEG2RAD
MUAS2RAD = 1e-6/3600*np.pi/180
MAXSHORTBASELINE = 2000000
KEYS = ['time', 'ant1', 'ant2', 'uv', 'vis', 'sigma']
FREQS = ['lo', 'hi']
......
......@@ -65,17 +65,16 @@ def read_data(prefix, day, freq, identify_short_baselines):
assert day in DAYS
fname = f'{prefix}_{day}_{freq}.csv'
d = read_data_raw(fname)
all_antennas = list(set(d['ant1']) | set(d['ant2']))
ant1, ant2 = d['ant1'], d['ant2']
all_antennas = list(set(ant1) | set(ant2))
if len(identify_short_baselines) > 0:
dct = {aa: ii for ii, aa in enumerate(all_antennas)}
if 'APAA' in identify_short_baselines:
dct['AP'] = dct['AA']
if 'SMJC' in identify_short_baselines:
dct['SM'] = dct['JC']
ant1 = [dct[ele] for ele in d['ant1']]
ant2 = [dct[ele] for ele in d['ant2']]
ant1 = [dct[ele] for ele in ant1]
ant2 = [dct[ele] for ele in ant2]
ddct = defaultdict(lambda: len(ddct))
d['ant1'] = np.array([ddct[ele] for ele in ant1])
d['ant2'] = np.array([ddct[ele] for ele in ant2])
......
......@@ -23,11 +23,10 @@ from .constants import SPEEDOFLIGHT
class RadioResponse(ift.LinearOperator):
def __init__(self, domain, uv, epsilon, j=1):
def __init__(self, domain, uv, epsilon):
ndata = uv.shape[0]
self._uvw = np.zeros((ndata, 3), dtype=uv.dtype)
self._uvw[:, 0:2] = uv
self._j = int(j)
self._eps = float(epsilon)
self._domain = ift.makeDomain(domain)
self._target = ift.makeDomain(ift.UnstructuredDomain(ndata))
......@@ -40,8 +39,8 @@ class RadioResponse(ift.LinearOperator):
nx, ny = self._domain[0].shape
f = np.array([SPEEDOFLIGHT])
if mode == self.TIMES:
res = ng.dirty2ms(self._uvw, f, x, None, dx, dy, self._eps, False,
nthreads=1)
res = ng.dirty2ms(self._uvw, f, x, None, dx, dy, self._eps,
False, nthreads=1)
res = res[:, 0]
else:
x = x[:, None]
......
......@@ -38,25 +38,27 @@ def save_state(sky, pos, pre, name, samples):
dom = sky.target['hi']
dt = dom[0].distances[0]
sc = ift.StatCalculator()
if len(samples) == 0:
samples = 2*(0*pos,)
if len(samples) == 1:
samples.append(samples[0])
for sam in samples:
sk = sky(pos+sam)
sk = 0.5*(sk['hi']+sk['lo'])
sc.add(sk)
vlim = [0., np.max(sc.mean.val)]
rel = sc.var.sqrt()/sc.mean
varlim = [0., np.max(rel.val)]
for i in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(i*24/dt),)).adjoint
for ii in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(ii*24/dt),)).adjoint
p.add(fi(sc.mean), zmin=vlim[0], zmax=vlim[1])
for i in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(i*24/dt),)).adjoint
for ii in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(ii*24/dt),)).adjoint
p.add(fi(rel), zmin=varlim[0], zmax=varlim[1])
for i in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(i*24/dt),)).adjoint
for ii in range(7):
fi = ift.DomainTupleFieldInserter(dom, 0, (int(ii*24/dt),)).adjoint
p.add(fi(sk), zmin=vlim[0], zmax=vlim[1])
p.output(nx=7, ny=3, name=f'{pre}{name}.png', xsize=15, ysize=5)
p.output(nx=7, ny=3, name=f'{pre}{name}.png', xsize=20, ysize=8)
def baselines(alst):
......
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