Commit f81ce5cf authored by Martin Reinecke's avatar Martin Reinecke

Merge branch 'master' into python3

parents 6be97691 e0f1f873
Pipeline #16928 passed with stage
in 20 minutes and 45 seconds
......@@ -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,
......
......@@ -481,7 +481,7 @@ class Field(Loggable, Versionable, object):
""" Yields a sampled field with `self`**2 as its power spectrum.
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
----------
......@@ -603,6 +603,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]
......@@ -617,9 +620,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:
......@@ -700,7 +711,7 @@ class Field(Loggable, Versionable, object):
# ---Properties---
def set_val(self, new_val=None, copy=False):
""" Sets the fields distributed_data_object.
""" Sets the field's distributed_data_object.
Parameters
----------
......@@ -873,7 +884,7 @@ class Field(Loggable, Versionable, object):
dtype : type
The datatype the output shall have. This can be used to override
the fields dtype.
the field's dtype.
Returns
-------
......
......@@ -119,7 +119,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
......
......@@ -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
......@@ -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))
......
# 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