Commit 91f06b35 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge master

parents 58378e58 2a7ee176
Pipeline #14424 passed with stage
in 5 minutes and 20 seconds
......@@ -115,7 +115,7 @@ if __name__ == "__main__":
# Solving the problem analytically
m0 = D0.inverse_times(j)
sample_variance = Field(sh.power_analyze(logarithmic=False).domain,val=0. + 0j)
sample_variance = Field(sh.domain,val=0. + 0j)
sample_mean = Field(sh.domain,val=0. + 0j)
# sampling the uncertainty map
......
......@@ -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
......@@ -1098,7 +1100,7 @@ class Field(Loggable, Versionable, object):
The Lq-norm of the field values.
"""
return np.sqrt(np.abs(self.dot(x=self)))
return np.sqrt(np.abs(self.vdot(x=self)))
def conjugate(self, inplace=False):
""" Retruns the complex conjugate of the field.
......
......@@ -75,9 +75,9 @@ class CriticalPowerEnergy(Energy):
@property
def value(self):
energy = exp(-self.position).dot(self.q + self.w / 2., bare= True)
energy += self.position.dot(self.alpha - 1. + self.rho / 2., bare=True)
energy += 0.5 * self.position.dot(self.T(self.position))
energy = exp(-self.position).vdot(self.q + self.w / 2., bare= True)
energy += self.position.vdot(self.alpha - 1. + self.rho / 2., bare=True)
energy += 0.5 * self.position.vdot(self.T(self.position))
return energy.real
@property
......
......@@ -36,8 +36,8 @@ class WienerFilterEnergy(Energy):
@property
def value(self):
residual = self._residual()
energy = 0.5 * self.position.dot(self.S.inverse_times(self.position))
energy += 0.5 * (residual).dot(self.N.inverse_times(residual))
energy = 0.5 * self.position.vdot(self.S.inverse_times(self.position))
energy += 0.5 * (residual).vdot(self.N.inverse_times(residual))
return energy.real
@property
......
......@@ -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,13 +121,13 @@ class ConjugateGradient(Loggable, object):
r = b - A(x0)
d = self.preconditioner(r)
previous_gamma = r.dot(d)
previous_gamma = (r.vdot(d)).real
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))
x = x0
norm_b = np.sqrt((b.vdot(b)).real)
x = x0.copy()
convergence = 0
iteration_number = 1
self.logger.info("Starting conjugate gradient.")
......@@ -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).real
if not np.isfinite(alpha):
self.logger.error("Alpha became infinite! Stopping.")
......@@ -146,7 +146,7 @@ class ConjugateGradient(Loggable, object):
x += d * alpha
reset = False
if alpha.real < 0:
if alpha < 0:
self.logger.warn("Positive definiteness of A violated!")
reset = True
if self.reset_count is not None:
......@@ -158,9 +158,9 @@ class ConjugateGradient(Loggable, object):
r -= q * alpha
s = self.preconditioner(r)
gamma = r.dot(s)
gamma = r.vdot(s).real
if gamma.real < 0:
if gamma < 0:
self.logger.warn("Positive definitness of preconditioner "
"violated!")
......@@ -170,10 +170,7 @@ class ConjugateGradient(Loggable, object):
self.logger.debug("Iteration : %08u alpha = %3.1E "
"beta = %3.1E delta = %3.1E" %
(iteration_number,
np.real(alpha),
np.real(beta),
np.real(delta)))
(iteration_number, alpha, beta, delta))
if gamma == 0:
convergence = self.convergence_level+1
......
......@@ -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---
......
......@@ -85,7 +85,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
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