Skip to content
Snippets Groups Projects
Commit a37d0e7d authored by Julian Rüstig's avatar Julian Rüstig
Browse files

multi_resolve: IRGSpace from binbounds

parent d4b5a253
No related branches found
No related tags found
1 merge request!50Draft: Mosaic imaging
Pipeline #232476 failed
......@@ -28,13 +28,44 @@ class IRGSpace(ift.StructuredDomain):
_needed_for_hash = ["_coordinates"]
def __init__(self, coordinates):
def __init__(self, coordinates, binbounds=None):
bb = np.array(coordinates)
if bb.ndim != 1:
raise TypeError
if np.any(np.diff(bb) <= 0.0):
raise ValueError("Coordinates must be sorted and strictly ascending")
raise ValueError(
"Coordinates must be sorted and strictly ascending")
self._coordinates = tuple(bb)
self._binbounds = binbounds
def binbounds(self):
if self._binbounds is not None:
return self._binbounds
if len(self._coordinates) == 1:
return np.array([-np.inf, np.inf])
c = np.array(self._coordinates)
bounds = np.empty(self.size + 1)
bounds[1:-1] = c[:-1] + 0.5*np.diff(c)
bounds[0] = c[0] - 0.5*(c[1] - c[0])
bounds[-1] = c[-1] + 0.5*(c[-1] - c[-2])
return bounds
@classmethod
def from_binbounds(cls, binbounds: list[float]):
'''Builds a IRGSpace such that the passed binbounds can be returned.'''
binbounds = np.array(binbounds)
if len(binbounds) < 2:
raise ValueError("binbounds must have at least 2 elements")
if np.any(np.diff(binbounds) <= 0):
raise ValueError("binbounds must be strictly ascending")
coordinates = np.array([(binbounds[ii+1]-binbounds[ii])/2+binbounds[ii]
for ii in range(len(binbounds)-1)])
return cls(coordinates, binbounds)
def __repr__(self):
return f"IRGSpace(coordinates={self._coordinates})"
......@@ -61,16 +92,6 @@ class IRGSpace(ift.StructuredDomain):
"""Assume that the coordinates are the center of symmetric pixels."""
return np.diff(self.binbounds())
def binbounds(self):
if len(self._coordinates) == 1:
return np.array([-np.inf, np.inf])
c = np.array(self._coordinates)
bounds = np.empty(self.size + 1)
bounds[1:-1] = c[:-1] + 0.5*np.diff(c)
bounds[0] = c[0] - 0.5*(c[1] - c[0])
bounds[-1] = c[-1] + 0.5*(c[-1] - c[-2])
return bounds
@property
def distances(self):
return np.diff(self._coordinates)
......
......@@ -48,7 +48,8 @@ def sky_model_diffuse(cfg, observations=[], nthreads=1):
op, aa = _multi_freq_logsky_cfm(cfg, sdom, pol_lbl)
elif cfg["freq mode"] == "iwp":
freq = _get_frequencies(cfg, observations)
op, aa = _multi_freq_logsky_integrated_wiener_process(cfg, sdom, pol_lbl, freq)
op, aa = _multi_freq_logsky_integrated_wiener_process(
cfg, sdom, pol_lbl, freq)
else:
raise RuntimeError
logsky[lbl] = op
......@@ -60,7 +61,8 @@ def sky_model_diffuse(cfg, observations=[], nthreads=1):
tgt = default_sky_domain(pdom=pdom, fdom=fdom, sdom=sdom)
logsky = reduce(add, (oo.ducktape_left(lbl) for lbl, oo in logsky.items()))
mexp = polarization_matrix_exponential_mf2f(logsky.target, nthreads=nthreads)
mexp = polarization_matrix_exponential_mf2f(
logsky.target, nthreads=nthreads)
sky = mexp @ logsky
sky = sky.ducktape_left(tgt)
......@@ -84,29 +86,37 @@ def sky_model_points(cfg, observations=[], nthreads=1):
alpha = cfg.getfloat("point sources alpha")
q = cfg.getfloat("point sources q")
inserter = PointInserter(default_sky_domain(pdom=pdom, sdom=sdom), ppos)
inserter = PointInserter(
default_sky_domain(pdom=pdom, sdom=sdom), ppos)
if pdom.labels_eq("I"):
points = ift.InverseGammaOperator(inserter.domain, alpha=alpha, q=q/sdom.scalar_dvol)
points = ift.InverseGammaOperator(
inserter.domain, alpha=alpha, q=q/sdom.scalar_dvol)
points = points.ducktape("points")
elif pdom.labels_eq(["I", "Q", "U"]) or pdom.labels_eq(["I", "Q", "U", "V"]):
points_domain = inserter.domain[-1]
npoints = points_domain.size
i = ift.InverseGammaOperator(points_domain, alpha=alpha, q=q/sdom.scalar_dvol).log().ducktape("points I")
q = ift.NormalTransform(cfg["point sources stokesq log mean"], cfg["point sources stokesq log stddev"], "points Q", npoints)
u = ift.NormalTransform(cfg["point sources stokesu log mean"], cfg["point sources stokesu log stddev"], "points U", npoints)
i = ift.InverseGammaOperator(
points_domain, alpha=alpha, q=q/sdom.scalar_dvol).log().ducktape("points I")
q = ift.NormalTransform(cfg["point sources stokesq log mean"],
cfg["point sources stokesq log stddev"], "points Q", npoints)
u = ift.NormalTransform(cfg["point sources stokesu log mean"],
cfg["point sources stokesu log stddev"], "points U", npoints)
i = i.ducktape_left("I")
q = q.ducktape_left("Q")
u = u.ducktape_left("U")
polsum = i + q + u
if pdom.labels_eq(["I", "Q", "U", "V"]):
v = ift.NormalTransform(cfg["point sources stokesv log mean"], cfg["point sources stokesv log stddev"], "points V", npoints)
v = ift.NormalTransform(cfg["point sources stokesv log mean"],
cfg["point sources stokesv log stddev"], "points V", npoints)
v = v.ducktape_left("V")
polsum = polsum + v
points = polarization_matrix_exponential_mf2f(polsum.target, nthreads=nthreads) @ polsum
points = polarization_matrix_exponential_mf2f(
polsum.target, nthreads=nthreads) @ polsum
points = points.ducktape_left(inserter.domain)
else:
raise NotImplementedError(f"single_frequency_sky does not support point sources on {pdom.labels} (yet?)")
raise NotImplementedError(
f"single_frequency_sky does not support point sources on {pdom.labels} (yet?)")
additional["point_list"] = points
sky = inserter @ points
......@@ -124,7 +134,8 @@ def sky_model_points(cfg, observations=[], nthreads=1):
inserter = PointInserter(sky.target, ppos)
udom = inserter.domain[-1]
p_i0 = ift.InverseGammaOperator(udom, alpha=alpha, q=q/sdom.scalar_dvol)
p_i0 = ift.InverseGammaOperator(
udom, alpha=alpha, q=q/sdom.scalar_dvol)
p_i0 = p_i0.ducktape("points")
p_alpha = _parse_or_none(cfg, "point sources alpha")
......@@ -142,9 +153,11 @@ def sky_model_points(cfg, observations=[], nthreads=1):
log_fdom = IRGSpace(np.sort(np.log(freq)))
nfreq = len(freq)
npoints = udom.size
p_xi = ift.ScalingOperator(ift.UnstructuredDomain(2*npoints*(nfreq - 1)), 1.).ducktape("points xi")
p_xi = ift.ScalingOperator(ift.UnstructuredDomain(
2*npoints*(nfreq - 1)), 1.).ducktape("points xi")
points = _integrated_wiener_process(p_i0, p_alpha, log_fdom, p_flex, p_asp, p_xi)
points = _integrated_wiener_process(
p_i0, p_alpha, log_fdom, p_flex, p_asp, p_xi)
points = points.ducktape_left(inserter.domain)
additional["point_list"] = points
......@@ -179,7 +192,8 @@ def _multi_freq_logsky_cfm(cfg, sdom, pol_label):
fdom = IRGSpace(freq0 + np.arange(fnpix)*df)
fdom_rg = ift.RGSpace(fnpix, df)
cfm = cfm_from_cfg(cfg, {"freq": fdom_rg, "space": sdom}, f"stokes{pol_label} diffuse")
cfm = cfm_from_cfg(
cfg, {"freq": fdom_rg, "space": sdom}, f"stokes{pol_label} diffuse")
op = cfm.finalize(0)
fampl, sampl = list(cfm.get_normalized_amplitudes())
......@@ -219,7 +233,8 @@ def _multi_freq_logsky_integrated_wiener_process(cfg, sdom, pol_label, freq):
flexibility = _parse_or_none(cfg, prefix + " wp flexibility")
if flexibility is None:
raise RuntimeError("freq flexibility cannot be None")
flexibility = ift.LognormalTransform(*flexibility, prefix + " wp flexibility", 0)
flexibility = ift.LognormalTransform(
*flexibility, prefix + " wp flexibility", 0)
asperity = _parse_or_none(cfg, prefix + " wp asperity")
asperity = ift.LognormalTransform(*asperity, prefix + " wp asperity", 0)
......@@ -251,16 +266,20 @@ def _integrated_wiener_process(i0, alpha, irg_space, flexibility, asperity, freq
shift = np.ones(dom.shape)
shift[0] = vol * vol / 12.0
if asperity is None:
shift = ift.DiagonalOperator(ift.makeField(dom, shift).sqrt(), intop.domain, 0)
shift = ift.DiagonalOperator(ift.makeField(
dom, shift).sqrt(), intop.domain, 0)
increments = shift @ (freq_xi * sig_flex)
else:
vasp = np.empty(dom.shape)
vasp[0] = 1
vasp[1] = 0
vasp = ift.DiagonalOperator(ift.makeField(dom, vasp), domain=broadcast.target, spaces=0)
vasp = ift.DiagonalOperator(ift.makeField(
dom, vasp), domain=broadcast.target, spaces=0)
sig_asp = broadcast_full @ vasp @ broadcast @ asperity
shift = ift.makeField(intop.domain, np.broadcast_to(shift[..., None, None], intop.domain.shape))
increments = freq_xi * sig_flex * (ift.Adder(shift) @ sig_asp).ptw("sqrt")
shift = ift.makeField(intop.domain, np.broadcast_to(
shift[..., None, None], intop.domain.shape))
increments = freq_xi * sig_flex * \
(ift.Adder(shift) @ sig_asp).ptw("sqrt")
return IntWProcessInitialConditions(i0, alpha, intop @ increments)
......@@ -268,7 +287,8 @@ def _integrated_wiener_process(i0, alpha, irg_space, flexibility, asperity, freq
def cfm_from_cfg(cfg, domain_dct, prefix, total_N=0, dofdex=None, override={}, domain_prefix=None):
assert len(prefix) > 0
product_spectrum = len(domain_dct) > 1
cfm = ift.CorrelatedFieldMaker(prefix if domain_prefix is None else domain_prefix, total_N=total_N)
cfm = ift.CorrelatedFieldMaker(
prefix if domain_prefix is None else domain_prefix, total_N=total_N)
for key_prefix, dom in domain_dct.items():
ll = _append_to_nonempty_string(key_prefix, " ")
kwargs = {kk: _parse_or_none(cfg, f"{prefix} {ll}{kk}", override)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment