Commit 542bf762 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

tweaks

parent ac0c8845
Pipeline #30505 failed with stages
in 42 seconds
......@@ -38,14 +38,17 @@ N_iter = 300
IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=N_iter)
inverter = ift.ConjugateGradient(IC)
sampling_inverter = ift.ConjugateGradient(IC)
D_inv_1 = ift.SamplingEnabler(ift.SandwichOperator.make(R_p, N.inverse), S.inverse, inverter, sampling_inverter)
D_inv_1 = ift.SamplingEnabler(ift.SandwichOperator.make(R_p, N.inverse),
S.inverse, sampling_inverter)
D_inv_1 = ift.InversionEnabler(D_inv_1, inverter)
D_inv_2 = ift.SamplingEnabler2(D_inv, inverter, sampling_inverter)
D_inv_2 = ift.SamplingEnabler2(D_inv, sampling_inverter)
samps_1 = [D_inv_1.draw_sample(from_inverse=True) for i in range(N_samps)] #GOOD
samps_2 = [D_inv_2.draw_sample(from_inverse=True) for i in range(N_samps)] #BAD
# GOOD
samps_1 = [D_inv_1.draw_sample(from_inverse=True) for i in range(N_samps)]
# BAD
samps_2 = [D_inv_2.draw_sample(from_inverse=True) for i in range(N_samps)]
m = D_inv_1.inverse_times(j)
m_x = sky(m)
......@@ -63,7 +66,8 @@ for i in range(N_samps):
else:
plt.plot(sky(samps_2[i]).to_global_data()[290:360], color='b',
**pltdict)
plt.plot(sky(samps_1[i]).to_global_data()[290:360], color='r', **pltdict)
plt.plot(sky(samps_1[i]).to_global_data()[290:360], color='r',
**pltdict)
plt.plot((s_x - m_x).to_global_data()[290:360], color='k',
label='signal - mean')
plt.title('Comparison of conservative vs default samples near area of low noise')
......
......@@ -18,4 +18,5 @@ __all__ = ["LinearOperator", "EndomorphicOperator", "ScalingOperator",
"DiagonalOperator", "HarmonicTransformOperator", "FFTOperator",
"FFTSmoothingOperator", "GeometryRemover",
"LaplaceOperator", "SmoothnessOperator", "PowerDistributor",
"InversionEnabler", "SandwichOperator",'SamplingEnabler', 'SamplingEnabler2']
"InversionEnabler", "SandwichOperator", "SamplingEnabler",
"SamplingEnabler2"]
......@@ -54,10 +54,6 @@ class InversionEnabler(EndomorphicOperator):
def domain(self):
return self._op.domain
@property
def target(self):
return self._op.target
@property
def capability(self):
return self._addInverse[self._op.capability]
......
......@@ -19,11 +19,11 @@
from ..minimization.quadratic_energy import QuadraticEnergy
from ..minimization.iteration_controller import IterationController
from ..logger import logger
from .inversion_enabler import InversionEnabler
from .endomorphic_operator import EndomorphicOperator
import numpy as np
class SamplingEnabler(InversionEnabler):
class SamplingEnabler(EndomorphicOperator):
"""Class which augments the capability of another operator object via
numerical inversion.
......@@ -44,20 +44,33 @@ class SamplingEnabler(InversionEnabler):
convergence.
"""
def __init__(self, likelihood, prior, application_inverter, sampling_inverter, approximation=None):
def __init__(self, likelihood, prior, sampling_inverter,
approximation=None):
self._op = likelihood + prior
super(SamplingEnabler, self).__init__(self._op, application_inverter, approximation=approximation )
self.likelihood = likelihood
self.prior = prior
self.sampling_inverter = sampling_inverter
super(SamplingEnabler, self).__init__()
self._likelihood = likelihood
self._prior = prior
self._sampling_inverter = sampling_inverter
def draw_sample(self, from_inverse=False, dtype=np.float64):
try:
return self._op.draw_sample(from_inverse, dtype)
except NotImplementedError:
s = self.prior.draw_sample()
sp = self.prior.inverse_times(s)
nj = self.likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj, _grad=self.likelihood(s) - nj)
energy, convergence = self.sampling_inverter(energy)
s = self._prior.draw_sample()
sp = self._prior.inverse_times(s)
nj = self._likelihood.draw_sample()
energy = QuadraticEnergy(s, self._op, sp + nj,
_grad=self._likelihood(s) - nj)
energy, convergence = self._sampling_inverter(energy)
return energy.position
@property
def domain(self):
return self._op.domain
@property
def capability(self):
return self._op.capability
def apply(self, x, mode):
return self._op.apply(x, mode)
......@@ -19,10 +19,12 @@
from ..minimization.quadratic_energy import QuadraticEnergy
from ..minimization.iteration_controller import IterationController
from ..logger import logger
from .endomorphic_operator import EndomorphicOperator
from .inversion_enabler import InversionEnabler
import numpy as np
class SamplingEnabler2(InversionEnabler):
class SamplingEnabler2(EndomorphicOperator):
"""Class which augments the capability of another operator object via
numerical inversion.
......@@ -43,13 +45,27 @@ class SamplingEnabler2(InversionEnabler):
convergence.
"""
def __init__(self, op, application_inverter, sampling_inverter, approximation=None):
super(SamplingEnabler2, self).__init__(op, application_inverter, approximation=approximation)
self.sampling_op = InversionEnabler(op, sampling_inverter,approximation=approximation)
def __init__(self, op, sampling_inverter, approximation=None):
super(SamplingEnabler2, self).__init__()
self._op = op
self._sampling_op = InversionEnabler(op, sampling_inverter,
approximation=approximation)
def draw_sample(self, from_inverse=False, dtype=np.float64):
try:
return self._op.draw_sample(from_inverse, dtype)
except NotImplementedError:
samp = self._op.draw_sample(not from_inverse, dtype)
return self.sampling_op.inverse_times(samp) if from_inverse else self(samp)
return self._sampling_op.inverse_times(samp) \
if from_inverse else self(samp)
@property
def domain(self):
return self._op.domain
@property
def capability(self):
return self._op.capability
def apply(self, x, mode):
return self._op.apply(x, mode)
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