Skip to content
Snippets Groups Projects
Commit 0a61a3de authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge

parents 151c4e69 2547e3d0
No related branches found
No related tags found
No related merge requests found
......@@ -41,11 +41,12 @@ if __name__ == '__main__':
power_space = A.target[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
dummy = ift.Field.from_random('normal', harmonic_space)
domain = ift.MultiDomain.union(
(A.domain, ift.MultiDomain.make({'xi': harmonic_space})))
domain = ift.MultiDomain.union((A.domain,
ift.MultiDomain.make({
'xi': harmonic_space
})))
correlated_field = ht(
power_distributor(A)*ift.FieldAdapter(domain, "xi"))
correlated_field = ht(power_distributor(A)*ift.FieldAdapter(domain, "xi"))
# alternatively to the block above one can do:
# correlated_field = ift.CorrelatedField(position_space, A)
......@@ -54,8 +55,7 @@ if __name__ == '__main__':
# Building the Line of Sight response
LOS_starts, LOS_ends = get_random_LOS(100)
R = ift.LOSResponse(position_space, starts=LOS_starts,
ends=LOS_ends)
R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
# build signal response model and model likelihood
signal_response = R(signal)
# specify noise
......@@ -68,8 +68,7 @@ if __name__ == '__main__':
data = signal_response(MOCK_POSITION) + N.draw_sample()
# set up model likelihood
likelihood = ift.GaussianEnergy(
mean=data, covariance=N)(signal_response)
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
# set up minimization and inversion schemes
ic_sampling = ift.GradientNormController(iteration_limit=100)
......@@ -84,17 +83,18 @@ if __name__ == '__main__':
position = INITIAL_POSITION
plot = ift.Plot()
plot.add(signal(MOCK_POSITION), title='ground truth')
plot.add(R.adjoint_times(data), title='data')
plot.add([A(MOCK_POSITION)], title='power')
plot.output(nx=3, xsize=16, ysize=5, title="setup", name="setup.png")
plot.add(signal(MOCK_POSITION), title='Ground Truth')
plot.add(R.adjoint_times(data), title='Data')
plot.add([A(MOCK_POSITION)], title='Power Spectrum')
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
# number of samples used to estimate the KL
N_samples = 20
for i in range(2):
metric = H(ift.Linearization.make_var(position)).metric
samples = [metric.draw_sample(from_inverse=True)
for _ in range(N_samples)]
samples = [
metric.draw_sample(from_inverse=True) for _ in range(N_samples)
]
KL = ift.SampledKullbachLeiblerDivergence(H, samples)
KL = ift.EnergyAdapter(position, KL)
......@@ -104,15 +104,17 @@ if __name__ == '__main__':
plot = ift.Plot()
plot.add(signal(position), title="reconstruction")
plot.add([A(position), A(MOCK_POSITION)], title="power")
plot.output(nx=2, xsize=12, ysize=6, title="loop", name="loop.png")
plot.output(ny=1, ysize=6, xsize=16, name="loop.png")
plot = ift.Plot()
sc = ift.StatCalculator()
for sample in samples:
sc.add(signal(sample+position))
plot.add(sc.mean, title="mean")
plot.add(ift.sqrt(sc.var), title="std deviation")
plot.add(sc.mean, title="Posterior Mean")
plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
powers = [A(s+position) for s in samples]
plot.add([A(position), A(MOCK_POSITION)]+powers, title="power")
plot.output(nx=3, xsize=16, ysize=5, title="results", name="results.png")
plot.add(
[A(position), A(MOCK_POSITION)]+powers,
title="Sampled Posterior Power Spectrum")
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment