Commit 2573b236 authored by Steininger, Theo (theos)'s avatar Steininger, Theo (theos)

Merge branch 'log_normal' into 'master'

Log normal

See merge request !177
parents 69db50ec 0eaf078f
Pipeline #15810 failed with stages
in 8 minutes and 3 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')
#
......@@ -23,7 +23,7 @@ from nifty.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):
......@@ -95,6 +95,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:
......
......@@ -831,6 +831,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(Loggable, object):
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 "
......
......@@ -111,7 +111,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.
......@@ -234,12 +236,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 xrange(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
......
......@@ -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):
......
......@@ -30,7 +30,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
......
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