From e7d4b8530bf185c7aa1324b5743bc6712ac4b6e8 Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Tue, 30 Jan 2018 10:25:35 +0100 Subject: [PATCH] avoid assert statements; simplify interface to GradientNormController --- demos/critical_filtering.py | 2 +- demos/log_normal_wiener_filter.py | 5 ++- demos/nonlinear_critical_filter.py | 5 ++- demos/paper_demos/cartesian_wiener_filter.py | 2 +- demos/paper_demos/wiener_filter.py | 31 ++++++++++++------- demos/wiener_filter_easy.py | 2 +- demos/wiener_filter_via_curvature.py | 2 +- demos/wiener_filter_via_hamiltonian.py | 4 +-- .../minimization/gradient_norm_controller.py | 10 ++---- .../minimization/line_search_strong_wolfe.py | 6 ++-- .../operators/harmonic_transform_operator.py | 10 +++--- nifty4/spaces/power_space.py | 18 +++++++---- test/test_minimization/test_minimizers.py | 2 +- 13 files changed, 55 insertions(+), 44 deletions(-) diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index 3fef05a5..7866f04d 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -53,7 +53,7 @@ if __name__ == "__main__": d_data = d.val 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) minimizer = ift.RelaxedNewton(IC1) diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py index 87563deb..5ab51353 100644 --- a/demos/log_normal_wiener_filter.py +++ b/demos/log_normal_wiener_filter.py @@ -47,9 +47,8 @@ if __name__ == "__main__": # Wiener filter m0 = ift.Field.zeros(harmonic_space) - ctrl = ift.GradientNormController(verbose=False, tol_abs_gradnorm=0.0001) - ctrl2 = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1, - name="outer") + ctrl = ift.GradientNormController(tol_abs_gradnorm=0.0001) + ctrl2 = ift.GradientNormController(tol_abs_gradnorm=0.1, name="outer") inverter = ift.ConjugateGradient(controller=ctrl) energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter) diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py index 7feef742..9a67ca09 100644 --- a/demos/nonlinear_critical_filter.py +++ b/demos/nonlinear_critical_filter.py @@ -72,13 +72,12 @@ if __name__ == "__main__": t0 = ift.Field(p_space, val=-4.) 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) LS = ift.LineSearchStrongWolfe(c2=0.02) minimizer = ift.RelaxedNewton(IC1, line_searcher=LS) - ICI = ift.GradientNormController(verbose=False, name="ICI", - iteration_limit=500, + ICI = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3) inverter = ift.ConjugateGradient(controller=ICI) diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py index 4fd696d0..9b466542 100644 --- a/demos/paper_demos/cartesian_wiener_filter.py +++ b/demos/paper_demos/cartesian_wiener_filter.py @@ -86,7 +86,7 @@ if __name__ == "__main__": # Wiener filter 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) wiener_curvature = ift.library.WienerFilterCurvature( S=S, N=N, R=R_harmonic, inverter=inverter) diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py index ddf5a63a..b0d682fc 100644 --- a/demos/paper_demos/wiener_filter.py +++ b/demos/paper_demos/wiener_filter.py @@ -19,7 +19,7 @@ if __name__ == "__main__": signal_space = ift.RGSpace([N_pixels, N_pixels], distances=L/N_pixels) 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( harmonic_space, binbounds=ift.PowerSpace.useful_binbounds( harmonic_space, logarithmic=True)) @@ -28,40 +28,47 @@ if __name__ == "__main__": S = ift.create_power_operator(harmonic_space, power_spectrum=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 mask = np.ones(signal_space.shape) N10 = int(N_pixels/10) mask[N10*5:N10*9, N10*5:N10*9] = 0. mask = ift.Field(signal_space, ift.dobj.from_global_data(mask)) - R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), - exposure=(mask,)) + R = ift.GeometryRemover(signal_space) + 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] - R_harmonic = R * fft # Setting up the noise covariance and drawing a random noise realization - ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise) - N = ift.DiagonalOperator(ndiag.weight(1)) + ndiag = 1e-8*ift.Field.full(data_domain, fft(mock_signal).var()/signal_to_noise) + N = ift.DiagonalOperator(ndiag) noise = ift.Field.from_random( domain=data_domain, random_type='normal', std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) data = R(mock_signal) + noise # Wiener filter - j = R_harmonic.adjoint_times(N.inverse_times(data)) - ctrl = ift.GradientNormController(verbose=True, tol_abs_gradnorm=0.1) + j = R.adjoint_times(N.inverse_times(data)) + ctrl = ift.GradientNormController(name="inverter", tol_abs_gradnorm=1e-6) inverter = ift.ConjugateGradient(controller=ctrl) 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 = 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 class Proby(ift.DiagonalProberMixin, ift.Prober): pass - proby = Proby(signal_space, probe_count=1, ncpu=1) - proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) + proby = Proby(harmonic_space, probe_count=1, ncpu=1) + proby(lambda z: wiener_curvature.inverse_times(z)) sm = ift.FFTSmoothingOperator(signal_space, sigma=0.03) variance = ift.sqrt(sm(proby.diagonal.weight(-1))) diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py index 5884b145..09f94632 100644 --- a/demos/wiener_filter_easy.py +++ b/demos/wiener_filter_easy.py @@ -54,7 +54,7 @@ if __name__ == "__main__": # Wiener filter 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) inverter = ift.ConjugateGradient(controller=IC) S_inv = fft.adjoint*Sh.inverse*fft diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py index 16ee1acb..0e3309b7 100644 --- a/demos/wiener_filter_via_curvature.py +++ b/demos/wiener_filter_via_curvature.py @@ -69,7 +69,7 @@ if __name__ == "__main__": j = R.adjoint_times(N.inverse_times(data)) 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) wiener_curvature = ift.library.WienerFilterCurvature( S=S, N=N, R=R, inverter=inverter) diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py index eeeff781..74f5ae90 100644 --- a/demos/wiener_filter_via_hamiltonian.py +++ b/demos/wiener_filter_via_hamiltonian.py @@ -50,9 +50,9 @@ if __name__ == "__main__": # 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) - 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) m0 = ift.Field.zeros(h_space) # Initializing the Wiener Filter energy diff --git a/nifty4/minimization/gradient_norm_controller.py b/nifty4/minimization/gradient_norm_controller.py index 2c9c9a4f..c06a9988 100644 --- a/nifty4/minimization/gradient_norm_controller.py +++ b/nifty4/minimization/gradient_norm_controller.py @@ -22,15 +22,13 @@ from .. import dobj class GradientNormController(IterationController): def __init__(self, tol_abs_gradnorm=None, tol_rel_gradnorm=None, - convergence_level=1, iteration_limit=None, - name=None, verbose=None): + convergence_level=1, iteration_limit=None, name=None): super(GradientNormController, self).__init__() self._tol_abs_gradnorm = tol_abs_gradnorm self._tol_rel_gradnorm = tol_rel_gradnorm self._convergence_level = convergence_level self._iteration_limit = iteration_limit self._name = name - self._verbose = verbose def start(self, energy): self._itcount = -1 @@ -56,10 +54,8 @@ class GradientNormController(IterationController): self._ccount = max(0, self._ccount-1) # report - if self._verbose: - msg = "" - if self._name is not None: - msg += self._name+":" + if self._name is not None: + msg = self._name+":" msg += " Iteration #" + str(self._itcount) msg += " energy=" + str(energy.value) msg += " gradnorm=" + str(energy.gradient_norm) diff --git a/nifty4/minimization/line_search_strong_wolfe.py b/nifty4/minimization/line_search_strong_wolfe.py index 1a139892..7987a627 100644 --- a/nifty4/minimization/line_search_strong_wolfe.py +++ b/nifty4/minimization/line_search_strong_wolfe.py @@ -209,8 +209,10 @@ class LineSearchStrongWolfe(LineSearch): alpha_recent = None phi_recent = None - assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0 - assert phiprime_lo*(alpha_hi-alpha_lo) < 0. + if phi_lo > phi_0 + self.c1*alpha_lo*phiprime_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): # assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0 # assert phiprime_lo*(alpha_hi-alpha_lo)<0. diff --git a/nifty4/operators/harmonic_transform_operator.py b/nifty4/operators/harmonic_transform_operator.py index c7e368cc..d3d8cf29 100644 --- a/nifty4/operators/harmonic_transform_operator.py +++ b/nifty4/operators/harmonic_transform_operator.py @@ -149,8 +149,9 @@ class HarmonicTransformOperator(LinearOperator): def _slice_p2h(self, inp): rr = self.sjob.alm2map_adjoint(inp) - assert len(rr) == ((self.mmax+1)*(self.mmax+2))//2 + \ - (self.mmax+1)*(self.lmax-self.mmax) + if len(rr) != ((self.mmax+1)*(self.mmax+2))//2 + \ + (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[0:self.lmax+1] = rr[0:self.lmax+1].real res[self.lmax+1::2] = np.sqrt(2)*rr[self.lmax+1:].real @@ -159,8 +160,9 @@ class HarmonicTransformOperator(LinearOperator): def _slice_h2p(self, inp): res = np.empty((len(inp)+self.lmax+1)//2, dtype=(inp[0]*1j).dtype) - assert len(res) == ((self.mmax+1)*(self.mmax+2))//2 + \ - (self.mmax+1)*(self.lmax-self.mmax) + if len(res) != ((self.mmax+1)*(self.mmax+2))//2 + \ + (self.mmax+1)*(self.lmax-self.mmax): + raise ValueError("array length mismatch") res[0:self.lmax+1] = inp[0:self.lmax+1] res[self.lmax+1:] = np.sqrt(0.5)*(inp[self.lmax+1::2] + 1j*inp[self.lmax+2::2]) diff --git a/nifty4/spaces/power_space.py b/nifty4/spaces/power_space.py index f5b9d38a..c4ec0b54 100644 --- a/nifty4/spaces/power_space.py +++ b/nifty4/spaces/power_space.py @@ -62,7 +62,8 @@ class PowerSpace(Space): units of the harmonic partner space. """ 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) @staticmethod @@ -81,7 +82,8 @@ class PowerSpace(Space): units of the harmonic partner space. """ 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)), np.log(float(last_bound)), nbin-1, base=np.e) @@ -107,7 +109,8 @@ class PowerSpace(Space): nbin_max = int((dists[-1]-dists[0])/binsz_min)+2 if nbin is None: 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: raise ValueError("nbin is too large") if logarithmic: @@ -141,16 +144,19 @@ class PowerSpace(Space): tbb = binbounds locdat = np.searchsorted(tbb, dobj.local_data(k_length_array.val)) 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 temp_rho = np.bincount(dobj.local_data(temp_pindex).ravel(), minlength=nbin) 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 # sometimes returns its result as an integer array, even when # 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(), minlength=nbin).astype(np.float64, copy=False) temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho diff --git a/test/test_minimization/test_minimizers.py b/test/test_minimization/test_minimizers.py index 30b54e33..31460d29 100644 --- a/test/test_minimization/test_minimizers.py +++ b/test/test_minimization/test_minimizers.py @@ -41,7 +41,7 @@ class Test_Minimizers(unittest.TestCase): covariance = ift.DiagonalOperator(covariance_diagonal) 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: minimizer = minimizer_class(controller=IC) energy = ift.QuadraticEnergy(A=covariance, b=required_result, -- GitLab