Commit f8e752b6 authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'master' into tweak_limited_exp

parents 104dc2d7 e0f1f873
Pipeline #16976 passed with stage
in 51 minutes and 32 seconds
...@@ -117,9 +117,9 @@ if __name__ == "__main__": ...@@ -117,9 +117,9 @@ if __name__ == "__main__":
convergence_level=1, convergence_level=1,
iteration_limit=5, iteration_limit=5,
callback=convergence_measure) callback=convergence_measure)
minimizer2 = VL_BFGS(convergence_tolerance=1e-4, minimizer2 = VL_BFGS(convergence_tolerance=1e-10,
convergence_level=1, convergence_level=1,
iteration_limit=20, iteration_limit=30,
callback=convergence_measure, callback=convergence_measure,
max_history_length=20) max_history_length=20)
minimizer3 = SteepestDescent(convergence_tolerance=1e-4, minimizer3 = SteepestDescent(convergence_tolerance=1e-4,
......
...@@ -477,7 +477,7 @@ class Field(Loggable, Versionable, object): ...@@ -477,7 +477,7 @@ class Field(Loggable, Versionable, object):
""" Yields a sampled field with `self`**2 as its power spectrum. """ Yields a sampled field with `self`**2 as its power spectrum.
This method draws a Gaussian random field in the harmonic partner This method draws a Gaussian random field in the harmonic partner
domain of this fields domains, using this field as power spectrum. domain of this field's domains, using this field as power spectrum.
Parameters Parameters
---------- ----------
...@@ -599,6 +599,9 @@ class Field(Loggable, Versionable, object): ...@@ -599,6 +599,9 @@ class Field(Loggable, Versionable, object):
if real_power: if real_power:
result = result_list[0] result = result_list[0]
if not issubclass(result_val_list[0].dtype.type,
np.complexfloating):
result = result.real
else: else:
result = result_list[0] + 1j*result_list[1] result = result_list[0] + 1j*result_list[1]
...@@ -613,9 +616,17 @@ class Field(Loggable, Versionable, object): ...@@ -613,9 +616,17 @@ class Field(Loggable, Versionable, object):
flipped_val = domain[space].hermitianize_inverter( flipped_val = domain[space].hermitianize_inverter(
x=flipped_val, x=flipped_val,
axes=domain_axes[space]) axes=domain_axes[space])
flipped_val = flipped_val.conjugate() # if no flips at all where performed `h` is a real field.
h = (val + flipped_val)/2. # if all spaces use the default implementation of doing nothing when
a = val - h # 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 # correct variance
if preserve_gaussian_variance: if preserve_gaussian_variance:
...@@ -696,7 +707,7 @@ class Field(Loggable, Versionable, object): ...@@ -696,7 +707,7 @@ class Field(Loggable, Versionable, object):
# ---Properties--- # ---Properties---
def set_val(self, new_val=None, copy=False): def set_val(self, new_val=None, copy=False):
""" Sets the fields distributed_data_object. """ Sets the field's distributed_data_object.
Parameters Parameters
---------- ----------
...@@ -869,7 +880,7 @@ class Field(Loggable, Versionable, object): ...@@ -869,7 +880,7 @@ class Field(Loggable, Versionable, object):
dtype : type dtype : type
The datatype the output shall have. This can be used to override The datatype the output shall have. This can be used to override
the fields dtype. the field's dtype.
Returns Returns
------- -------
......
...@@ -116,7 +116,7 @@ class LineSearchStrongWolfe(LineSearch): ...@@ -116,7 +116,7 @@ class LineSearchStrongWolfe(LineSearch):
if self.preferred_initial_step_size is not None: if self.preferred_initial_step_size is not None:
alpha1 = self.preferred_initial_step_size 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) alpha1 = min(1.0, 1.01*2*(phi_0 - old_phi_0)/phiprime_0)
if alpha1 < 0: if alpha1 < 0:
alpha1 = 1.0 alpha1 = 1.0
......
...@@ -20,17 +20,17 @@ import numpy as np ...@@ -20,17 +20,17 @@ import numpy as np
from nifty.field import Field from nifty.field import Field
from nifty.spaces.power_space import PowerSpace from nifty.spaces.power_space import PowerSpace
from nifty.operators.endomorphic_operator import EndomorphicOperator from nifty.operators.endomorphic_operator import EndomorphicOperator
from nifty import sqrt
import nifty.nifty_utilities as utilities import nifty.nifty_utilities as utilities
class LaplaceOperator(EndomorphicOperator): class LaplaceOperator(EndomorphicOperator):
"""A irregular LaplaceOperator with free boundary and excluding monopole. """A irregular LaplaceOperator with free boundary and excluding monopole.
This LaplaceOperator implements the second derivative of a Field in PowerSpace This LaplaceOperator implements the second derivative of a Field in
on logarithmic or linear scale with vanishing curvature at the boundary, starting PowerSpace on logarithmic or linear scale with vanishing curvature at the
at the second entry of the Field. The second derivative of the Field on the irregular grid boundary, starting at the second entry of the Field. The second derivative
is calculated using finite differences. of the Field on the irregular grid is calculated using finite differences.
Parameters Parameters
---------- ----------
...@@ -50,15 +50,19 @@ class LaplaceOperator(EndomorphicOperator): ...@@ -50,15 +50,19 @@ class LaplaceOperator(EndomorphicOperator):
self._logarithmic = bool(logarithmic) self._logarithmic = bool(logarithmic)
pos = self.domain[0].kindex.copy()
if self.logarithmic: if self.logarithmic:
self.positions = self.domain[0].kindex.copy() pos[1:] = np.log(pos[1:])
self.positions[1:] = np.log(self.positions[1:]) pos[0] = pos[1]-1.
self.positions[0] = -1.
else: self._dpos = pos[1:]-pos[:-1] # defined between points
self.positions = self.domain[0].kindex.copy() # centered distances (also has entries for the first and last point
self.positions[0] = -1 # for convenience, but they will never affect the result)
self._dposc = np.empty_like(pos)
self.fwd_dist = self.positions[1:] - self.positions[:-1] self._dposc[:-1] = self._dpos
self._dposc[-1] = 0.
self._dposc[1:] += self._dpos
self._dposc *= 0.5
@property @property
def target(self): def target(self):
...@@ -94,24 +98,20 @@ class LaplaceOperator(EndomorphicOperator): ...@@ -94,24 +98,20 @@ class LaplaceOperator(EndomorphicOperator):
else: else:
axes = x.domain_axes[spaces[0]] axes = x.domain_axes[spaces[0]]
axis = axes[0] axis = axes[0]
nval = len(self._dposc)
prefix = (slice(None),) * axis prefix = (slice(None),) * axis
fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape) sl_l = prefix + (slice(None, -1),) # "left" slice
positions = self.positions.reshape((1,)*axis + self.positions.shape) 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() ret = x.val.copy_empty()
x = x.val ret[sl_l] = deriv
ret[prefix + (slice(1, -1),)] = \ ret[prefix + (-1,)] = 0.
(-((x[prefix + (slice(1, -1),)] - x[prefix + (slice(0, -2),)]) / ret[sl_r] -= deriv
fwd_dist[prefix + (slice(0, -1),)]) + ret /= sqrt(dposc)
((x[prefix + (slice(2, None),)] - x[prefix + (slice(1, -1),)]) / ret[prefix + (slice(None, 2),)] = 0.
fwd_dist[prefix + (slice(1, None),)])) ret[prefix + (-1,)] = 0.
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),)]
return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces) return Field(self.domain, val=ret).weight(power=-0.5, spaces=spaces)
def _adjoint_times(self, x, spaces): def _adjoint_times(self, x, spaces):
...@@ -124,42 +124,19 @@ class LaplaceOperator(EndomorphicOperator): ...@@ -124,42 +124,19 @@ class LaplaceOperator(EndomorphicOperator):
else: else:
axes = x.domain_axes[spaces[0]] axes = x.domain_axes[spaces[0]]
axis = axes[0] axis = axes[0]
nval = len(self._dposc)
prefix = (slice(None),) * axis prefix = (slice(None),) * axis
fwd_dist = self.fwd_dist.reshape((1,)*axis + self.fwd_dist.shape) sl_l = prefix + (slice(None, -1),) # "left" slice
positions = self.positions.reshape((1,)*axis + self.positions.shape) 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 = x.copy().weight(power=0.5).val
y[prefix + (slice(2, None),)] *= \ y /= sqrt(dposc)
np.sqrt(fwd_dist)[prefix + (slice(1, None),)] y[prefix + (slice(None, 2),)] = 0.
y[prefix + (slice(0, 2),)] = 0 y[prefix + (-1,)] = 0.
y[prefix + (slice(-1, -1),)] = 0 deriv = (y[sl_r]-y[sl_l])/dpos # defined between points
ret = y.copy_empty() ret = x.val.copy_empty()
y[prefix + (slice(1, -1),)] /= \ ret[sl_l] = deriv
(positions[prefix + (slice(2, None),)] - ret[prefix + (-1,)] = 0.
positions[prefix + (slice(None, -2),)]) ret[sl_r] -= deriv
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),)]
return Field(self.domain, val=ret).weight(-1, spaces=spaces) 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
...@@ -111,10 +111,12 @@ def generate_posterior_sample(mean, covariance): ...@@ -111,10 +111,12 @@ def generate_posterior_sample(mean, covariance):
power = S.diagonal().power_analyze()**.5 power = S.diagonal().power_analyze()**.5
mock_signal = power.power_synthesize(real_signal=True) 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, 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_data = R(mock_signal) + mock_noise
mock_j = R.adjoint_times(N.inverse_times(mock_data)) mock_j = R.adjoint_times(N.inverse_times(mock_data))
......
# 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 <http://www.gnu.org/licenses/>.
#
# 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))
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