Commit 7107323b authored by Martin Reinecke's avatar Martin Reinecke
Browse files

updates

parent 7fc8256f
import numpy as np import numpy as np
from nifty import (VL_BFGS, DiagonalOperator, FFTOperator, Field, import nifty2go as ift
LinearOperator, PowerSpace, RelaxedNewton, RGSpace,
SteepestDescent, create_power_operator, exp, log, sqrt)
from nifty.library.critical_filter import CriticalPowerEnergy
from nifty.library.wiener_filter import WienerFilterEnergy
import plotly.graph_objs as go import plotly.graph_objs as go
import plotly.offline as pl import plotly.offline as pl
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
np.random.seed(42) np.random.seed(42)
def plot_parameters(m, t, p, p_d): def plot_parameters(m, t, p, p_d):
x = ift.log(t.domain[0].kindex)
x = log(t.domain[0].kindex)
m = fft.adjoint_times(m) m = fft.adjoint_times(m)
m = m.val.get_full_data().real m = m.val.real
t = t.val.get_full_data().real t = t.val.real
p = p.val.get_full_data().real p = p.val.real
p_d = p_d.val.get_full_data().real p_d = p_d.val.real
pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False) pl.plot([go.Heatmap(z=m)], filename='map.html', auto_open=False)
pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p), pl.plot([go.Scatter(x=x, y=t), go.Scatter(x=x, y=p),
go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False) go.Scatter(x=x, y=p_d)], filename="t.html", auto_open=False)
class AdjointFFTResponse(LinearOperator): class AdjointFFTResponse(ift.LinearOperator):
def __init__(self, FFT, R, default_spaces=None): def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces) super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target self._domain = FFT.target
...@@ -56,34 +47,28 @@ class AdjointFFTResponse(LinearOperator): ...@@ -56,34 +47,28 @@ class AdjointFFTResponse(LinearOperator):
if __name__ == "__main__": if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space # Set up position space
s_space = RGSpace([128, 128]) s_space = ift.RGSpace([128, 128])
# s_space = HPSpace(32) # s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space # Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space) fft = ift.FFTOperator(s_space)
h_space = fft.target[0] h_space = fft.target[0]
# Set up power space # Set up power space
p_space = PowerSpace(h_space, logarithmic=True, p_space = ift.PowerSpace(h_space)
distribution_strategy=distribution_strategy)
# Choose the prior correlation structure and defining correlation operator # Choose the prior correlation structure and defining correlation operator
p_spec = (lambda k: (.5 / (k + 1) ** 3)) p_spec = (lambda k: (.5 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec, S = ift.create_power_operator(h_space, power_spectrum=p_spec)
distribution_strategy=distribution_strategy)
# Draw a sample sh from the prior distribution in harmonic space # Draw a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec, sp = ift.Field(p_space, val=p_spec(p_space.kindex))
distribution_strategy=distribution_strategy)
sh = sp.power_synthesize(real_signal=True) sh = sp.power_synthesize(real_signal=True)
# Choose the measurement instrument # Choose the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01) # Instrument = SmoothingOperator(s_space, sigma=0.01)
Instrument = DiagonalOperator(s_space, diagonal=1.) Instrument = ift.DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0 # Instrument._diagonal.val[200:400, 200:400] = 0
# Instrument._diagonal.val[64:512-64, 64:512-64] = 0 # Instrument._diagonal.val[64:512-64, 64:512-64] = 0
...@@ -91,10 +76,10 @@ if __name__ == "__main__": ...@@ -91,10 +76,10 @@ if __name__ == "__main__":
R = AdjointFFTResponse(fft, Instrument) R = AdjointFFTResponse(fft, Instrument)
noise = 1. noise = 1.
N = DiagonalOperator(s_space, diagonal=noise, bare=True) N = ift.DiagonalOperator(s_space, diagonal=noise, bare=True)
n = Field.from_random(domain=s_space, n = ift.Field.from_random(domain=s_space,
random_type='normal', random_type='normal',
std=sqrt(noise), std=ift.sqrt(noise),
mean=0) mean=0)
# Create mock data # Create mock data
...@@ -102,10 +87,9 @@ if __name__ == "__main__": ...@@ -102,10 +87,9 @@ if __name__ == "__main__":
# The information source # The information source
j = R.adjoint_times(N.inverse_times(d)) j = R.adjoint_times(N.inverse_times(d))
realized_power = log(sh.power_analyze(binbounds=p_space.binbounds)) realized_power = ift.log(sh.power_analyze(binbounds=p_space.binbounds))
data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds)) data_power = ift.log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real d_data = d.val.real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False) pl.plot([go.Heatmap(z=d_data)], filename='data.html', auto_open=False)
# Minimization strategy # Minimization strategy
...@@ -113,36 +97,33 @@ if __name__ == "__main__": ...@@ -113,36 +97,33 @@ if __name__ == "__main__":
x = a_energy.value x = a_energy.value
print(x, iteration) print(x, iteration)
minimizer1 = RelaxedNewton(convergence_tolerance=1e-4, IC1 = ift.DefaultIterationController(verbose=True,iteration_limit=5)
convergence_level=1, minimizer1 = ift.RelaxedNewton(IC1)
iteration_limit=5, IC2 = ift.DefaultIterationController(verbose=True,iteration_limit=30)
callback=convergence_measure) minimizer2 = ift.VL_BFGS(IC2, max_history_length=20)
minimizer2 = VL_BFGS(convergence_tolerance=1e-10, IC3 = ift.DefaultIterationController(verbose=True,iteration_limit=100)
convergence_level=1, minimizer3 = ift.SteepestDescent(IC3)
iteration_limit=30,
callback=convergence_measure,
max_history_length=20)
minimizer3 = SteepestDescent(convergence_tolerance=1e-4,
iteration_limit=100,
callback=convergence_measure)
# Set starting position # Set starting position
flat_power = Field(p_space, val=1e-8) flat_power = ift.Field(p_space, val=1e-8)
m0 = flat_power.power_synthesize(real_signal=True) m0 = flat_power.power_synthesize(real_signal=True)
t0 = Field(p_space, val=log(1./(1+p_space.kindex)**2)) def ps0(k):
return (1./(1.+k)**2)
t0 = ift.Field(p_space, val=ift.log(1./(1+p_space.kindex)**2))
for i in range(500): for i in range(500):
S0 = create_power_operator(h_space, power_spectrum=exp(t0), S0 = ift.create_power_operator(h_space, power_spectrum=ps0)
distribution_strategy=distribution_strategy)
# Initialize non-linear Wiener Filter energy # Initialize non-linear Wiener Filter energy
map_energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0) ICI = ift.DefaultIterationController(verbose=True,iteration_limit=60)
inverter = ift.ConjugateGradient(controller=ICI,preconditioner=S0.times)
map_energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S0, inverter=inverter)
# Solve the Wiener Filter analytically # Solve the Wiener Filter analytically
D0 = map_energy.curvature D0 = map_energy.curvature
m0 = D0.inverse_times(j) m0 = D0.inverse_times(j)
# Initialize power energy with updated parameters # Initialize power energy with updated parameters
power_energy = CriticalPowerEnergy(position=t0, m=m0, D=D0, power_energy = ift.library.CriticalPowerEnergy(position=t0, m=m0, D=D0,
smoothness_prior=10., samples=3) smoothness_prior=10., samples=3)
(power_energy, convergence) = minimizer2(power_energy) (power_energy, convergence) = minimizer2(power_energy)
...@@ -153,4 +134,4 @@ if __name__ == "__main__": ...@@ -153,4 +134,4 @@ if __name__ == "__main__":
# Plot current estimate # Plot current estimate
print(i) print(i)
if i % 5 == 0: if i % 5 == 0:
plot_parameters(m0, t0, log(sp), data_power) plot_parameters(m0, t0, ift.log(sp), data_power)
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from nifty import * import nifty2go as ift
import nifty as ift import numpy as np
if __name__ == "__main__": if __name__ == "__main__":
np.random.seed(42)
# Setting up parameters |\label{code:wf_parameters}| # Setting up parameters |\label{code:wf_parameters}|
correlation_length = 1. # Typical distance over which the field is correlated correlation_length = 1. # Typical distance over which the field is correlated
field_variance = 2. # Variance of field in position space field_variance = 2. # Variance of field in position space
...@@ -19,74 +19,73 @@ if __name__ == "__main__": ...@@ -19,74 +19,73 @@ if __name__ == "__main__":
L = 2. # Total side-length of the domain L = 2. # Total side-length of the domain
N_pixels = 128 # Grid resolution (pixels per axis) N_pixels = 128 # Grid resolution (pixels per axis)
#signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels) #signal_space = RGSpace([N_pixels, N_pixels], distances=L/N_pixels)
signal_space = HPSpace(16) signal_space = ift.HPSpace(16)
harmonic_space = FFTOperator.get_default_codomain(signal_space) harmonic_space = signal_space.get_default_codomain()
fft = FFTOperator(harmonic_space, target=signal_space) fft = ift.FFTOperator(harmonic_space, target=signal_space)
power_space = PowerSpace(harmonic_space) power_space = ift.PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}| # Creating the mock signal |\label{code:wf_mock_signal}|
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum) S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum) mock_power = ift.Field(power_space, val=power_spectrum(power_space.kindex))
mock_signal = fft(mock_power.power_synthesize(real_signal=True)) mock_signal = fft(mock_power.power_synthesize(real_signal=True))
# Setting up an exemplary response # Setting up an exemplary response
mask = Field(signal_space, val=1.) mask = ift.Field(signal_space, val=1.)
N10 = int(N_pixels/10) N10 = int(N_pixels/10)
#mask.val[N10*5:N10*9, N10*5:N10*9] = 0. #mask.val[N10*5:N10*9, N10*5:N10*9] = 0.
R = ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}| R = ift.ResponseOperator(signal_space, sigma=(response_sigma,), exposure=(mask,)) #|\label{code:wf_response}|
data_domain = R.target[0] data_domain = R.target[0]
R_harmonic = ComposedOperator([fft, R], default_spaces=[0, 0]) R_harmonic = ift.ComposedOperator([fft, R], default_spaces=[0, 0])
# Setting up the noise covariance and drawing a random noise realization # Setting up the noise covariance and drawing a random noise realization
N = DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True) N = ift.DiagonalOperator(data_domain, diagonal=mock_signal.var()/signal_to_noise, bare=True)
noise = Field.from_random(domain=data_domain, random_type='normal', noise = ift.Field.from_random(domain=data_domain, random_type='normal',
std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0) std=mock_signal.std()/np.sqrt(signal_to_noise), mean=0)
data = R(exp(mock_signal)) + noise #|\label{code:wf_mock_data}| data = R(ift.exp(mock_signal)) + noise #|\label{code:wf_mock_data}|
# Wiener filter # Wiener filter
m0 = Field(harmonic_space, val=0.j) m0 = ift.Field(harmonic_space, val=0.)
ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=1) ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=1)
ctrl2 = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1, name="outer") ctrl2 = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl, preconditioner=S.times) inverter = ift.ConjugateGradient(controller=ctrl, preconditioner=S.times)
energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter) energy = ift.library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter)
minimizer1 = VL_BFGS(controller=ctrl2,max_history_length=20) minimizer1 = ift.VL_BFGS(controller=ctrl2,max_history_length=20)
minimizer2 = RelaxedNewton(controller=ctrl2) minimizer2 = ift.RelaxedNewton(controller=ctrl2)
minimizer3 = SteepestDescent(controller=ctrl2) minimizer3 = ift.SteepestDescent(controller=ctrl2)
print type(energy.value)
# me1 = minimizer1(energy) me1 = minimizer1(energy)
# me2 = minimizer2(energy) me2 = minimizer2(energy)
me3 = minimizer3(energy) me3 = minimizer3(energy)
# m1 = fft(me1[0].position) m1 = fft(me1[0].position)
# m2 = fft(me2[0].position) m2 = fft(me2[0].position)
# m3 = fft(me3[0].position) m3 = fft(me3[0].position)
#
# # Probing the variance # Probing the variance
# class Proby(DiagonalProberMixin, Prober): pass class Proby(ift.DiagonalProberMixin, ift.Prober): pass
# proby = Proby(signal_space, probe_count=100) proby = Proby(signal_space, probe_count=100)
# proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z)))) proby(lambda z: fft(wiener_curvature.inverse_times(fft.inverse_times(z))))
#
# sm = SmoothingOperator(signal_space, sigma=0.02) sm = SmoothingOperator(signal_space, sigma=0.02)
# variance = sm(proby.diagonal.weight(-1)) variance = sm(proby.diagonal.weight(-1))
#Plotting #|\label{code:wf_plotting}| #Plotting #|\label{code:wf_plotting}|
#plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap()) #plotter = plotting.RG2DPlotter(color_map=plotting.colormaps.PlankCmap())
plotter = plotting.HealpixPlotter(color_map=plotting.colormaps.PlankCmap()) plotter = ift.plotting.HealpixPlotter(color_map=ift.plotting.colormaps.PlankCmap())
plotter.figure.xaxis = plotting.Axis(label='Pixel Index') plotter.figure.xaxis = ift.plotting.Axis(label='Pixel Index')
plotter.figure.yaxis = plotting.Axis(label='Pixel Index') plotter.figure.yaxis = ift.plotting.Axis(label='Pixel Index')
plotter.plot.zmax = 5; plotter.plot.zmin = -5 plotter.plot.zmax = 5; plotter.plot.zmin = -5
## plotter(variance, path = 'variance.html') plotter(variance, path = 'variance.html')
# #plotter.plot.zmin = exp(mock_signal.min()); # #plotter.plot.zmin = exp(mock_signal.min());
# plotter(mock_signal.real, path='mock_signal.html') # 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)), # plotter(Field(signal_space, val=np.log(data.val.get_full_data().real).reshape(signal_space.shape)),
# path = 'log_of_data.html') # path = 'log_of_data.html')
# #
# plotter(m1.real, path='m_LBFGS.html') plotter(m1.real, path='m_LBFGS.html')
# plotter(m2.real, path='m_Newton.html') plotter(m2.real, path='m_Newton.html')
# plotter(m3.real, path='m_SteepestDescent.html') plotter(m3.real, path='m_SteepestDescent.html')
# #
...@@ -37,7 +37,7 @@ class Field(object): ...@@ -37,7 +37,7 @@ class Field(object):
In NIFTY, Fields are used to store data arrays and carry all the needed In NIFTY, Fields are used to store data arrays and carry all the needed
metainformation (i.e. the domain) for operators to be able to work on them. metainformation (i.e. the domain) for operators to be able to work on them.
In addition Field has methods to work with power-spectra. In addition, Field has methods to work with power spectra.
Parameters Parameters
---------- ----------
...@@ -52,7 +52,7 @@ class Field(object): ...@@ -52,7 +52,7 @@ class Field(object):
array's dimensions must match the domain's. array's dimensions must match the domain's.
dtype : type dtype : type
A numpy.type. Most common are int, float and complex. A numpy.type. Most common are float and complex.
copy: boolean copy: boolean
...@@ -265,11 +265,9 @@ class Field(object): ...@@ -265,11 +265,9 @@ class Field(object):
for part in parts] for part in parts]
if keep_phase_information: if keep_phase_information:
result_field = parts[0] + 1j*parts[1] return parts[0] + 1j*parts[1]
else: else:
result_field = parts[0] return parts[0]
return result_field
@classmethod @classmethod
def _single_power_analyze(cls, work_field, space_index, binbounds): def _single_power_analyze(cls, work_field, space_index, binbounds):
...@@ -295,21 +293,14 @@ class Field(object): ...@@ -295,21 +293,14 @@ class Field(object):
# create the result field and put power_spectrum into it # create the result field and put power_spectrum into it
result_domain = list(work_field.domain) result_domain = list(work_field.domain)
result_domain[space_index] = power_domain result_domain[space_index] = power_domain
result_dtype = power_spectrum.dtype
result_field = work_field.copy_empty(
domain=result_domain,
dtype=result_dtype)
result_field.set_val(new_val=power_spectrum, copy=False)
return result_field return Field(domain=result_domain,val=power_spectrum,
dtype=power_spectrum.dtype)
@classmethod @classmethod
def _calculate_power_spectrum(cls, field_val, pdomain, axes=None): def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
pindex = pdomain.pindex pindex = pdomain.pindex
# MR FIXME: how about iterating over slices, instead of replicating
# pindex? Would save memory and probably isn't slower.
if axes is not None: if axes is not None:
pindex = cls._shape_up_pindex( pindex = cls._shape_up_pindex(
pindex=pindex, pindex=pindex,
...@@ -711,8 +702,6 @@ class Field(object): ...@@ -711,8 +702,6 @@ class Field(object):
The weighted field. The weighted field.
""" """
new_field = self if inplace else self.copy_empty()
new_val = self.get_val(copy=False) new_val = self.get_val(copy=False)
spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain)) spaces = utilities.cast_axis_to_tuple(spaces, len(self.domain))
...@@ -725,9 +714,10 @@ class Field(object): ...@@ -725,9 +714,10 @@ class Field(object):
power=power, power=power,
axes=self.domain_axes[ind], axes=self.domain_axes[ind],
inplace=inplace) inplace=inplace)
# we need at most one copy, the rest can happen in place
inplace = True
new_field.set_val(new_val=new_val, copy=False) return Field(self.domain, new_val, self.dtype)
return new_field
def vdot(self, x=None, spaces=None, bare=False): def vdot(self, x=None, spaces=None, bare=False):
""" Computes the volume-factor-aware dot product of 'self' with x. """ Computes the volume-factor-aware dot product of 'self' with x.
...@@ -794,7 +784,6 @@ class Field(object): ...@@ -794,7 +784,6 @@ class Field(object):
The complex conjugated field. The complex conjugated field.
""" """
if inplace: if inplace:
self.imag*=-1 self.imag*=-1
return self return self
...@@ -837,10 +826,7 @@ class Field(object): ...@@ -837,10 +826,7 @@ class Field(object):
for i in range(len(self.domain)) for i in range(len(self.domain))
if i not in spaces) if i not in spaces)
return_field = Field(domain=return_domain, return Field(domain=return_domain, val=data, copy=False)
val=data,
copy=False)
return return_field
def sum(self, spaces=None): def sum(self, spaces=None):
return self._contraction_helper('sum', spaces) return self._contraction_helper('sum', spaces)
......
...@@ -224,9 +224,6 @@ class LinearOperator(with_metaclass( ...@@ -224,9 +224,6 @@ class LinearOperator(with_metaclass(
If the operator has an `inverse` then the inverse adjoint is identical If the operator has an `inverse` then the inverse adjoint is identical
to the adjoint inverse. We provide both names for convenience. to the adjoint inverse. We provide both names for convenience.
See Also
--------
""" """
spaces = self._check_input_compatibility(x, spaces) spaces = self._check_input_compatibility(x, spaces)
...@@ -261,8 +258,7 @@ class LinearOperator(with_metaclass( ...@@ -261,8 +258,7 @@ class LinearOperator(with_metaclass(
def _check_input_compatibility(self, x, spaces, inverse=False): def _check_input_compatibility(self, x, spaces, inverse=False):
if not isinstance(x, Field): if not isinstance(x, Field):
raise ValueError( raise ValueError("supplied object is not a `Field`.")
"supplied object is not a `Field`.")
if spaces is None and self.default_spaces is not None: if spaces is None and self.default_spaces is not None:
if not inverse: if not inverse:
...@@ -281,10 +277,7 @@ class LinearOperator(with_metaclass( ...@@ -281,10 +277,7 @@ class LinearOperator(with_metaclass(
# 2. Case: # 2. Case:
# The domains of self and x match completely. # The domains of self and x match completely.
if not inverse: self_domain = self.target if inverse else self.domain
self_domain = self.domain
else:
self_domain = self.target
if spaces is None: if spaces is None:
if self_domain != x.domain: if self_domain != x.domain:
......
...@@ -50,28 +50,22 @@ class ResponseOperator(LinearOperator): ...@@ -50,28 +50,22 @@ class ResponseOperator(LinearOperator):
default_spaces=None): default_spaces=None):
super(ResponseOperator, self).__init__(default_spaces) super(ResponseOperator, self).__init__(default_spaces)
self._domain = self._parse_domain(domain)
kernel_smoothing = len(self._domain)*[None]
kernel_exposure = len(self._domain)*[None]
if len(sigma) != len(exposure): if len(sigma) != len(exposure):
raise ValueError("Length of smoothing kernel and length of" raise ValueError("Length of smoothing kernel and length of"
"exposure do not match") "exposure do not match")
nsigma = len(sigma)
for ii in range(len(kernel_smoothing)): self._domain = self._parse_domain(domain)
kernel_smoothing[ii] = FFTSmoothingOperator(self._domain[ii],
sigma=sigma[ii]) kernel_smoothing = [FFTSmoothingOperator(self._domain[x], sigma[x])
kernel_exposure[ii] = DiagonalOperator(self._domain[ii], for x in range(nsigma)]
diagonal=exposure[ii]) kernel_exposure = [DiagonalOperator(self._domain[x],
diagonal=exposure[x]) for x in range(nsigma)]
self._composed_kernel = ComposedOperator(kernel_smoothing) self._composed_kernel = ComposedOperator(kernel_smoothing)
self._composed_exposure = ComposedOperator(kernel_exposure) self._composed_exposure = ComposedOperator(kernel_exposure)
target_list = [] target_list = [FieldArray(x.shape) for x in self.domain]
for space in self.domain:
target_list += [FieldArray(space.shape)]
self._target = self._parse_domain(target_list) self._target = self._parse_domain(target_list)