Commit d43638c1 authored by Martin Reinecke's avatar Martin Reinecke

introduce low-level support library; tweaks

parent 8d294b9b
Pipeline #18438 passed with stage
in 6 minutes and 24 seconds
......@@ -93,10 +93,14 @@ class DomainTuple(object):
def __eq__(self, x):
if not isinstance(x, DomainTuple):
x = DomainTuple(x)
x = DomainTuple.make(x)
if self is x:
return True
return self._dom == x._dom
def __ne__(self, x):
if not isinstance(x, DomainTuple):
x = DomainTuple(x)
x = DomainTuple.make(x)
if self is x:
return False
return self._dom != x._dom
......@@ -97,7 +97,7 @@ class Field(object):
if isinstance(val, Field):
return val.domain
if np.isscalar(val):
return DomainTuple(()) # empty domain tuple
return DomainTuple.make(()) # empty domain tuple
raise TypeError("could not infer domain from value")
return DomainTuple.make(domain)
......
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2017 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
import numpy as np
__all__ = ["hartley", "general_axpy"]
special_hartley = False
special_fill_array = False
use_numba = False
if special_hartley or special_fill_array:
import hartley as extmod
if not special_hartley:
def _fill_upper_half(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
_fill_upper_half(tmp[dim1], res[dim1], axes2)
def _fill_array(tmp, res, axes):
if axes is None:
axes = range(a.ndim)
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
np.add(tmp.real, tmp.imag, out=res[slice1])
_fill_upper_half(tmp, res, axes)
return res
def hartley(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes does not match array shape")
if issubclass(a.dtype.type, np.complexfloating):
raise TypeError("Hartley tansform requires real-valued arrays.")
if special_hartley:
return extmod.hartley(a, np.empty_like(a), axes)
else:
from pyfftw.interfaces.numpy_fft import rfftn
tmp = rfftn(a, axes=axes)
if special_fill_array:
return extmod.fill_array(tmp, np.empty_like(a), axes)
else:
return _fill_array(tmp, np.empty_like(a), axes)
if use_numba:
from numba import complex128 as ncplx, float64 as nflt, vectorize as nvct
@nvct([nflt(nflt,nflt,nflt), ncplx(nflt,ncplx,ncplx)])
def _general_axpy(a, x, y):
return a*x + y
def general_axpy(a, x, y, out):
if x.domain != y.domain or x.domain != out.domain:
raise ValueError ("Incompatible domains")
return _general_axpy(a, x.val, y.val, out.val)
else:
def general_axpy(a,x,y,out):
if x.domain != y.domain or x.domain != out.domain:
raise ValueError ("Incompatible domains")
x = x.val
y = y.val
out = out.val
if out is x:
if a != 1.:
out*=a
out += y
elif out is y:
if a!=1.:
out += a*x
else:
out += x
else:
out[()] = y
if a != 1.:
out += a*x
else:
out += x
return out
......@@ -19,7 +19,8 @@
from __future__ import division
from .minimizer import Minimizer
import numpy as np
from .. import Field
from ..low_level_library import general_axpy
class ConjugateGradient(Minimizer):
""" Implementation of the Conjugate Gradient scheme.
......@@ -73,7 +74,7 @@ class ConjugateGradient(Minimizer):
return energy, status
norm_b = energy.norm_b
r = -energy.gradient
r = energy.gradient
if preconditioner is not None:
d = preconditioner(r)
else:
......@@ -82,6 +83,7 @@ class ConjugateGradient(Minimizer):
if previous_gamma == 0:
return energy, controller.CONVERGED
tpos = Field(d.domain,dtype=d.dtype) # temporary buffer
while True:
q = energy.curvature(d)
ddotq = d.vdot(q).real
......@@ -92,14 +94,10 @@ class ConjugateGradient(Minimizer):
if alpha < 0:
return energy, controller.ERROR
# MR: BLAS candidate
q *= alpha
r -= q
general_axpy(-alpha, q, r, out=r)
# MR: BLAS candidate
tpos = d*alpha
tpos += energy.position
energy = energy.at_with_grad(tpos, -r)
general_axpy(-alpha, d, energy.position, out=tpos)
energy = energy.at_with_grad(tpos, r)
if preconditioner is not None:
s = preconditioner(r)
......@@ -119,8 +117,6 @@ class ConjugateGradient(Minimizer):
if status != controller.CONTINUE:
return energy, status
# MR: BLAS candidate
d *= max(0, gamma/previous_gamma)
d += s
general_axpy(max(0, gamma/previous_gamma), d, s, out=d)
previous_gamma = gamma
......@@ -19,7 +19,7 @@
from __future__ import division
import numpy as np
from ... import nifty_utilities as utilities
from ...low_level_library import hartley
class Transformation(object):
def __init__(self, domain, codomain):
......@@ -43,47 +43,6 @@ class RGRGTransformation(Transformation):
def unitary(self):
return True
@staticmethod
def _fill_upper_half(tmp, res, axes):
lastaxis = axes[-1]
nlast = res.shape[lastaxis]
ntmplast = tmp.shape[lastaxis]
nrem = nlast - ntmplast
slice1 = [slice(None)]*lastaxis + [slice(ntmplast, None)]
slice2 = [slice(None)]*lastaxis + [slice(nrem, 0, -1)]
for i in axes[:-1]:
slice1[i] = slice(1, None)
slice2[i] = slice(None, 0, -1)
np.subtract(tmp[slice2].real, tmp[slice2].imag, out=res[slice1])
for i, ax in enumerate(axes[:-1]):
dim1 = [slice(None)]*ax + [slice(0, 1)]
axes2 = axes[:i] + axes[i+1:]
RGRGTransformation._fill_upper_half(tmp[dim1], res[dim1], axes2)
@staticmethod
def _fill_array(tmp, res, axes):
if axes is None:
axes = range(a.ndim)
lastaxis = axes[-1]
ntmplast = tmp.shape[lastaxis]
slice1 = [slice(None)]*lastaxis + [slice(0, ntmplast)]
np.add(tmp.real, tmp.imag, out=res[slice1])
RGRGTransformation._fill_upper_half(tmp, res, axes)
return res
@staticmethod
def _hartley(a, axes=None):
# Check if the axes provided are valid given the shape
if axes is not None and \
not all(axis < len(a.shape) for axis in axes):
raise ValueError("Provided axes does not match array shape")
from pyfftw.interfaces.numpy_fft import rfftn
if issubclass(a.dtype.type, np.complexfloating):
raise TypeError("Hartley tansform requires real-valued arrays.")
tmp = rfftn(a, axes=axes)
return RGRGTransformation._fill_array(tmp, np.empty_like(a), axes)
def transform(self, val, axes=None):
"""
RG -> RG transform method.
......@@ -109,10 +68,10 @@ class RGRGTransformation(Transformation):
# Perform the transformation
if issubclass(val.dtype.type, np.complexfloating):
Tval = self._hartley(val.real, axes) \
+ 1j*self._hartley(val.imag, axes)
Tval = hartley(val.real, axes) \
+ 1j*hartley(val.imag, axes)
else:
Tval = self._hartley(val, axes)
Tval = hartley(val, axes)
return Tval, fct
......
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