Commit 34ef9abd authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge master

parents 8d6c5a9e 47567c7e
Pipeline #15813 failed with stage
in 7 minutes and 48 seconds
# -*- coding: utf-8 -*-
from nifty import *
if __name__ == "__main__":
nifty_configuration['default_distribution_strategy'] = 'fftw'
# Setting up parameters |\label{code:wf_parameters}|
correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space
response_sigma = 0.02 # Smoothing length of response (in same unit as L)
signal_to_noise = 100 # The signal to noise ratio
np.random.seed(43) # Fixing the random seed
def power_spectrum(k): # Defining the power spectrum
a = 4 * correlation_length * field_variance**2
return a / (1 + k * correlation_length) ** 4
# Setting up the geometry |\label{code:wf_geometry}|
L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis)
#signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = HPSpace(16)
harmonic_space = FFTOperator.get_default_codomain(signal_space)
fft = FFTOperator(harmonic_space, target=signal_space, target_dtype=np.float)
power_space = PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum)
mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response
mask = Field(signal_space, val=1.)
N10 = int(N_pixels/10)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization
N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
# Wiener filter
m0 = Field(harmonic_space, val=0.j)
energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S)
minimizer1 = VL_BFGS(convergence_tolerance=1e-5,
iteration_limit=3000,
#callback=convergence_measure,
max_history_length=20)
minimizer2 = RelaxedNewton(convergence_tolerance=1e-5,
iteration_limit=10,
#callback=convergence_measure
)
minimizer3 = SteepestDescent(convergence_tolerance=1e-5, iteration_limit=1000)
# me1 = minimizer1(energy)
# me2 = minimizer2(energy)
# me3 = minimizer3(energy)
# m1 = fft(me1[0].position)
# m2 = fft(me2[0].position)
# m3 = fft(me3[0].position)
#
# # Probing the variance
# class Proby(DiagonalProberMixin, Prober): pass
# proby = Proby(signal_space, probe_count=100)
# proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
#
# sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1))
#Plotting #|\label{code:wf_plotting}|
#plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap())
plotter.figure.xaxis = plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = plotting.Axis(label='Pixel Index')
plotter.plot.zmax = 5; plotter.plot.zmin = -5
## plotter(variance, path = 'variance.html')
# #plotter.plot.zmin = exp(mock_signal.min());
# plotter(mock_signal.real, path='mock_signal.html')
# plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
# path = 'log_of_data.html')
#
# plotter(m1.real, path='m_LBFGS.html')
# plotter(m2.real, path='m_Newton.html')
# plotter(m3.real, path='m_SteepestDescent.html')
#
......@@ -24,7 +24,7 @@ from .field import Field
__all__ = ['cos', 'sin', 'cosh', 'sinh', 'tan', 'tanh', 'arccos', 'arcsin',
'arccosh', 'arcsinh', 'arctan', 'arctanh', 'sqrt', 'exp', 'log',
'conjugate']
'conjugate', 'clipped_exp']
def _math_helper(x, function):
......@@ -96,6 +96,10 @@ def exp(x):
return _math_helper(x, np.exp)
def clipped_exp(x):
return _math_helper(x, lambda z: np.exp(np.minimum(200, z)))
def log(x, base=None):
result = _math_helper(x, np.log)
if base is not None:
......
......@@ -131,7 +131,6 @@ class Field(Loggable, Versionable, object):
else:
self.set_val(new_val=val, copy=copy)
def _parse_domain(self, domain, val=None):
if domain is None:
if isinstance(val, Field):
......@@ -244,6 +243,15 @@ class Field(Loggable, Versionable, object):
# random number generator to it
sample = f.get_val(copy=False)
generator_function = getattr(Random, random_type)
comm = sample.comm
size = comm.size
if (sample.distribution_strategy in DISTRIBUTION_STRATEGIES['not'] and
size > 1):
seed = np.random.randint(10000000)
seed = comm.bcast(seed, root=0)
np.random.seed(seed)
sample.apply_generator(
lambda shape: generator_function(dtype=f.dtype,
shape=shape,
......@@ -835,6 +843,24 @@ class Field(Loggable, Versionable, object):
except TypeError:
return 0.
@property
def real(self):
""" The real part of the field (data is not copied).
"""
real_part = self.val.real
result = self.copy_empty(dtype=real_part.dtype)
result.set_val(new_val=real_part, copy=False)
return result
@property
def imag(self):
""" The imaginary part of the field (data is not copied).
"""
real_part = self.val.imag
result = self.copy_empty(dtype=real_part.dtype)
result.set_val(new_val=real_part, copy=False)
return result
# ---Special unary/binary operations---
def cast(self, x=None, dtype=None):
......
from .critical_filter import *
from .log_normal_wiener_filter import *
from .wiener_filter import *
import numpy as np
from nifty.energies.energy import Energy
from nifty.operators.smoothness_operator import SmoothnessOperator
......
# -*- coding: utf-8 -*-
from log_normal_wiener_filter_curvature import LogNormalWienerFilterCurvature
from log_normal_wiener_filter_energy import LogNormalWienerFilterEnergy
from nifty.operators import EndomorphicOperator,\
InvertibleOperatorMixin
from nifty.energies.memoization import memo
from nifty.basic_arithmetics import clipped_exp
from nifty.sugar import create_composed_fft_operator
class LogNormalWienerFilterCurvature(InvertibleOperatorMixin,
EndomorphicOperator):
"""The curvature of the LogNormalWienerFilterEnergy.
This operator implements the second derivative of the
LogNormalWienerFilterEnergy used in some minimization algorithms or for
error estimates of the posterior maps. It is the inverse of the propagator
operator.
Parameters
----------
R: LinearOperator,
The response operator of the Wiener filter measurement.
N : EndomorphicOperator
The noise covariance.
S: DiagonalOperator,
The prior signal covariance
"""
def __init__(self, R, N, S, d, position, inverter=None,
preconditioner=None, fft4exp=None, **kwargs):
self._cache = {}
self.R = R
self.N = N
self.S = S
self.d = d
self.position = position
if preconditioner is None:
preconditioner = self.S.times
self._domain = self.S.domain
if fft4exp is None:
self._fft = create_composed_fft_operator(self.domain,
all_to='position')
else:
self._fft = fft4exp
super(LogNormalWienerFilterCurvature, self).__init__(
inverter=inverter,
preconditioner=preconditioner,
**kwargs)
@property
def domain(self):
return self._domain
@property
def self_adjoint(self):
return True
@property
def unitary(self):
return False
# ---Added properties and methods---
def _times(self, x, spaces):
part1 = self.S.inverse_times(x)
# part2 = self._exppRNRexppd * x
part3 = self._fft.adjoint_times(self._expp_sspace * self._fft(x))
part3 = self._fft.adjoint_times(
self._expp_sspace *
self._fft(self.R.adjoint_times(
self.N.inverse_times(self.R(part3)))))
return part1 + part3 # + part2
@property
@memo
def _expp_sspace(self):
return clipped_exp(self._fft(self.position))
@property
@memo
def _expp(self):
return self._fft.adjoint_times(self._expp_sspace)
@property
@memo
def _Rexppd(self):
return self.R(self._expp) - self.d
@property
@memo
def _NRexppd(self):
return self.N.inverse_times(self._Rexppd)
@property
@memo
def _exppRNRexppd(self):
return self._fft.adjoint_times(
self._expp_sspace *
self._fft(self.R.adjoint_times(self._NRexppd)))
from nifty.energies.energy import Energy
from nifty.energies.memoization import memo
from nifty.library.log_normal_wiener_filter import \
LogNormalWienerFilterCurvature
from nifty.sugar import create_composed_fft_operator
class LogNormalWienerFilterEnergy(Energy):
"""The Energy for the log-normal Wiener filter.
It covers the case of linear measurement with
Gaussian noise and Gaussain signal prior with known covariance.
Parameters
----------
position: Field,
The current position.
d : Field,
the data.
R : Operator,
The response operator, describtion of the measurement process.
N : EndomorphicOperator,
The noise covariance in data space.
S : EndomorphicOperator,
The prior signal covariance in harmonic space.
"""
def __init__(self, position, d, R, N, S, fft4exp=None):
super(LogNormalWienerFilterEnergy, self).__init__(position=position)
self.d = d
self.R = R
self.N = N
self.S = S
if fft4exp is None:
self._fft = create_composed_fft_operator(self.S.domain,
all_to='position')
else:
self._fft = fft4exp
def at(self, position):
return self.__class__(position=position, d=self.d, R=self.R, N=self.N,
S=self.S, fft4exp=self._fft)
@property
@memo
def value(self):
return 0.5*(self.position.vdot(self._Sp) +
self._Rexppd.vdot(self._NRexppd))
@property
@memo
def gradient(self):
return self._Sp + self._exppRNRexppd
@property
@memo
def curvature(self):
return LogNormalWienerFilterCurvature(R=self.R, N=self.N, S=self.S,
d=self.d, position=self.position,
fft4exp=self._fft)
@property
def _expp(self):
return self.curvature._expp
@property
def _Rexppd(self):
return self.curvature._Rexppd
@property
def _NRexppd(self):
return self.curvature._NRexppd
@property
def _exppRNRexppd(self):
return self.curvature._exppRNRexppd
@property
@memo
def _Sp(self):
return self.S.inverse_times(self.position)
......@@ -149,17 +149,23 @@ class DescentMinimizer(with_metaclass(NiftyMeta, type('NewBase', (Loggable, obje
descent_direction = self.get_descent_direction(energy)
# compute the step length, which minimizes energy.value along the
# search direction
step_length, f_k, new_energy = \
self.line_searcher.perform_line_search(
energy=energy,
pk=descent_direction,
f_k_minus_1=f_k_minus_1)
try:
step_length, f_k, new_energy = \
self.line_searcher.perform_line_search(
energy=energy,
pk=descent_direction,
f_k_minus_1=f_k_minus_1)
except RuntimeError:
self.logger.warn(
"Stopping because of RuntimeError in line-search")
break
if f_k_minus_1 is None:
delta=1e30
delta = 1e30
else:
delta = abs(f_k -f_k_minus_1)/max(abs(f_k),abs(f_k_minus_1),1.)
delta = (abs(f_k-f_k_minus_1) /
max(abs(f_k), abs(f_k_minus_1), 1.))
f_k_minus_1 = energy.value
tx1=energy.position-new_energy.position
# check if new energy value is bigger than old energy value
if (new_energy.value - energy.value) > 0:
self.logger.info("Line search algorithm returned a new energy "
......
......@@ -114,7 +114,9 @@ class LineSearchStrongWolfe(LineSearch):
le_0 = self.line_energy.at(0)
phi_0 = le_0.value
phiprime_0 = le_0.directional_derivative
assert phiprime_0<0, "input direction must be a descent direction"
if phiprime_0 >= 0:
self.logger.error("Input direction must be a descent direction")
raise RuntimeError
# set alphas
alpha0 = 0.
......@@ -237,12 +239,14 @@ class LineSearchStrongWolfe(LineSearch):
cubic_delta = 0.2 # cubic
quad_delta = 0.1 # quadratic
phiprime_alphaj = 0.
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.
assert phiprime_lo*(alpha_hi-alpha_lo) < 0.
for i in range(max_iterations):
#assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
#assert phiprime_lo*(alpha_hi-alpha_lo)<0.
# assert phi_lo <= phi_0 + self.c1*alpha_lo*phiprime_0
# assert phiprime_lo*(alpha_hi-alpha_lo)<0.
delta_alpha = alpha_hi - alpha_lo
if delta_alpha < 0:
a, b = alpha_hi, alpha_lo
......
......@@ -264,6 +264,10 @@ class MPIFFT(Transform):
if p.has_output:
result = p.output_array
if result.shape != val.shape:
raise ValueError("Output shape is different than input shape. "
"Maybe fftw tries to optimize the "
"bit-alignment? Try a different array-size.")
else:
return None
......
......@@ -32,7 +32,9 @@ class GLMollweide(Heatmap):
xsize=self.xsize,
color_map=self.color_map,
webgl=self.webgl,
smoothing=self.smoothing)
smoothing=self.smoothing,
zmin=self.zmin,
zmax=self.zmax)
@staticmethod
def _find_closest(A, target):
......
......@@ -31,7 +31,9 @@ class HPMollweide(Heatmap):
xsize=self.xsize,
color_map=self.color_map,
webgl=self.webgl,
smoothing=self.smoothing)
smoothing=self.smoothing,
zmin=self.zmin,
zmax=self.zmax)
def _mollview(self, x):
xsize = self.xsize
......
import unittest
import numpy as np
from numpy.testing import assert_approx_equal
from nifty import Field,\
......@@ -10,7 +11,7 @@ from itertools import product
from test.common import expand
class ResponseOperator_Tests(unittest.TestCase):
spaces = [RGSpace(100)]
spaces = [RGSpace(128)]
@expand(product(spaces, [0., 5., 1.], [0., 1., .33] ))
def test_property(self, space, sigma, exposure):
......
......@@ -37,7 +37,7 @@ def _get_rtol(tp):
return 1e-5
class SmoothingOperator_Tests(unittest.TestCase):
spaces = [RGSpace(100)]
spaces = [RGSpace(128)]
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_property(self, space, sigma, log_distances):
......@@ -83,7 +83,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
tt2 = rand2.vdot(op.inverse_adjoint_times(rand1))
assert_approx_equal(tt1, tt2)
@expand(product([100, 200], [1, 0.4], [0., 1., 3.7],
@expand(product([128, 256], [1, 0.4], [0., 1., 3.7],
[np.float64, np.complex128]))
def test_smooth_regular1(self, sz, d, sigma, tp):
tol = _get_rtol(tp)
......
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