Commit d356d9bb authored by Philipp Frank's avatar Philipp Frank
Browse files

cleanup

parent 0db16464
Pipeline #91846 failed with stages
in 4 minutes and 26 seconds
......@@ -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))
......
......@@ -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})
......
......@@ -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'])
......
......@@ -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)
......
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