Commit e7d4b853 authored by Martin Reinecke's avatar Martin Reinecke

avoid assert statements; simplify interface to GradientNormController

parent a5859e54
Pipeline #24186 passed with stage
in 4 minutes and 32 seconds
...@@ -53,7 +53,7 @@ if __name__ == "__main__": ...@@ -53,7 +53,7 @@ if __name__ == "__main__":
d_data = d.val d_data = d.val
ift.plot(d, name="data.png") ift.plot(d, name="data.png")
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100, IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
minimizer = ift.RelaxedNewton(IC1) minimizer = ift.RelaxedNewton(IC1)
......
...@@ -47,9 +47,8 @@ if __name__ == "__main__": ...@@ -47,9 +47,8 @@ if __name__ == "__main__":
# Wiener filter # Wiener filter
m0 = ift.Field.zeros(harmonic_space) m0 = ift.Field.zeros(harmonic_space)
ctrl = ift.GradientNormController(verbose=False, tol_abs_gradnorm=0.0001) ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001)
ctrl2 = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1, ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer")
name="outer")
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic,
N, S, inverter=inverter) N, S, inverter=inverter)
......
...@@ -72,13 +72,12 @@ if __name__ == "__main__": ...@@ -72,13 +72,12 @@ if __name__ == "__main__":
t0 = ift.Field(p_space, val=-4.) t0 = ift.Field(p_space, val=-4.)
power0 = Projection.adjoint_times(ift.exp(0.5 * t0)) power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100, IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
tol_abs_gradnorm=1e-3) tol_abs_gradnorm=1e-3)
LS = ift.LineSearchStrongWolfe(c2=0.02) LS = ift.LineSearchStrongWolfe(c2=0.02)
minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
ICI = ift.GradientNormController(verbose=False, name="ICI", ICI = ift.GradientNormController(iteration_limit=500,
iteration_limit=500,
tol_abs_gradnorm=1e-3) tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI) inverter = ift.ConjugateGradient(controller=ICI)
......
...@@ -86,7 +86,7 @@ if __name__ == "__main__": ...@@ -86,7 +86,7 @@ if __name__ == "__main__":
# Wiener filter # Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R_harmonic.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter) S=S, N=N, R=R_harmonic, inverter=inverter)
......
...@@ -19,7 +19,7 @@ if __name__ == "__main__": ...@@ -19,7 +19,7 @@ if __name__ == "__main__":
signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
harmonic_space = signal_space.get_default_codomain() harmonic_space = signal_space.get_default_codomain()
fft = ift.FFTOperator(harmonic_space, target=signal_space) fft = ift.HarmonicTransformOperator(harmonic_space, target=signal_space)
power_space = ift.PowerSpace( power_space = ift.PowerSpace(
harmonic_space, binbounds=ift.PowerSpace.useful_binbounds( harmonic_space, binbounds=ift.PowerSpace.useful_binbounds(
harmonic_space, logarithmic=True)) harmonic_space, logarithmic=True))
...@@ -28,40 +28,47 @@ if __name__ == "__main__": ...@@ -28,40 +28,47 @@ if __name__ == "__main__":
S = ift.create_power_operator(harmonic_space, S = ift.create_power_operator(harmonic_space,
power_spectrum=power_spectrum) power_spectrum=power_spectrum)
mock_power = ift.PS_field(power_space, power_spectrum) mock_power = ift.PS_field(power_space, power_spectrum)
mock_signal = fft(ift.power_synthesize(mock_power, real_signal=True)) mock_signal = ift.power_synthesize(mock_power, real_signal=True)
# Setting up an exemplary response # Setting up an exemplary response
mask = np.ones(signal_space.shape) mask = np.ones(signal_space.shape)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
mask[N10*5:N10*9, N10*5:N10*9] = 0. mask[N10*5:N10*9, N10*5:N10*9] = 0.
mask = ift.Field(signal_space, ift.dobj.from_global_data(mask)) mask = ift.Field(signal_space, ift.dobj.from_global_data(mask))
R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), R = ift.GeometryRemover(signal_space)
exposure=(mask,)) R = R*ift.DiagonalOperator(mask)
R = R*fft
R = R * ift.create_harmonic_smoothing_operator((harmonic_space,),0,response_sigma)
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = R * fft
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) ndiag = 1e-8*ift.Field.full(data_domain, fft(mock_signal).var()/signal_to_noise)
N = ift.DiagonalOperator(ndiag.weight(1)) N = ift.DiagonalOperator(ndiag)
noise = ift.Field.from_random( noise = ift.Field.from_random(
domain=data_domain, random_type='normal', domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(mock_signal) + noise data = R(mock_signal) + noise
# Wiener filter # Wiener filter
j = R_harmonic.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-6)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R_harmonic, inverter=inverter) S=S, N=N, R=R, inverter=inverter)
m_k = wiener_curvature.inverse_times(j) m_k = wiener_curvature.inverse_times(j)
m = fft(m_k) m = fft(m_k)
plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
"colormap": "Planck-like"}
ift.plot(mock_signal, name="mock_signal.png", **plotdict)
ift.plot(ift.Field(signal_space, val=data.val),
name="data.png", **plotdict)
ift.plot(m, name="map.png", **plotdict)
# Probing the uncertainty # Probing the uncertainty
class Proby(ift.DiagonalProberMixin, ift.Prober): class Proby(ift.DiagonalProberMixin, ift.Prober):
pass pass
proby = Proby(signal_space, probe_count=1, ncpu=1) proby = Proby(harmonic_space, probe_count=1, ncpu=1)
proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) proby(lambda z: wiener_curvature.inverse_times(z))
sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03)
variance = ift.sqrt(sm(proby.diagonal.weight(-1))) variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
......
...@@ -54,7 +54,7 @@ if __name__ == "__main__": ...@@ -54,7 +54,7 @@ if __name__ == "__main__":
# Wiener filter # Wiener filter
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
IC = ift.GradientNormController(verbose=True, iteration_limit=500, IC = ift.GradientNormController(name="inverter", iteration_limit=500,
tol_abs_gradnorm=0.1) tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=IC) inverter = ift.ConjugateGradient(controller=IC)
S_inv = fft.adjoint*Sh.inverse*fft S_inv = fft.adjoint*Sh.inverse*fft
......
...@@ -69,7 +69,7 @@ if __name__ == "__main__": ...@@ -69,7 +69,7 @@ if __name__ == "__main__":
j = R.adjoint_times(N.inverse_times(data)) j = R.adjoint_times(N.inverse_times(data))
ctrl = ift.GradientNormController( ctrl = ift.GradientNormController(
verbose=True, tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality))) name="inverter", tol_abs_gradnorm=1e-5/(nu.K*(nu.m**dimensionality)))
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
wiener_curvature = ift.library.WienerFilterCurvature( wiener_curvature = ift.library.WienerFilterCurvature(
S=S, N=N, R=R, inverter=inverter) S=S, N=N, R=R, inverter=inverter)
......
...@@ -50,9 +50,9 @@ if __name__ == "__main__": ...@@ -50,9 +50,9 @@ if __name__ == "__main__":
# Choosing the minimization strategy # Choosing the minimization strategy
ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl) inverter = ift.ConjugateGradient(controller=ctrl)
controller = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) controller = ift.GradientNormController(name="min", tol_abs_gradnorm=0.1)
minimizer = ift.RelaxedNewton(controller=controller) minimizer = ift.RelaxedNewton(controller=controller)
m0 = ift.Field.zeros(h_space) m0 = ift.Field.zeros(h_space)
# Initializing the Wiener Filter energy # Initializing the Wiener Filter energy
......
...@@ -22,15 +22,13 @@ from .. import dobj ...@@ -22,15 +22,13 @@ from .. import dobj
class GradientNormController(IterationController): class GradientNormController(IterationController):
def __init__(self, tol_abs_gradnorm=None, tol_rel_gradnorm=None, def __init__(self, tol_abs_gradnorm=None, tol_rel_gradnorm=None,
convergence_level=1, iteration_limit=None, convergence_level=1, iteration_limit=None, name=None):
name=None, verbose=None):
super(GradientNormController, self).__init__() super(GradientNormController, self).__init__()
self._tol_abs_gradnorm = tol_abs_gradnorm self._tol_abs_gradnorm = tol_abs_gradnorm
self._tol_rel_gradnorm = tol_rel_gradnorm self._tol_rel_gradnorm = tol_rel_gradnorm
self._convergence_level = convergence_level self._convergence_level = convergence_level
self._iteration_limit = iteration_limit self._iteration_limit = iteration_limit
self._name = name self._name = name
self._verbose = verbose
def start(self, energy): def start(self, energy):
self._itcount = -1 self._itcount = -1
...@@ -56,10 +54,8 @@ class GradientNormController(IterationController): ...@@ -56,10 +54,8 @@ class GradientNormController(IterationController):
self._ccount = max(0, self._ccount-1) self._ccount = max(0, self._ccount-1)
# report # report
if self._verbose: if self._name is not None:
msg = "" msg = self._name+":"
if self._name is not None:
msg += self._name+":"
msg += " Iteration #" + str(self._itcount) msg += " Iteration #" + str(self._itcount)
msg += " energy=" + str(energy.value) msg += " energy=" + str(energy.value)
msg += " gradnorm=" + str(energy.gradient_norm) msg += " gradnorm=" + str(energy.gradient_norm)
......
...@@ -209,8 +209,10 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -209,8 +209,10 @@ class LineSearchStrongWolfe(LineSearch):
alpha_recent = None alpha_recent = None
phi_recent = None phi_recent = None
assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0 if phi_lo > phi_0 + self.c1*alpha_lo*phiprime_0:
assert phiprime_lo*(alpha_hi-alpha_lo) < 0. raise ValueError("inconsistent data")
if phiprime_lo*(alpha_hi-alpha_lo) >= 0.:
raise ValueError("inconsistent data")
for i in range(self.max_zoom_iterations): for i in range(self.max_zoom_iterations):
# assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0 # assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
# assert phiprime_lo*(alpha_hi-alpha_lo)<0. # assert phiprime_lo*(alpha_hi-alpha_lo)<0.
......
...@@ -149,8 +149,9 @@ class HarmonicTransformOperator(LinearOperator): ...@@ -149,8 +149,9 @@ class HarmonicTransformOperator(LinearOperator):
def _slice_p2h(self, inp): def _slice_p2h(self, inp):
rr = self.sjob.alm2map_adjoint(inp) rr = self.sjob.alm2map_adjoint(inp)
assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \ if len(rr) != ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax) (self.mmax+1)*(self.lmax-self.mmax):
raise ValueError("array length mismatch")
res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype) res = np.empty(2*len(rr)-self.lmax-1, dtype=rr[0].real.dtype)
res[0:self.lmax+1] = rr[0:self.lmax+1].real res[0:self.lmax+1] = rr[0:self.lmax+1].real
res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real
...@@ -159,8 +160,9 @@ class HarmonicTransformOperator(LinearOperator): ...@@ -159,8 +160,9 @@ class HarmonicTransformOperator(LinearOperator):
def _slice_h2p(self, inp): def _slice_h2p(self, inp):
res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype) res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype)
assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \ if len(res) != ((self.mmax+1)*(self.mmax+2))//2 + \
(self.mmax+1)*(self.lmax-self.mmax) (self.mmax+1)*(self.lmax-self.mmax):
raise ValueError("array length mismatch")
res[0:self.lmax+1] = inp[0:self.lmax+1] res[0:self.lmax+1] = inp[0:self.lmax+1]
res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] + res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] +
1j*inp[self.lmax+2::2]) 1j*inp[self.lmax+2::2])
......
...@@ -62,7 +62,8 @@ class PowerSpace(Space): ...@@ -62,7 +62,8 @@ class PowerSpace(Space):
units of the harmonic partner space. units of the harmonic partner space.
""" """
nbin = int(nbin) nbin = int(nbin)
assert nbin >= 3, "nbin must be at least 3" if nbin < 3:
raise ValueError("nbin must be at least 3")
return np.linspace(float(first_bound), float(last_bound), nbin-1) return np.linspace(float(first_bound), float(last_bound), nbin-1)
@staticmethod @staticmethod
...@@ -81,7 +82,8 @@ class PowerSpace(Space): ...@@ -81,7 +82,8 @@ class PowerSpace(Space):
units of the harmonic partner space. units of the harmonic partner space.
""" """
nbin = int(nbin) nbin = int(nbin)
assert nbin >= 3, "nbin must be at least 3" if nbin < 3:
raise ValueError("nbin must be at least 3")
return np.logspace(np.log(float(first_bound)), return np.logspace(np.log(float(first_bound)),
np.log(float(last_bound)), np.log(float(last_bound)),
nbin-1, base=np.e) nbin-1, base=np.e)
...@@ -107,7 +109,8 @@ class PowerSpace(Space): ...@@ -107,7 +109,8 @@ class PowerSpace(Space):
nbin_max = int((dists[-1]-dists[0])/binsz_min)+2 nbin_max = int((dists[-1]-dists[0])/binsz_min)+2
if nbin is None: if nbin is None:
nbin = nbin_max nbin = nbin_max
assert nbin >= 3, "nbin must be at least 3" if nbin < 3:
raise ValueError("nbin must be at least 3")
if nbin > nbin_max: if nbin > nbin_max:
raise ValueError("nbin is too large") raise ValueError("nbin is too large")
if logarithmic: if logarithmic:
...@@ -141,16 +144,19 @@ class PowerSpace(Space): ...@@ -141,16 +144,19 @@ class PowerSpace(Space):
tbb = binbounds tbb = binbounds
locdat = np.searchsorted(tbb, dobj.local_data(k_length_array.val)) locdat = np.searchsorted(tbb, dobj.local_data(k_length_array.val))
temp_pindex = dobj.from_local_data( temp_pindex = dobj.from_local_data(
k_length_array.val.shape, locdat, dobj.distaxis(k_length_array.val)) k_length_array.val.shape, locdat,
dobj.distaxis(k_length_array.val))
nbin = len(tbb)+1 nbin = len(tbb)+1
temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(), temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(),
minlength=nbin) minlength=nbin)
temp_rho = dobj.np_allreduce_sum(temp_rho) temp_rho = dobj.np_allreduce_sum(temp_rho)
assert not (temp_rho == 0).any(), "empty bins detected" if (temp_rho == 0).any():
raise ValueError("empty bins detected")
# The explicit conversion to float64 is necessary because bincount # The explicit conversion to float64 is necessary because bincount
# sometimes returns its result as an integer array, even when # sometimes returns its result as an integer array, even when
# floating-point weights are present ... # floating-point weights are present ...
temp_k_lengths = np.bincount(dobj.local_data(temp_pindex).ravel(), temp_k_lengths = np.bincount(
dobj.local_data(temp_pindex).ravel(),
weights=dobj.local_data(k_length_array.val).ravel(), weights=dobj.local_data(k_length_array.val).ravel(),
minlength=nbin).astype(np.float64, copy=False) minlength=nbin).astype(np.float64, copy=False)
temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho
......
...@@ -41,7 +41,7 @@ class Test_Minimizers(unittest.TestCase): ...@@ -41,7 +41,7 @@ class Test_Minimizers(unittest.TestCase):
covariance = ift.DiagonalOperator(covariance_diagonal) covariance = ift.DiagonalOperator(covariance_diagonal)
required_result = ift.Field.ones(space, dtype=np.float64) required_result = ift.Field.ones(space, dtype=np.float64)
IC = ift.GradientNormController(verbose=True,tol_abs_gradnorm=1e-5, iteration_limit=1000) IC = ift.GradientNormController(tol_abs_gradnorm=1e-5, iteration_limit=1000)
try: try:
minimizer = minimizer_class(controller=IC) minimizer = minimizer_class(controller=IC)
energy = ift.QuadraticEnergy(A=covariance, b=required_result, energy = ift.QuadraticEnergy(A=covariance, b=required_result,
......
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