Commit f23bb2b9 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'master' into index_games

parents 619ec79f 0c8a2610
Pipeline #14177 passed with stage
in 5 minutes and 7 seconds
......@@ -24,7 +24,7 @@ class WienerFilterEnergy(Energy):
@property
def value(self):
D_inv_x = self.D_inverse_x()
H = 0.5 * D_inv_x.dot(self.position) - self.j.dot(self.position)
H = 0.5 * D_inv_x.vdot(self.position) - self.j.dot(self.position)
return H.real
@property
......
......@@ -109,7 +109,7 @@ class LineEnergy(Energy):
@property
def gradient(self):
return self.energy.gradient.dot(self.line_direction)
return self.energy.gradient.vdot(self.line_direction)
@property
def curvature(self):
......
......@@ -170,6 +170,8 @@ class Field(Loggable, Versionable, object):
else:
dtype = np.dtype(dtype)
dtype = np.result_type(dtype, np.float)
return dtype
def _parse_distribution_strategy(self, distribution_strategy, val):
......@@ -1038,7 +1040,7 @@ class Field(Loggable, Versionable, object):
new_field.set_val(new_val=new_val, copy=False)
return new_field
def dot(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.
Parameters
......
......@@ -70,7 +70,7 @@ class ConjugateGradient(Loggable, object):
References
----------
Thomas V. Mikosch et al., "Numerical Optimization", Second Edition,
Jorge Nocedal & Stephen Wright, "Numerical Optimization", Second Edition,
2006, Springer-Verlag New York
"""
......@@ -121,12 +121,12 @@ class ConjugateGradient(Loggable, object):
r = b - A(x0)
d = self.preconditioner(r)
previous_gamma = r.dot(d)
previous_gamma = r.vdot(d)
if previous_gamma == 0:
self.logger.info("The starting guess is already perfect solution "
"for the inverse problem.")
return x0, self.convergence_level+1
norm_b = np.sqrt(b.dot(b))
norm_b = np.sqrt(b.vdot(b))
x = x0
convergence = 0
iteration_number = 1
......@@ -137,7 +137,7 @@ class ConjugateGradient(Loggable, object):
self.callback(x, iteration_number)
q = A(d)
alpha = previous_gamma/d.dot(q)
alpha = previous_gamma/d.vdot(q)
if not np.isfinite(alpha):
self.logger.error("Alpha became infinite! Stopping.")
......@@ -158,7 +158,7 @@ class ConjugateGradient(Loggable, object):
r -= q * alpha
s = self.preconditioner(r)
gamma = r.dot(s)
gamma = r.vdot(s)
if gamma.real < 0:
self.logger.warn("Positive definitness of preconditioner "
......
......@@ -137,7 +137,7 @@ class DescentMinimizer(Loggable, object):
# compute the the gradient for the current location
gradient = energy.gradient
gradient_norm = gradient.dot(gradient)
gradient_norm = gradient.vdot(gradient)
# check if position is at a flat point
if gradient_norm == 0:
......
......@@ -261,7 +261,7 @@ class InformationStore(object):
"""
key = tuple(sorted((i, j)))
if key not in self._ss_store:
self._ss_store[key] = self.s[i].dot(self.s[j])
self._ss_store[key] = self.s[i].vdot(self.s[j])
return self._ss_store[key]
def sy_store(self, i, j):
......@@ -284,7 +284,7 @@ class InformationStore(object):
"""
key = (i, j)
if key not in self._sy_store:
self._sy_store[key] = self.s[i].dot(self.y[j])
self._sy_store[key] = self.s[i].vdot(self.y[j])
return self._sy_store[key]
def yy_store(self, i, j):
......@@ -307,7 +307,7 @@ class InformationStore(object):
"""
key = tuple(sorted((i, j)))
if key not in self._yy_store:
self._yy_store[key] = self.y[i].dot(self.y[j])
self._yy_store[key] = self.y[i].vdot(self.y[j])
return self._yy_store[key]
def sgrad_store(self, i):
......@@ -319,7 +319,7 @@ class InformationStore(object):
Scalar product.
"""
return self.s[i].dot(self.last_gradient)
return self.s[i].vdot(self.last_gradient)
def ygrad_store(self, i):
"""Returns scalar product between y_i and gradient on initial position.
......@@ -330,7 +330,7 @@ class InformationStore(object):
Scalar product.
"""
return self.y[i].dot(self.last_gradient)
return self.y[i].vdot(self.last_gradient)
def gradgrad_store(self):
"""Returns scalar product of gradient on initial position with itself.
......@@ -341,7 +341,7 @@ class InformationStore(object):
Scalar product.
"""
return self.last_gradient.dot(self.last_gradient)
return self.last_gradient.vdot(self.last_gradient)
def add_new_point(self, x, gradient):
"""Updates the s list and y list.
......
......@@ -51,7 +51,7 @@ class ComposedOperator(LinearOperator):
Notes
-----
Very usefull in case one has to transform a Field living over a product
Very useful in case one has to transform a Field living over a product
space (see example below).
Examples
......@@ -96,6 +96,8 @@ class ComposedOperator(LinearOperator):
is not easily forecasteable what the output of an operator-call
will look like.
"""
if spaces is None:
spaces = self.default_spaces
return spaces
# ---Mandatory properties and methods---
......
......@@ -169,89 +169,8 @@ class DiagonalOperator(EndomorphicOperator):
out : Field
The inverse of the diagonal of the Operator.
"""
return 1./self.diagonal(bare=bare, copy=False)
def trace(self, bare=False):
""" Returns the trace the operator.
Parameters
----------
bare : boolean
Whether the returned Field values should be bare or not.
Returns
-------
out : scalar
The trace of the Operator.
"""
return self.diagonal(bare=bare, copy=False).sum()
def inverse_trace(self, bare=False):
""" Returns the inverse-trace of the operator.
Parameters
----------
bare : boolean
Whether the returned Field values should be bare or not.
Returns
-------
out : scalar
The inverse of the trace of the Operator.
"""
return self.inverse_diagonal(bare=bare).sum()
def trace_log(self):
""" Returns the trave-log of the operator.
Returns
-------
out : scalar
the trace of the logarithm of the Operator.
"""
log_diagonal = nifty_log(self.diagonal(copy=False))
return log_diagonal.sum()
def determinant(self):
""" Returns the determinant of the operator.
Returns
-------
out : scalar
out : scalar
the determinant of the Operator
"""
return self.diagonal(copy=False).val.prod()
def inverse_determinant(self):
""" Returns the inverse-determinant of the operator.
Returns
-------
out : scalar
the inverse-determinant of the Operator
"""
return 1/self.determinant()
def log_determinant(self):
""" Returns the log-eterminant of the operator.
Returns
-------
out : scalar
the log-determinant of the Operator
"""
return np.log(self.determinant())
return 1./self.diagonal(bare=bare, copy=False)
# ---Mandatory properties and methods---
......
......@@ -84,7 +84,7 @@ class ResponseOperator(LinearOperator):
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"
"exposure do not match")
......
......@@ -53,7 +53,7 @@ class DirectSmoothingOperator(SmoothingOperator):
wgt = []
expfac = 1. / (2. * sigma*sigma)
for i in range(x.size):
if nval[i]>0:
if nval[i] > 0:
t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
t = np.exp(-t*t*expfac)
t *= 1./np.sum(t)
......
......@@ -25,10 +25,6 @@ class Figure2D(FigureFromPlot):
xaxis = False if (xaxis is None) else xaxis
yaxis = False if (yaxis is None) else yaxis
else:
width = None
height = None
super(Figure2D, self).__init__(plots, title, width, height)
self.xaxis = xaxis
self.yaxis = yaxis
......
......@@ -2,4 +2,5 @@
from healpix_plotter import HealpixPlotter
from gl_plotter import GLPlotter
from power_plotter import PowerPlotter
from rg1d_plotter import RG1DPlotter
from rg2d_plotter import RG2DPlotter
......@@ -87,8 +87,8 @@ class PlotterBase(Loggable, object):
for slice_list in utilities.get_slice_list(data_list[0].shape, axes):
plots_list += \
[[self.plot.at(self._parse_data(current_data,
field,
spaces))
field,
spaces))
for (current_data, field) in zip(data_list, fields)]]
figures = [self.figure.at(plots) for plots in plots_list]
......
# -*- coding: utf-8 -*-
import numpy as np
from nifty.spaces import RGSpace
from nifty.plotting.descriptors import Axis
from nifty.plotting.figures import Figure2D
from nifty.plotting.plots import Cartesian2D
from .plotter_base import PlotterBase
class RG1DPlotter(PlotterBase):
def __init__(self, interactive=False, path='.', title="",
line=None, marker=None):
super(RG1DPlotter, self).__init__(interactive, path, title)
self.line = line
self.marker = marker
@property
def domain_classes(self):
return (RGSpace, )
def _initialize_plot(self):
return Cartesian2D(data=None)
def _initialize_figure(self):
xaxis = Axis()
yaxis = Axis()
return Figure2D(plots=None, xaxis=xaxis, yaxis=yaxis)
def _parse_data(self, data, field, spaces):
y_data = data
rgspace = field.domain[spaces[0]]
xy_data = np.empty((2, y_data.shape[0]))
xy_data[1] = y_data
num = rgspace.shape[0]
length = rgspace.distances[0]*num
xy_data[0] = np.linspace(start=0,
stop=length,
num=num,
endpoint=False)
if rgspace.zerocenter[0]:
xy_data[0] -= np.floor(length/2)
return xy_data
......@@ -30,7 +30,7 @@ class TraceProberMixin(object):
super(TraceProberMixin, self).reset()
def finish_probe(self, probe, pre_result):
result = probe[1].dot(pre_result, bare=True)
result = probe[1].vdot(pre_result, bare=True)
self.__sum_of_probings += result
if self.compute_variance:
self.__sum_of_squares += result.conjugate() * result
......
......@@ -184,7 +184,9 @@ class LMSpace(Space):
return res
def get_fft_smoothing_kernel_function(self, sigma):
# FIXME why x(x+1) ? add reference to paper!
# cf. "All-sky convolution for polarimetry experiments"
# by Challinor et al.
# http://arxiv.org/abs/astro-ph/0008228
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma*sigma)
# ---Added properties and methods---
......
......@@ -314,7 +314,7 @@ class RGSpace(Space):
return 0.5*(tmp[:-1]+tmp[1:])
def get_fft_smoothing_kernel_function(self, sigma):
return lambda x: np.exp(-0.5 * np.pi*np.pi * x*x * sigma*sigma)
return lambda x: np.exp(-2. * np.pi*np.pi * x*x * sigma*sigma)
# ---Added properties and methods---
......
......@@ -40,8 +40,8 @@ class ComposedOperator_Tests(unittest.TestCase):
rand1 = Field.from_random('normal', domain=(space1,space2))
rand2 = Field.from_random('normal', domain=(space1,space2))
tt1 = rand2.dot(op.times(rand1))
tt2 = rand1.dot(op.adjoint_times(rand2))
tt1 = rand2.vdot(op.times(rand1))
tt2 = rand1.vdot(op.adjoint_times(rand2))
assert_approx_equal(tt1, tt2)
@expand(product(spaces, spaces))
......
......@@ -33,8 +33,8 @@ class DiagonalOperator_Tests(unittest.TestCase):
rand2 = Field.from_random('normal', domain=space)
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
tt1 = rand1.dot(D.times(rand2))
tt2 = rand2.dot(D.times(rand1))
tt1 = rand1.vdot(D.times(rand2))
tt2 = rand2.vdot(D.times(rand1))
assert_approx_equal(tt1, tt2)
@expand(product(spaces, [True, False], [True, False]))
......@@ -89,48 +89,5 @@ class DiagonalOperator_Tests(unittest.TestCase):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
diag_op = D.inverse_diagonal()
assert_allclose(1./diag.val.get_full_data(), diag_op.val.get_full_data())
@expand(product(spaces, [True, False]))
def test_trace(self, space, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
trace_op = D.trace()
assert_allclose(trace_op, np.sum(diag.val.get_full_data()))
@expand(product(spaces, [True, False]))
def test_inverse_trace(self, space, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
trace_op = D.inverse_trace()
assert_allclose(trace_op, np.sum(1./diag.val.get_full_data()))
@expand(product(spaces, [True, False]))
#MR FIXME: what if any diagonal element <=0?
def test_trace_log(self, space, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
trace_log = D.trace_log()
assert_allclose(trace_log, np.sum(np.log(diag.val.get_full_data())))
@expand(product(spaces, [True, False]))
def test_determinant(self, space, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, copy=copy)
det = D.determinant()
assert_allclose(det, np.prod(diag.val.get_full_data()))
@expand(product(spaces, [True, False], [True, False]))
def test_inverse_determinant(self, space, bare, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
inv_det = D.inverse_determinant()
assert_allclose(inv_det, 1./D.determinant())
@expand(product(spaces, [True, False], [True, False]))
#MR FIXME: what if determinant <=0?
def test_log_determinant(self, space, bare, copy):
diag = Field.from_random('normal', domain=space)
D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
log_det = D.log_determinant()
assert_allclose(log_det, np.log(D.determinant()))
assert_allclose(1./diag.val.get_full_data(),
diag_op.val.get_full_data())
......@@ -143,8 +143,8 @@ class FFTOperatorTests(unittest.TestCase):
inp = Field.from_random(domain=a, random_type='normal', std=1, mean=0,
dtype=tp)
out = fft.times(inp)
v1 = np.sqrt(out.dot(out))
v2 = np.sqrt(inp.dot(fft.adjoint_times(out)))
v1 = np.sqrt(out.vdot(out))
v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
assert_allclose(v1, v2, rtol=tol, atol=tol)
@expand(product([128, 256],
......@@ -159,6 +159,6 @@ class FFTOperatorTests(unittest.TestCase):
inp = Field.from_random(domain=a, random_type='normal', std=1, mean=0,
dtype=tp)
out = fft.times(inp)
v1 = np.sqrt(out.dot(out))
v2 = np.sqrt(inp.dot(fft.adjoint_times(out)))
v1 = np.sqrt(out.vdot(out))
v2 = np.sqrt(inp.vdot(fft.adjoint_times(out)))
assert_allclose(v1, v2, rtol=tol, atol=tol)
Supports Markdown
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