diff --git a/config.py b/config.py
index 29a052a2b886de9eeeb12ff2bebec868302ee011..a53493f64da5e7045cf7b416522b43e9e33b325c 100644
--- a/config.py
+++ b/config.py
@@ -75,119 +75,6 @@ timeparams = {
'prefix': 'time'
}
-zm_single = 0, 1, 0.3, ''
-spaceparams_single = {
- 'target_subdomain': doms,
- 'fluctuations_mean': 0.75,
- 'fluctuations_stddev': 0.3,
- 'flexibility_mean': 0.001,
- 'flexibility_stddev': 0.0001,
- 'asperity_mean': 0.0001,
- 'asperity_stddev': 0.00001,
- 'loglogavgslope_mean': -1.5,
- 'loglogavgslope_stddev': 0.5,
- 'prefix': 'space'
-}
-
-cfm_single = ift.CorrelatedFieldMaker.make(*zm_single)
-cfm_single.add_fluctuations(**spaceparams_single)
-
-logsky = cfm_single.finalize(prior_info=0)
-
-
-sky_single = vlbi.normalize(doms) @ logsky.exp()
-def smoothed_sky_model(sigma):
- cfm_single = ift.CorrelatedFieldMaker.make(*zm_single)
- cfm_single.add_fluctuations(**spaceparams_single)
-
- logsky = cfm_single.finalize(prior_info=0)
-
- SMO = ift.HarmonicSmoothingOperator(logsky.target,sigma*vlbi.MUAS2RAD)
-
- sky_single = vlbi.normalize(doms) @ SMO @ (SMO @ logsky).exp()
- return sky_single
-
-def excitation_smoother(sigma,sigma_old, xi):
- codomain = xi.domain[0]
- kernel = codomain.get_k_length_array()
- old_sig = codomain.get_fft_smoothing_kernel_function(sigma_old*vlbi.MUAS2RAD)
- new_sig = codomain.get_fft_smoothing_kernel_function(sigma*vlbi.MUAS2RAD)
- weights = old_sig(kernel).clip(1e-200,1e300)/new_sig(kernel).clip(1e-200,1e300)
-
- return weights * xi
-
-def smoothed_movie_model(sigma_x,sigma_t):
-
- domt = ift.RGSpace(npixt, dt)
- dom = ift.makeDomain((domt, doms))
- domt_zeropadded = ift.RGSpace(int(oversampling_factor * npixt), dt)
-
- cfm = ift.CorrelatedFieldMaker.make(*zm)
- cfm.add_fluctuations(domt_zeropadded, **timeparams)
- cfm.add_fluctuations(**spaceparams)
- 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
-
- 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)
-
- ls_lo = (FA_hi + FA_lo)
- ls_hi = (FA_hi - FA_lo)
- t_SMO = ift.HarmonicSmoothingOperator(ls_lo.target,sigma_t,space=0)
- x_SMO = ift.HarmonicSmoothingOperator(ls_lo.target,sigma_x*vlbi.MUAS2RAD,space=1)
-
- SMO = t_SMO @ x_SMO
- logsky_lo_smo = SMO @ ls_lo
- logsky_hi_smo = SMO @ ls_hi
-
- sky_lo_smo = SMO @ logsky_lo_smo.exp()
- sky_hi_smo = SMO @ logsky_hi_smo.exp()
-
- sky_lo = FA_lo.adjoint @ sky_lo_smo
- sky_hi = FA_hi.adjoint @ sky_hi_smo
-
- sky_mf = (sky_hi + sky_lo) @ (logsky_1 + logsky_2)
- smooth_movie = vlbi.normalize(logsky_mf.target) @ sky_mf
- return smooth_movie
-
-def movie_excitation_smoother(sigma_x, sigma_t, sigma_x_old, sigma_t_old, xi_hi,xi_lo):
-
- def get_kernel(sigma_x, sigma_t,domain):
-
- def get_single_kernel(domain, sigma):
- k = domain.get_k_length_array()
- kernel = domain.get_fft_smoothing_kernel_function(sigma)
- sig = kernel(k)
- return sig
-
- sig_t = get_single_kernel(domain[0],sigma_t).val
- sig_x = get_single_kernel(domain[1],sigma_x*vlbi.MUAS2RAD).val
- sig = np.outer(sig_t, sig_x)
- sig = sig.reshape(sig_t.shape + sig_x.shape )
- sig = ift.makeField(domain, sig)
- return sig
- sig_old = get_kernel(sigma_x_old, sigma_t_old, xi_hi.domain)
- sig_new = get_kernel(sigma_x, sigma_t, xi_hi.domain)
-
- weights = sig_old.clip(1e-200,1e300)/sig_new.clip(1e-200,1e300)
-
- return xi_hi * weights, xi_lo * weights
-
-
-
-
-
-
-
-pspec_single = cfm_single.amplitude
domt = ift.RGSpace(npixt, dt)
dom = ift.makeDomain((domt, doms))
diff --git a/movie_start.py b/movie_start.py
index 27a7544c485ea987717adb040d7528d0f3ee67e4..c5a58e88458e51c860a52cc956790744e68415b4 100644
--- a/movie_start.py
+++ b/movie_start.py
@@ -21,8 +21,7 @@ import nifty6 as ift
import src as vlbi
from config import comm, rank, master
from config import oversampling_factor as ofac
-from config import sky_movie_mf
-from config import sky_movie_mf as full_sky
+from config import sky_movie_mf as sky
@@ -34,8 +33,8 @@ def main():
p = ift.Plot()
n = 5
for _ in range(5): # FIXME: should this be 'range(n)'?
- pos = ift.from_random(sky_movie_mf.domain, 'normal')
- ss = sky_movie_mf(pos)['hi']
+ pos = ift.from_random(sky.domain, 'normal')
+ ss = sky(pos)['hi']
mi, ma = 0, np.max(ss.val)
for ii in range(0, ss.shape[0], 4):
extr = ift.DomainTupleFieldInserter(ss.domain, 0, (ii,)).adjoint
@@ -44,10 +43,9 @@ def main():
print('Start inf check')
for ii in range(20):
- pp = ift.from_random(sky_movie_mf.domain, 'normal') #* 2.5
- sky_movie_mf(pp)
+ pp = ift.from_random(sky.domain, 'normal') #* 2.5
+ sky(pp)
- sky = full_sky
fld = vlbi.gaussian_profile(sky.target['hi'][1], 30 *vlbi.MUAS2RAD)
fld = ift.ContractionOperator(sky.target['hi'], 0).adjoint(fld).val
multi = ift.makeField(sky.target, {'lo': fld, 'hi': fld})
diff --git a/src/sugar.py b/src/sugar.py
index de0dafde1601881960f14c8b06ad6c7a76d8a9c6..9782d80609d414a9f7ffd8ad1f5e4cbb32605370 100644
--- a/src/sugar.py
+++ b/src/sugar.py
@@ -167,17 +167,6 @@ def gaussian_profile(dom, rad):
profile = 1/(2*np.pi*rad**2)*np.exp(-0.5*(xx**2 + yy**2)/rad**2)
return ift.makeField(dom, profile)
-#FIXME unused
-def lognormal_to_normal(mean, sig):
- tmp = np.log((sig/mean)**2 + 1)
- return np.log(mean)-0.5*tmp, np.sqrt(tmp)
-
-#FIXME unused
-def normal_to_lognormal(mu, sig):
- tmp = np.exp(sig**2)
- logmu = np.exp(mu) * np.sqrt(tmp)
- return logmu, logmu*np.sqrt(tmp-1)
-
class DomainTuple2MultiField(ift.LinearOperator):
def __init__(self, domain, active_inds):
@@ -207,28 +196,6 @@ class DomainTuple2MultiField(ift.LinearOperator):
res = {k: ift.Field.from_raw(self.dom0, v) for k, v in res.items()}
return ift.MultiField.from_dict(res)
-#FIXME unused
-class DomainTuple2MultiField_sf(ift.LinearOperator):
- def __init__(self, domain, active_inds):
- self._domain = ift.DomainTuple.make(domain)
- tgt = {str(ii): self._domain[1:] for ii in set(active_inds)}
- self._target = ift.makeDomain(tgt)
- self._capability = self.TIMES | self.ADJOINT_TIMES
-
- def apply(self, x, mode):
- self._check_input(x, mode)
- x = x.val
- if mode == self.TIMES:
- return ift.MultiField.from_dict({
- ii: ift.makeField(self._domain[1:], x[int(ii)])
- for ii in self._target.keys()
- })
- res = np.zeros(self._domain.shape,
- dtype=x[list(self._target.keys())[0]].dtype)
- for ii in self._target.keys():
- res[int(ii)] = x[ii]
- return ift.Field.from_raw(self._domain, res)
-
def freq_avg(fld):
return 0.5*(fld['hi'] + fld['lo'])
diff --git a/test.py b/test.py
index 6e9b314cdf92009999fe65c2e76862a6826007c3..7be4bb9993151ed733b3f2de6149cfc2e1f6ca13 100644
--- a/test.py
+++ b/test.py
@@ -22,7 +22,7 @@ import pytest
import nifty6 as ift
import src as vlbi
-from config import sky_movie_mf, sky_single
+from config import sky_movie_mf
def setup_function():
@@ -103,7 +103,7 @@ def test_normalize(seed, shape):
np.testing.assert_allclose(s_normalized.s_integrate(), 1.)
-@pmp('sky', [sky_single, sky_movie_mf])
+@pmp('sky', [sky_movie_mf])
def test_normalize2(sky):
fld = sky(ift.from_random(sky.domain, 'normal'))
if isinstance(fld, ift.Field):
@@ -264,10 +264,6 @@ def test_DomainTupleMultiField(ddtype, tdtype):
active_inds = ["0_lo", "4_hi"]
foo = vlbi.DomainTuple2MultiField(dom, active_inds)
ift.extra.consistency_check(foo, domain_dtype=ddtype, target_dtype=tdtype)
- dom = ift.RGSpace(5, 3)
- active_inds = [0, 4]
- foo = vlbi.DomainTuple2MultiField_sf(dom, active_inds)
- ift.extra.consistency_check(foo, domain_dtype=ddtype, target_dtype=tdtype)
@pmp('seed', seeds)
@@ -299,14 +295,6 @@ def test_position_states(op, pre, iteration):
ift.extra.assert_allclose(ss0, ss1, 0, 1e-7)
-@pmp('mean', [0, 1.2])
-@pmp('sig', [1, 1.3])
-def test_normal_lognormal(mean, sig):
- mean1, sig1 = vlbi.lognormal_to_normal(*vlbi.normal_to_lognormal(mean, sig))
- np.testing.assert_allclose(mean, mean1)
- np.testing.assert_allclose(sig, sig1)
-
-
def _std_check(arr):
tol = np.std(arr) if arr.size > 3 else 0.5
np.testing.assert_allclose(np.mean(arr), 1, atol=tol)