Skip to content

GitLab

Commit 77ee7dcc authored by Martin Reinecke's avatar Martin Reinecke
Browse files

simplify field class

parent 3c829eb0
Hide whitespace changes
Inline Side-by-side
Showing with 81 additions and 239 deletions
demos/log_normal_wiener_filter.py
View file @ 77ee7dcc
# -*- coding: utf-8 -*-
from nifty import *
import nifty as ift
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
......@@ -22,7 +21,7 @@ if __name__ == "__main__":
#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)
fft = FFTOperator(harmonic_space, target=signal_space)
power_space = PowerSpace(harmonic_space)
# Creating the mock signal |\label{code:wf_mock_signal}|
......@@ -46,24 +45,18 @@ if __name__ == "__main__":
# 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)
ctrl = ift.DefaultIterationController(verbose=False,tol_abs_gradnorm=1)
ctrl2 = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1, name="outer")
inverter = ift.ConjugateGradient(controller=ctrl, preconditioner=S.times)
energy = library.LogNormalWienerFilterEnergy(m0, data, R_harmonic, N, S, inverter=inverter)
minimizer1 = VL_BFGS(controller=ctrl2,max_history_length=20)
minimizer2 = RelaxedNewton(controller=ctrl2)
minimizer3 = SteepestDescent(controller=ctrl2)
# me1 = minimizer1(energy)
# me2 = minimizer2(energy)
# me3 = minimizer3(energy)
me3 = minimizer3(energy)
# m1 = fft(me1[0].position)
# m2 = fft(me2[0].position)
......
demos/paper_demos/cartesian_wiener_filter.py
View file @ 77ee7dcc
......@@ -26,7 +26,7 @@ if __name__ == "__main__":
fft_1 = ift.FFTOperator(harmonic_space_1, target=signal_space_1)
power_space_1 = ift.PowerSpace(harmonic_space_1)
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1)
mock_power_1 = ift.Field(power_space_1, val=power_spectrum_1(power_space_1.kindex))
......@@ -49,7 +49,7 @@ if __name__ == "__main__":
fft_2 = ift.FFTOperator(harmonic_space_2, target=signal_space_2)
power_space_2 = ift.PowerSpace(harmonic_space_2)
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2)
mock_power_2 = ift.Field(power_space_2, val=power_spectrum_2(power_space_2.kindex))
fft = ift.ComposedOperator((fft_1, fft_2))
......
demos/paper_demos/wiener_filter.py
View file @ 77ee7dcc
......@@ -26,7 +26,7 @@ if __name__ == "__main__":
# Creating the mock signal |\label{code:wf_mock_signal}|
S = ift.create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = ift.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))
# Setting up an exemplary response
......
demos/wiener_filter_via_curvature.py
View file @ 77ee7dcc
......@@ -41,7 +41,7 @@ if __name__ == "__main__":
# Creating the mock data
S = create_power_operator(harmonic_space, power_spectrum=power_spectrum)
mock_power = Field(power_space, val=power_spectrum)
mock_power = Field(power_space, val=power_spectrum(power_space.kindex))
np.random.seed(43)
mock_harmonic = mock_power.power_synthesize(real_signal=True)
mock_harmonic = mock_harmonic.real
......
demos/wiener_filter_via_hamiltonian.py
View file @ 77ee7dcc
from nifty import *
import plotly.offline as pl
import plotly.graph_objs as go
from nifty.library.wiener_filter import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
rank = comm.rank
import nifty as ift
import numpy as np
np.random.seed(42)
class AdjointFFTResponse(LinearOperator):
class AdjointFFTResponse(ift.LinearOperator):
def __init__(self, FFT, R, default_spaces=None):
super(AdjointFFTResponse, self).__init__(default_spaces)
self._domain = FFT.target
......@@ -39,43 +32,38 @@ class AdjointFFTResponse(LinearOperator):
if __name__ == "__main__":
distribution_strategy = 'not'
# Set up position space
s_space = RGSpace([128, 128])
s_space = ift.RGSpace([128, 128])
# s_space = HPSpace(32)
# Define harmonic transformation and associated harmonic space
fft = FFTOperator(s_space)
fft = ift.FFTOperator(s_space)
h_space = fft.target[0]
# Setting up power space
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
p_space = ift.PowerSpace(h_space)
# Choosing the prior correlation structure and defining
# correlation operator
p_spec = (lambda k: (42 / (k + 1) ** 3))
S = create_power_operator(h_space, power_spectrum=p_spec,
distribution_strategy=distribution_strategy)
S = ift.create_power_operator(h_space, power_spectrum=p_spec)
# Drawing a sample sh from the prior distribution in harmonic space
sp = Field(p_space, val=p_spec,
distribution_strategy=distribution_strategy)
sp = ift.Field(p_space, val=p_spec(p_space.kindex))
sh = sp.power_synthesize(real_signal=True)
ss = fft.adjoint_times(sh)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.05)
Instrument = DiagonalOperator(s_space, diagonal=1.)
Instrument = ift.DiagonalOperator(s_space, diagonal=1.)
# Instrument._diagonal.val[200:400, 200:400] = 0
# Adding a harmonic transformation to the instrument
R = AdjointFFTResponse(fft, Instrument)
signal_to_noise = 1.
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
n = Field.from_random(domain=s_space,
N = ift.DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
n = ift.Field.from_random(domain=s_space,
random_type='normal',
std=ss.std()/np.sqrt(signal_to_noise),
mean=0)
......@@ -86,44 +74,26 @@ if __name__ == "__main__":
# Choosing the minimization strategy
def convergence_measure(energy, iteration): # returns current energy
x = energy.value
print(x, iteration)
# minimizer = SteepestDescent(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure)
minimizer = RelaxedNewton(convergence_tolerance=0,
iteration_limit=1,
callback=convergence_measure)
#
# minimizer = VL_BFGS(convergence_tolerance=0,
# iteration_limit=50,
# callback=convergence_measure,
# max_history_length=3)
#
inverter = ConjugateGradient(convergence_level=3,
convergence_tolerance=1e-5,
preconditioner=None)
ctrl = ift.DefaultIterationController(verbose=True,tol_abs_gradnorm=0.1)
inverter = ift.ConjugateGradient(controller=ctrl, preconditioner=S.times)
# Setting starting position
m0 = Field(h_space, val=.0)
m0 = ift.Field(h_space, val=.0)
# Initializing the Wiener Filter energy
energy = WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S)
energy = ift.library.WienerFilterEnergy(position=m0, d=d, R=R, N=N, S=S,
inverter=inverter)
D0 = energy.curvature
# Solving the problem analytically
m0 = D0.inverse_times(j)
sample_variance = Field(sh.domain, val=0.)
sample_mean = Field(sh.domain, val=0.)
sample_variance = ift.Field(sh.domain, val=0.)
sample_mean = ift.Field(sh.domain, val=0.)
# sampling the uncertainty map
n_samples = 10
n_samples = 50
for i in range(n_samples):
sample = fft(sugar.generate_posterior_sample(0., D0))
sample = fft(ift.sugar.generate_posterior_sample(0., D0))
sample_variance += sample**2
sample_mean += sample
variance = (sample_variance - sample_mean**2)/n_samples
nifty/domain_object.py
View file @ 77ee7dcc
......@@ -167,57 +167,3 @@ class DomainObject(with_metaclass(
"""
raise NotImplementedError(
"There is no generic weight-method for DomainObject.")
def pre_cast(self, x, axes):
""" Casts input for Field.val before Field performs the cast.
Parameters
----------
x : {array-like, castable}
an array-like object or anything that can be cast to arrays.
axes : tuple of ints
Specifies the axes of x which correspond to this domain.
Returns
-------
{array-like, castable}
Processed input where casting that needs Space-specific knowledge
(for example location of pixels on the manifold) was performed.
See Also
--------
post_cast
Notes
-----
Usually returns x, except if a power spectrum is given to a
PowerSpace, where this spectrum is evaluated at the power indices.
"""
return x
def post_cast(self, x, axes):
""" Performs casting operations that are done after Field's cast.
Parameters
----------
x : {array-like, castable}
an array-like object or anything that can be cast to arrays.
axes : tuple of ints
Specifies the axes of x which correspond to this domain.
See Also
--------
pre_cast
Returns
-------
numpy.ndarray
Processed input where casting that needs Space-specific knowledge
(for example location of pixels on the manifold) was performed.
"""
return x
nifty/field.py
View file @ 77ee7dcc
......@@ -82,19 +82,24 @@ class Field(object):
def __init__(self, domain=None, val=None, dtype=None, copy=False):
self.domain = self._parse_domain(domain=domain, val=val)
self.domain_axes = self._get_axes_tuple(self.domain)
self.dtype = self._infer_dtype(dtype=dtype, val=val)
if val is None:
self._val = None
shape_tuple = tuple(sp.shape for sp in self.domain)
if len(shape_tuple)==0:
global_shape = ()
else:
self.set_val(new_val=val, copy=copy)
global_shape = reduce(lambda x, y: x + y, shape_tuple)
dtype = self._infer_dtype(dtype=dtype, val=val)
self._val = np.empty(global_shape,dtype=dtype)
self.set_val(new_val=val, copy=copy)
def _parse_domain(self, domain, val=None):
if domain is None:
if isinstance(val, Field):
domain = val.domain
else:
elif np.isscalar(val):
domain = ()
else:
raise TypeError("could not infer domain from value")
elif isinstance(domain, DomainObject):
domain = (domain,)
elif not isinstance(domain, tuple):
......@@ -528,10 +533,26 @@ class Field(object):
"""
new_val = self.cast(new_val)
if copy:
new_val = new_val.copy()
self._val = new_val
if new_val is None:
pass
elif isinstance(new_val, Field):
if self.domain!=new_val.domain:
raise ValueError("Domain mismatch")
if copy:
self._val[()] = new_val.val
else:
self._val = new_val.val
elif (np.isscalar(new_val)):
self._val[()]=new_val
elif isinstance(new_val, np.ndarray):
if copy:
self._val[()] = new_val
else:
if self.shape!=new_val.shape:
raise ValueError("Shape mismatch")
self._val = new_val
else:
raise TypeError("unknown source type")
return self
def get_val(self, copy=False):
......@@ -553,13 +574,7 @@ class Field(object):
"""
if self._val is None:
self.set_val(None)
if copy:
return self._val.copy()
else:
return self._val
return self._val.copy() if copy else self._val
@property
def val(self):
......@@ -575,12 +590,16 @@ class Field(object):
"""
return self.get_val(copy=False)
return self._val
@val.setter
def val(self, new_val):
self.set_val(new_val=new_val, copy=False)
@property
def dtype(self):
return self._val.dtype
@property
def shape(self):
""" Returns the total shape of the Field's data array.
......@@ -596,14 +615,7 @@ class Field(object):
dim
"""
if not hasattr(self, '_shape'):
shape_tuple = tuple(sp.shape for sp in self.domain)
try:
global_shape = reduce(lambda x, y: x + y, shape_tuple)
except TypeError:
global_shape = ()
self._shape = global_shape
return self._shape
return self._val.shape
@property
def dim(self):
......@@ -672,60 +684,6 @@ class Field(object):
# ---Special unary/binary operations---
def cast(self, x=None, dtype=None):
""" Transforms x to an object with the correct dtype and shape.
Parameters
----------
x : scalar, numpy.ndarray, Field, array_like
The input that shall be casted on a numpy.ndarray of the same shape
like the domain.
dtype : type
The datatype the output shall have. This can be used to override
the field's dtype.
Returns
-------
out : numpy.ndarray
The output object.
See Also
--------
_actual_cast
"""
if dtype is None:
dtype = self.dtype
else:
dtype = np.dtype(dtype)
casted_x = x
for ind, sp in enumerate(self.domain):
casted_x = sp.pre_cast(casted_x,
axes=self.domain_axes[ind])
casted_x = self._actual_cast(casted_x, dtype=dtype)
for ind, sp in enumerate(self.domain):
casted_x = sp.post_cast(casted_x,
axes=self.domain_axes[ind])
return casted_x
def _actual_cast(self, x, dtype=None):
if isinstance(x, Field):
x = x.get_val()
if dtype is None:
dtype = self.dtype
if x is not None:
if np.isscalar(x):
return np.full(self.shape,x, dtype=dtype)
return np.asarray(x, dtype=dtype).reshape(self.shape)
else:
return np.empty(self.shape, dtype=dtype)
def copy(self, domain=None, dtype=None):
""" Returns a full copy of the Field.
......
nifty/spaces/power_space/power_space.py
View file @ 77ee7dcc
......@@ -153,31 +153,6 @@ class PowerSpace(Space):
binbounds = harmonic_partner.get_natural_binbounds()
return np.searchsorted(binbounds, distance_array)
def pre_cast(self, x, axes):
""" Casts power spectrum functions to discretized power spectra.
This function takes an array or a function. If it is an array it does
nothing, otherwise it interpretes the function as power spectrum and
evaluates it at every k-mode.
Parameters
----------
x : {array-like, function array-like -> array-like}
power spectrum given either in discretized form or implicitly as a
function
axes : tuple of ints
Specifies the axes of x which correspond to this space. For
explicifying the power spectrum function, this is ignored.
Returns
-------
array-like
discretized power spectrum
"""
return x(self.kindex) if callable(x) else x
# ---Mandatory properties and methods---
def __repr__(self):
......
nifty/sugar.py
View file @ 77ee7dcc
......@@ -40,9 +40,9 @@ def create_power_operator(domain, power_spectrum, dtype=None):
----------
domain : DomainObject
Domain over which the power operator shall live.
power_spectrum : (array-like, method)
An array-like object, or a method that implements the square root
of a power spectrum as a function of k.
power_spectrum : callable
A method that implements the square root of a power spectrum as a
function of k.
dtype : type *optional*
dtype that the field holding the power spectrum shall use
(default : None).
......@@ -54,12 +54,12 @@ def create_power_operator(domain, power_spectrum, dtype=None):
"""
if isinstance(power_spectrum, Field):
power_domain = power_spectrum.domain
else:
power_domain = PowerSpace(domain)
if not callable(power_spectrum):
raise TypeError("power_spectrum must be callable")
power_domain = PowerSpace(domain)
fp = Field(power_domain, val=power_spectrum, dtype=dtype)
fp = Field(power_domain,
val=power_spectrum(power_domain.kindex), dtype=dtype)
f = fp.power_synthesize(mean=1, std=0, real_signal=False)
if not issubclass(fp.dtype.type, np.complexfloating):
......
test/test_field.py
View file @ 77ee7dcc
......@@ -91,11 +91,11 @@ class Test_Functionality(unittest.TestCase):
p1 = PowerSpace(space1)
spec1 = lambda k: 42/(1+k)**2
fp1 = Field(p1, val=spec1)
fp1 = Field(p1, val=spec1(p1.kindex))
p2 = PowerSpace(space2)
spec2 = lambda k: 42/(1+k)**3
fp2 = Field(p2, val=spec2)
fp2 = Field(p2, val=spec2(p2.kindex))
outer = np.outer(fp1.val, fp2.val)
fp = Field((p1, p2), val=outer)
......
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