Commit 378af2c5 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'work_on_demos' into 'NIFTy_4'

Work on demos

See merge request ift/NIFTy!218
parents cfe6879c 4e062f6b
Pipeline #24402 passed with stage
in 6 minutes and 11 seconds
......@@ -494,12 +494,7 @@
},
"outputs": [],
"source": [
"sc = ift.probing.utils.StatCalculator()\n",
"for i in range(200):\n",
" print i\n",
" sc.add(HT(curv.generate_posterior_sample()))\n",
"\n",
"m_var = sc.var"
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 200)\n"
]
},
{
......@@ -568,7 +563,7 @@
"outputs": [],
"source": [
"fig = plt.figure(figsize=(15,10))\n",
"plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=1)\n",
"plt.plot(s_data, 'g', label=\"Signal\", alpha=1, linewidth=4)\n",
"plt.plot(d_data, 'k+', label=\"Data\", alpha=.5)\n",
"plt.plot(m_data, 'r', label=\"Reconstruction\")\n",
"plt.axvspan(l, h, facecolor='0.8', alpha=.5)\n",
......@@ -649,18 +644,7 @@
"m = D(j)\n",
"\n",
"# Uncertainty\n",
"sc = ift.probing.utils.StatCalculator()\n",
"\n",
"IC = ift.GradientNormController(iteration_limit=50000,\n",
" tol_abs_gradnorm=0.1)\n",
"inverter = ift.ConjugateGradient(controller=IC)\n",
"curv = ift.library.wiener_filter_curvature.WienerFilterCurvature(R,N,Sh,inverter)\n",
"\n",
"for i in range(20):\n",
" print i\n",
" sc.add(HT(curv.generate_posterior_sample()))\n",
"\n",
"m_var = sc.var\n",
"m_mean, m_var = ift.probe_with_posterior_samples(curv, m, HT, 20)\n",
"\n",
"# Get data\n",
"s_power = ift.power_analyze(sh)\n",
......
......@@ -108,5 +108,5 @@ if __name__ == "__main__":
ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict)
ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict)
# sampling the uncertainty map
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
mean, variance = ift.probe_with_posterior_samples(wiener_curvature, m_k, ht, 10)
ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)
......@@ -72,7 +72,7 @@ if __name__ == "__main__":
sample_mean = ift.Field.zeros(signal_space)
n_samples = 10
for i in range(n_samples):
sample = ht(wiener_curvature.generate_posterior_sample()) + m
sample = ht(wiener_curvature.draw_sample()) + m
sample_variance += sample**2
sample_mean += sample
variance = sample_variance/n_samples - (sample_mean/n_samples)**2
......
......@@ -62,9 +62,12 @@ if __name__ == "__main__":
D = ift.InversionEnabler(D, inverter)
m = D(j)
# Plotting
d_field = ift.Field(s_space, val=d.val)
zmax = max(ht(sh).max(), d_field.max(), ht(m).max())
zmin = min(ht(sh).min(), d_field.min(), ht(m).min())
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
"colormap": "Planck-like", "zmax": zmax, "zmin": zmin}
ift.plot(ht(sh), name="mock_signal.png", **plotdict)
ift.plot(ift.Field(s_space, val=d.val),
name="data.png", **plotdict)
ift.plot(ht(m), name="map.png", **plotdict)
ift.plot(d_field, name="data.png", **plotdict)
ift.plot(ht(m), name="reconstruction.png", **plotdict)
......@@ -65,7 +65,6 @@ if __name__ == "__main__":
domain=data_domain, random_type='normal',
std=noise_amplitude, mean=0)
data = noiseless_data + noise
# Wiener filter
j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(
......@@ -79,10 +78,10 @@ if __name__ == "__main__":
sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
ift.plot(ift.Field(sspace2, ht(mock_signal).val)/nu.K, name="mock_signal.png")
#data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)
#data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))
zmax = max(m_s.max(), ht(mock_signal).max())
zmin = min(m_s.min(), ht(mock_signal).min())
plotdict = {"zmax": zmax/nu.K, "zmin": zmin/nu.K}
ift.plot(ift.Field(sspace2, ht(mock_signal).val)/nu.K, name="mock_signal.png", **plotdict)
ift.plot(ift.Field(sspace2, val=data.val), name="data.png")
print "msig",np.min(ht(mock_signal).val)/nu.K, np.max(ht(mock_signal).val)/nu.K
print "map",np.min(m_s.val)/nu.K, np.max(m_s.val)/nu.K
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
ift.plot(ift.Field(sspace2, m_s.val)/nu.K, name="reconstruction.png", **plotdict)
......@@ -6,33 +6,38 @@ np.random.seed(42)
if __name__ == "__main__":
# Set up position space
s_space = ift.RGSpace([128, 128])
# s_space = ift.HPSpace(32)
# s_space = ift.RGSpace([128, 128])
s_space = ift.HPSpace(32)
# Define associated harmonic space and harmonic transformation
h_space = s_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(h_space, s_space)
# Setting up power space
# Set up power space
p_space = ift.PowerSpace(h_space)
# Choosing the prior correlation structure and defining
# correlation operator
# Choose prior correlation structure and define correlation operator
p_spec = (lambda k: (42/(k+1)**3))
S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space
# Draw sample sh from the prior distribution in harmonic space
sp = ift.PS_field(p_space, p_spec)
sh = ift.power_synthesize(sp, real_signal=True)
# Choosing the measurement instrument
# Choose measurement instrument
diag = np.ones(s_space.shape)
diag[20:80, 20:80] = 0
if len(s_space.shape) == 1:
diag[3000:7000] = 0
elif len(s_space.shape) == 2:
diag[20:80, 20:80] = 0
else:
raise NotImplementedError
diag = ift.Field(s_space, ift.dobj.from_global_data(diag))
Instrument = ift.DiagonalOperator(diag)
# Adding a harmonic transformation to the instrument
# Add harmonic transformation to the instrument
R = Instrument*ht
noiseless_data = R(sh)
signal_to_noise = 1.
......@@ -43,31 +48,36 @@ if __name__ == "__main__":
std=noise_amplitude,
mean=0)
# Creating the mock data
# Create mock data
d = noiseless_data + n
j = R.adjoint_times(N.inverse_times(d))
# Choosing the minimization strategy
# Choose minimization strategy
ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl)
controller = ift.GradientNormController(name="min", tol_abs_gradnorm=0.1)
minimizer = ift.RelaxedNewton(controller=controller)
m0 = ift.Field.zeros(h_space)
# Initializing the Wiener Filter energy
# Initialize Wiener filter energy
energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
inverter=inverter)
energy, convergence = minimizer(energy)
m = energy.position
D = energy.curvature
ift.plot(ht(sh), name="signal.png", colormap="Planck-like")
ift.plot(ht(m), name="m.png", colormap="Planck-like")
# sampling the uncertainty map
# Generate plots
zmax = max(ht(sh).max(), ht(m).max())
zmin = min(ht(sh).min(), ht(m).min())
plotdict = {"zmax": zmax, "zmin": zmin, "colormap": "Planck-like"}
ift.plot(ht(sh), name="mock_signal.png", **plotdict)
ift.plot(ht(m), name="reconstruction.png", **plotdict)
# Sample uncertainty map
sample_variance = ift.Field.zeros(s_space)
sample_mean = ift.Field.zeros(s_space)
mean, variance = ift.probe_with_posterior_samples(D, ht, 50)
ift.plot(variance, name="variance.png", colormap="Planck-like")
ift.plot(mean, name="mean.png", colormap="Planck-like")
mean, variance = ift.probe_with_posterior_samples(D, sh, ht, 50)
ift.plot(variance, name="posterior_variance.png", colormap="Planck-like")
ift.plot(mean, name="posterior_mean.png", **plotdict)
......@@ -91,7 +91,7 @@ class CriticalPowerEnergy(Energy):
if self.D is not None:
w = Field.zeros(self.position.domain, dtype=self.m.dtype)
for i in range(self.samples):
sample = self.D.generate_posterior_sample() + self.m
sample = self.D.draw_sample() + self.m
w += P(abs(sample)**2)
w *= 1./self.samples
......
......@@ -46,7 +46,7 @@ class NoiseEnergy(Energy):
if samples is None or samples == 0:
xi_sample_list = [xi]
else:
xi_sample_list = [D.generate_posterior_sample() + xi
xi_sample_list = [D.draw_sample() + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
......
......@@ -70,7 +70,7 @@ class NonlinearPowerEnergy(Energy):
if samples is None or samples == 0:
xi_sample_list = [xi]
else:
xi_sample_list = [D.generate_posterior_sample() + xi
xi_sample_list = [D.draw_sample() + xi
for _ in range(samples)]
self.xi_sample_list = xi_sample_list
self.inverter = inverter
......
......@@ -59,23 +59,22 @@ class WienerFilterCurvature(EndomorphicOperator):
def apply(self, x, mode):
return self._op.apply(x, mode)
def generate_posterior_sample(self):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
def draw_sample(self):
""" Generates a sample from a Gaussian distribution with
covariance given by the operator.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
right correlation.
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
Returns the a sample from the Gaussian of given covariance.
"""
mock_signal = self.S.generate_posterior_sample()
mock_noise = self.N.generate_posterior_sample()
mock_signal = self.S.draw_sample()
mock_noise = self.N.draw_sample()
mock_data = self.R(mock_signal) + mock_noise
......
......@@ -141,19 +141,18 @@ class DiagonalOperator(EndomorphicOperator):
return DiagonalOperator(self._diagonal.conjugate(), self._domain,
self._spaces)
def generate_posterior_sample(self):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
def draw_sample(self):
""" Generates a sample from a Gaussian distribution with
covariance given by the operator.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
right correlation.
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
Returns the a sample from the Gaussian of given covariance.
"""
if self._spaces is not None:
......
......@@ -89,19 +89,18 @@ class ScalingOperator(EndomorphicOperator):
return (self.TIMES | self.ADJOINT_TIMES |
self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
def generate_posterior_sample(self):
""" Generates a posterior sample from a Gaussian distribution with
given mean and covariance.
def draw_sample(self):
""" Generates a sample from a Gaussian distribution with
covariance given by the operator.
This method generates samples by setting up the observation and
reconstruction of a mock signal in order to obtain residuals of the
right correlation which are added to the given mean.
right correlation.
Returns
-------
sample : Field
Returns the a sample from the Gaussian of given mean and
covariance.
Returns the a sample from the Gaussian of given covariance.
"""
return Field.from_random(random_type="normal",
......
......@@ -47,9 +47,12 @@ class StatCalculator(object):
return self._M2 * (1./(self._count-1))
def probe_with_posterior_samples(op, post_op, nprobes):
def probe_with_posterior_samples(op, m, post_op, nprobes):
sc = StatCalculator()
for i in range(nprobes):
sample = post_op(op.generate_posterior_sample())
sample = post_op(op.draw_sample() + m)
sc.add(sample)
if nprobes == 1:
return sc.mean, None
return sc.mean, sc.var
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