diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py index 60cab0f6f3c18ba605ddaca5fff4d1c3bb81ef3e..79d445882b0bb9203a884fcfc841c73f69394880 100644 --- a/demos/critical_filtering.py +++ b/demos/critical_filtering.py @@ -117,9 +117,9 @@ if __name__ == "__main__": convergence_level=1, iteration_limit=5, callback=convergence_measure) - minimizer2 = VL_BFGS(convergence_tolerance=1e-4, + minimizer2 = VL_BFGS(convergence_tolerance=1e-10, convergence_level=1, - iteration_limit=20, + iteration_limit=30, callback=convergence_measure, max_history_length=20) minimizer3 = SteepestDescent(convergence_tolerance=1e-4, diff --git a/nifty/field.py b/nifty/field.py index d24b1af3470cbd427318cb66033d91f87b1aa7e1..a870dae37fd5e00720b025655e7a54eb87801ddf 100644 --- a/nifty/field.py +++ b/nifty/field.py @@ -599,6 +599,9 @@ class Field(Loggable, Versionable, object): if real_power: result = result_list[0] + if not issubclass(result_val_list[0].dtype.type, + np.complexfloating): + result = result.real else: result = result_list[0] + 1j*result_list[1] @@ -613,9 +616,17 @@ class Field(Loggable, Versionable, object): flipped_val = domain[space].hermitianize_inverter( x=flipped_val, axes=domain_axes[space]) - flipped_val = flipped_val.conjugate() - h = (val + flipped_val)/2. - a = val - h + # if no flips at all where performed `h` is a real field. + # if all spaces use the default implementation of doing nothing when + # no flips are applied, one can use `is` to infer this case. + + if flipped_val is val: + h = flipped_val.real + a = 1j * flipped_val.imag + else: + flipped_val = flipped_val.conjugate() + h = (val + flipped_val)/2. + a = val - h # correct variance if preserve_gaussian_variance: diff --git a/nifty/minimization/line_searching/line_search_strong_wolfe.py b/nifty/minimization/line_searching/line_search_strong_wolfe.py index bdff43d3703a09a388ed3bc2fc1e58b866c79028..eb5589817e4877483674204a7aa0002add26db73 100644 --- a/nifty/minimization/line_searching/line_search_strong_wolfe.py +++ b/nifty/minimization/line_searching/line_search_strong_wolfe.py @@ -116,7 +116,7 @@ class LineSearchStrongWolfe(LineSearch): if self.preferred_initial_step_size is not None: alpha1 = self.preferred_initial_step_size - elif old_phi_0 is not None and phiprime_0 != 0: + elif old_phi_0 is not None: alpha1 = min(1.0, 1.01*2*(phi_0 - old_phi_0)/phiprime_0) if alpha1 < 0: alpha1 = 1.0 diff --git a/nifty/operators/laplace_operator/laplace_operator.py b/nifty/operators/laplace_operator/laplace_operator.py index 00f6eba8fca1464454e08b8cdf8c0ec456a66086..e3b40d93ae0c074e122392916be445711e752730 100644 --- a/nifty/operators/laplace_operator/laplace_operator.py +++ b/nifty/operators/laplace_operator/laplace_operator.py @@ -20,17 +20,17 @@ import numpy as np from nifty.field import Field from nifty.spaces.power_space import PowerSpace from nifty.operators.endomorphic_operator import EndomorphicOperator - +from nifty import sqrt import nifty.nifty_utilities as utilities class LaplaceOperator(EndomorphicOperator): """A irregular LaplaceOperator with free boundary and excluding monopole. - This LaplaceOperator implements the second derivative of a Field in PowerSpace - on logarithmic or linear scale with vanishing curvature at the boundary, starting - at the second entry of the Field. The second derivative of the Field on the irregular grid - is calculated using finite differences. + This LaplaceOperator implements the second derivative of a Field in + PowerSpace on logarithmic or linear scale with vanishing curvature at the + boundary, starting at the second entry of the Field. The second derivative + of the Field on the irregular grid is calculated using finite differences. Parameters ---------- @@ -50,15 +50,19 @@ class LaplaceOperator(EndomorphicOperator): self._logarithmic = bool(logarithmic) + pos = self.domain[0].kindex.copy() if self.logarithmic: - self.positions = self.domain[0].kindex.copy() - self.positions[1:] = np.log(self.positions[1:]) - self.positions[0] = -1. - else: - self.positions = self.domain[0].kindex.copy() - self.positions[0] = -1 - - self.fwd_dist = self.positions[1:] - self.positions[:-1] + pos[1:] = np.log(pos[1:]) + pos[0] = pos[1]-1. + + self._dpos = pos[1:]-pos[:-1] # defined between points + # centered distances (also has entries for the first and last point + # for convenience, but they will never affect the result) + self._dposc = np.empty_like(pos) + self._dposc[:-1] = self._dpos + self._dposc[-1] = 0. + self._dposc[1:] += self._dpos + self._dposc *= 0.5 @property def target(self): @@ -94,24 +98,20 @@ class LaplaceOperator(EndomorphicOperator): else: axes = x.domain_axes[spaces[0]] axis = axes[0] + nval = len(self._dposc) prefix = (slice(None),) * axis - fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape) - positions = self.positions.reshape((1,)*axis + self.positions.shape) + sl_l = prefix + (slice(None, -1),) # "left" slice + sl_r = prefix + (slice(1, None),) # "right" slice + dpos = self._dpos.reshape((1,)*axis + (nval-1,)) + dposc = self._dposc.reshape((1,)*axis + (nval,)) + deriv = (x.val[sl_r]-x.val[sl_l])/dpos # defined between points ret = x.val.copy_empty() - x = x.val - ret[prefix + (slice(1, -1),)] = \ - (-((x[prefix + (slice(1, -1),)] - x[prefix + (slice(0, -2),)]) / - fwd_dist[prefix + (slice(0, -1),)]) + - ((x[prefix + (slice(2, None),)] - x[prefix + (slice(1, -1),)]) / - fwd_dist[prefix + (slice(1, None),)])) - ret[prefix + (slice(1, -1),)] /= \ - (positions[prefix + (slice(2, None),)] - - positions[prefix + (slice(None, -2),)]) - ret *= 2. - ret[prefix + (slice(0, 2),)] = 0 - ret[prefix + (slice(-1, -1),)] = 0 - ret[prefix + (slice(2, None),)] *= \ - np.sqrt(fwd_dist)[prefix + (slice(1, None),)] + ret[sl_l] = deriv + ret[prefix + (-1,)] = 0. + ret[sl_r] -= deriv + ret /= sqrt(dposc) + ret[prefix + (slice(None, 2),)] = 0. + ret[prefix + (-1,)] = 0. return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces) def _adjoint_times(self, x, spaces): @@ -124,42 +124,19 @@ class LaplaceOperator(EndomorphicOperator): else: axes = x.domain_axes[spaces[0]] axis = axes[0] + nval = len(self._dposc) prefix = (slice(None),) * axis - fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape) - positions = self.positions.reshape((1,)*axis + self.positions.shape) + sl_l = prefix + (slice(None, -1),) # "left" slice + sl_r = prefix + (slice(1, None),) # "right" slice + dpos = self._dpos.reshape((1,)*axis + (nval-1,)) + dposc = self._dposc.reshape((1,)*axis + (nval,)) y = x.copy().weight(power=0.5).val - y[prefix + (slice(2, None),)] *= \ - np.sqrt(fwd_dist)[prefix + (slice(1, None),)] - y[prefix + (slice(0, 2),)] = 0 - y[prefix + (slice(-1, -1),)] = 0 - ret = y.copy_empty() - y[prefix + (slice(1, -1),)] /= \ - (positions[prefix + (slice(2, None),)] - - positions[prefix + (slice(None, -2),)]) - y *= 2 - ret[prefix + (slice(1, -1),)] = \ - (-y[prefix + (slice(1, -1),)]/fwd_dist[prefix + (slice(0, -1),)] - - y[prefix + (slice(1, -1),)]/fwd_dist[prefix + (slice(1, None),)]) - ret[prefix + (slice(0, -2),)] += \ - y[prefix + (slice(1, -1),)] / fwd_dist[prefix + (slice(0, -1),)] - ret[prefix + (slice(2, None),)] += \ - y[prefix + (slice(1, -1),)] / fwd_dist[prefix + (slice(1, None),)] + y /= sqrt(dposc) + y[prefix + (slice(None, 2),)] = 0. + y[prefix + (-1,)] = 0. + deriv = (y[sl_r]-y[sl_l])/dpos # defined between points + ret = x.val.copy_empty() + ret[sl_l] = deriv + ret[prefix + (-1,)] = 0. + ret[sl_r] -= deriv return Field(self.domain, val=ret).weight(-1, spaces=spaces) - - def _irregular_laplace(self, x): - ret = np.zeros_like(x) - ret[1:-1] = (-(x[1:-1] - x[0:-2]) / self.fwd_dist[:-1] + - (x[2:] - x[1:-1]) / self.fwd_dist[1:]) - ret[1:-1] /= self.positions[2:] - self.positions[:-2] - ret *= 2. - return ret - - def _irregular_adj_laplace(self, x): - ret = np.zeros_like(x) - y = x.copy() - y[1:-1] /= self.positions[2:] - self.positions[:-2] - y *= 2 - ret[1:-1] = -y[1:-1] / self.fwd_dist[:-1] - y[1:-1] / self.fwd_dist[1:] - ret[0:-2] += y[1:-1] / self.fwd_dist[:-1] - ret[2:] += y[1:-1] / self.fwd_dist[1:] - return ret diff --git a/nifty/sugar.py b/nifty/sugar.py index 3048cc58300d9f44d6799ac23bcc711fbd4f158f..3a77bfcb95d754eb18f92d663fb598da01d25084 100644 --- a/nifty/sugar.py +++ b/nifty/sugar.py @@ -111,10 +111,12 @@ def generate_posterior_sample(mean, covariance): power = S.diagonal().power_analyze()**.5 mock_signal = power.power_synthesize(real_signal=True) - noise = N.diagonal(bare=True).val + noise = N.diagonal(bare=True) mock_noise = Field.from_random(random_type="normal", domain=N.domain, - std=sqrt(noise), dtype=noise.dtype) + dtype=noise.dtype) + mock_noise *= sqrt(noise) + mock_data = R(mock_signal) + mock_noise mock_j = R.adjoint_times(N.inverse_times(mock_data)) diff --git a/test/test_operators/test_laplace_operator.py b/test/test_operators/test_laplace_operator.py new file mode 100644 index 0000000000000000000000000000000000000000..a2591dccaf75df24f2693c2e8a3a431a7e3b05e8 --- /dev/null +++ b/test/test_operators/test_laplace_operator.py @@ -0,0 +1,35 @@ +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# Copyright(C) 2013-2017 Max-Planck-Society +# +# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik +# and financially supported by the Studienstiftung des deutschen Volkes. + +import unittest +import numpy as np +import nifty as ift +from numpy.testing import assert_allclose +from itertools import product +from test.common import expand + + +class LaplaceOperatorTests(unittest.TestCase): + @expand(product([None, False, True], [False, True], [10, 100, 1000])) + def test_Laplace(self, log1, log2, sz): + s = ift.RGSpace(sz, harmonic=True) + p = ift.PowerSpace(s, logarithmic=log1) + L = ift.LaplaceOperator(p, logarithmic=log2) + arr = np.random.random(p.shape[0]) + fp = ift.Field(p, val=arr) + assert_allclose(L(fp).vdot(L(fp)), L.adjoint_times(L(fp)).vdot(fp))