Commit a50ec021 authored by Reimar H Leike's avatar Reimar H Leike
Browse files

Merge branch 'performance_pa' into 'NIFTy_6'

Performance pa

See merge request !417
parents 25808c2e f2975f76
Pipeline #70486 passed with stages
in 16 minutes and 32 seconds
......@@ -123,6 +123,44 @@ def _domain_check(op):
'be instances of either DomainTuple or MultiDomain.')
def _performance_check(op, pos, raise_on_fail):
class CountingOp(LinearOperator):
def __init__(self, domain):
from .sugar import makeDomain
self._domain = self._target = makeDomain(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES
self._count = 0
def apply(self, x, mode):
self._count += 1
return x
@property
def count(self):
return self._count
for wm in [False, True]:
cop = CountingOp(op.domain)
op = op @ cop
op(pos)
cond = [cop.count != 1]
lin = op(2*Linearization.make_var(pos, wm))
cond.append(cop.count != 2)
lin.jac(pos)
cond.append(cop.count != 3)
lin.jac.adjoint(lin.val)
cond.append(cop.count != 4)
if wm and op.target is DomainTuple.scalar_domain():
lin.metric(pos)
cond.append(cop.count != 6)
if any(cond):
s = 'The operator has a performance problem (want_metric={}).'.format(wm)
from .logger import logger
logger.error(s)
logger.info(cond)
if raise_on_fail:
raise RuntimeError(s)
def consistency_check(op, domain_dtype=np.float64, target_dtype=np.float64,
atol=0, rtol=1e-7, only_r_linear=False):
"""
......@@ -199,7 +237,15 @@ def _get_acceptable_location(op, loc, lin):
return loc2, lin2
def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100):
def _linearization_value_consistency(op, loc):
for wm in [False, True]:
lin = Linearization.make_var(loc, wm)
fld0 = op(loc)
fld1 = op(lin).val
assert_allclose(fld0, fld1, 0, 1e-7)
def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100, perf_check=True):
"""
Checks the Jacobian of an operator against its finite difference
approximation.
......@@ -216,9 +262,13 @@ def check_jacobian_consistency(op, loc, tol=1e-8, ntries=100):
as op. The location at which the gradient is checked
tol : float
Tolerance for the check.
perf_check : Boolean
Do performance check. May be disabled for very unimportant operators.
"""
_domain_check(op)
_actual_domain_check_nonlinear(op, loc)
_performance_check(op, loc, bool(perf_check))
_linearization_value_consistency(op, loc)
for _ in range(ntries):
lin = op(Linearization.make_var(loc))
loc2, lin2 = _get_acceptable_location(op, loc, lin)
......
......@@ -35,9 +35,9 @@ from ..operators.endomorphic_operator import EndomorphicOperator
from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.linear_operator import LinearOperator
from ..operators.operator import Operator
from ..operators.simple_linear_operators import ducktape
from ..operators.simple_linear_operators import FieldAdapter, ducktape
from ..probing import StatCalculator
from ..sugar import makeField, full, makeDomain
from ..sugar import full, makeDomain, makeField, makeOp
def _reshaper(x, N):
......@@ -207,7 +207,7 @@ class _TwoLogIntegrations(LinearOperator):
class _Normalization(Operator):
def __init__(self, domain, space=0):
self._domain = self._target = makeDomain(domain)
self._domain = self._target = DomainTuple.make(domain)
assert isinstance(self._domain[space], PowerSpace)
hspace = list(self._domain)
hspace[space] = hspace[space].harmonic_partner
......@@ -218,16 +218,17 @@ class _Normalization(Operator):
mode_multiplicity = pd.adjoint(full(pd.target, 1.)).val_rw()
zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
mode_multiplicity[zero_mode] = 0
self._mode_multiplicity = makeField(self._domain, mode_multiplicity)
self._mode_multiplicity = makeOp(makeField(self._domain, mode_multiplicity))
self._specsum = _SpecialSum(self._domain, space)
def apply(self, x):
self._check_input(x)
amp = x.exp()
spec = (2*x).exp()
fa = FieldAdapter(self._domain, 'foo')
amp = fa.exp()
spec = (2*fa).exp()
# FIXME This normalizes also the zeromode which is supposed to be left
# untouched by this operator
return self._specsum(self._mode_multiplicity*spec)**(-0.5)*amp
return (self._specsum(self._mode_multiplicity(spec))**(-0.5)*amp)(fa.adjoint(x))
class _SpecialSum(EndomorphicOperator):
......
......@@ -343,7 +343,8 @@ for op in ["__iadd__", "__isub__", "__imul__", "__idiv__",
setattr(MultiField, op, func(op))
for f in ["sqrt", "exp", "log", "log1p", "expm1", "tanh", "one_over"]:
for f in ["sqrt", "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "tanh",
"absolute", "sinc", "sign", "log10", "log1p", "expm1"]:
def func(f):
def func2(self):
fu = getattr(Field, f)
......
......@@ -15,9 +15,12 @@
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..field import Field
from ..multi_field import MultiField
from .operator import Operator
from ..sugar import makeDomain
class Adder(Operator):
......@@ -25,18 +28,22 @@ class Adder(Operator):
Parameters
----------
field : Field or MultiField
a : Field or MultiField or Scalar
The field by which the input is shifted.
"""
def __init__(self, field, neg=False):
if not isinstance(field, (Field, MultiField)):
def __init__(self, a, neg=False, domain=None):
self._a = a
if isinstance(a, (Field, MultiField)):
dom = a.domain
elif np.isscalar(a):
dom = makeDomain(domain)
else:
raise TypeError
self._field = field
self._domain = self._target = field.domain
self._domain = self._target = dom
self._neg = bool(neg)
def apply(self, x):
self._check_input(x)
if self._neg:
return x - self._field
return x + self._field
return x - self._a
return x + self._a
......@@ -19,17 +19,17 @@ import numpy as np
from .. import utilities
from ..domain_tuple import DomainTuple
from ..multi_domain import MultiDomain
from ..field import Field
from ..multi_field import MultiField
from ..linearization import Linearization
from ..sugar import makeDomain, makeOp, full
from ..multi_domain import MultiDomain
from ..multi_field import MultiField
from ..sugar import makeDomain, makeOp
from .linear_operator import LinearOperator
from .operator import Operator
from .sampling_enabler import SamplingEnabler
from .sandwich_operator import SandwichOperator
from .scaling_operator import ScalingOperator
from .simple_linear_operators import VdotOperator, FieldAdapter
from .simple_linear_operators import FieldAdapter, VdotOperator
class EnergyOperator(Operator):
......@@ -130,17 +130,13 @@ class VariableCovarianceGaussianEnergy(EnergyOperator):
def apply(self, x):
self._check_input(x)
from .contraction_operator import ContractionOperator
lin = isinstance(x, Linearization)
r = FieldAdapter(self._domain[self._r], self._r)
icov = FieldAdapter(self._domain[self._icov], self._icov)
res0 = r.vdot(r*icov).real
res1 = icov.log().sum()
res = 0.5*(res0-res1)
res = res(x)
if not lin:
return Field.scalar(res)
if not x.want_metric:
res = (res0-res1).scale(0.5)(x)
if not lin or not x.want_metric:
return res
mf = {self._r: x.val[self._icov], self._icov: .5*x.val[self._icov]**(-2)}
metric = makeOp(MultiField.from_dict(mf))
......@@ -242,15 +238,9 @@ class PoissonianEnergy(EnergyOperator):
def apply(self, x):
self._check_input(x)
res = x.sum()
tmp = res.val.val if isinstance(res, Linearization) else res
# if we have no infinity here, we can continue with the calculation;
# otherwise we know that the result must also be infinity
if not np.isinf(tmp):
res = res - x.log().vdot(self._d)
if not isinstance(x, Linearization):
return Field.scalar(res)
if not x.want_metric:
fa = FieldAdapter(self._domain, 'foo')
res = (fa.sum() - fa.log().vdot(self._d))(fa.adjoint(x))
if not isinstance(x, Linearization) or not x.want_metric:
return res
metric = SandwichOperator.make(x.jac, makeOp(1./x.val))
return res.add_metric(metric)
......@@ -280,20 +270,20 @@ class InverseGammaLikelihood(EnergyOperator):
def __init__(self, beta, alpha=-0.5):
if not isinstance(beta, Field):
raise TypeError
self._domain = DomainTuple.make(beta.domain)
self._beta = beta
if np.isscalar(alpha):
alpha = Field(beta.domain, np.full(beta.shape, alpha))
elif not isinstance(alpha, Field):
raise TypeError
self._alphap1 = alpha+1
self._domain = DomainTuple.make(beta.domain)
def apply(self, x):
self._check_input(x)
res = x.log().vdot(self._alphap1) + (1./x).vdot(self._beta)
if not isinstance(x, Linearization):
return Field.scalar(res)
if not x.want_metric:
fa = FieldAdapter(self._domain, 'foo')
x = fa.adjoint(x)
res = (fa.log().vdot(self._alphap1) + fa.one_over().vdot(self._beta))(x)
if not isinstance(x, Linearization) or not x.want_metric:
return res
metric = SandwichOperator.make(x.jac, makeOp(self._alphap1/(x.val**2)))
return res.add_metric(metric)
......@@ -359,10 +349,11 @@ class BernoulliEnergy(EnergyOperator):
def apply(self, x):
self._check_input(x)
v = -(x.log().vdot(self._d) + (1. - x).log().vdot(1. - self._d))
if not isinstance(x, Linearization):
return Field.scalar(v)
if not x.want_metric:
iden = FieldAdapter(self._domain, 'foo')
from .adder import Adder
v = -iden.log().vdot(self._d) + (Adder(1, domain=self._domain) @ iden.scale(-1)).log().vdot(self._d-1.)
v = v(iden.adjoint(x))
if not isinstance(x, Linearization) or not x.want_metric:
return v
met = makeOp(1./(x.val*(1. - x.val)))
met = SandwichOperator.make(x.jac, met)
......@@ -413,13 +404,11 @@ class StandardHamiltonian(EnergyOperator):
def apply(self, x):
self._check_input(x)
if (self._ic_samp is None or not isinstance(x, Linearization)
or not x.want_metric):
return self._lh(x) + self._prior(x)
if (self._ic_samp is None or not isinstance(x, Linearization) or not x.want_metric):
return (self._lh + self._prior)(x)
else:
lhx, prx = self._lh(x), self._prior(x)
mtr = SamplingEnabler(lhx.metric, prx.metric,
self._ic_samp)
mtr = SamplingEnabler(lhx.metric, prx.metric, self._ic_samp)
return (lhx + prx).add_metric(mtr)
def __repr__(self):
......@@ -461,5 +450,11 @@ class AveragedEnergy(EnergyOperator):
def apply(self, x):
self._check_input(x)
mymap = map(lambda v: self._h(x + v), self._res_samples)
return utilities.my_sum(mymap)*(1./len(self._res_samples))
if isinstance(self._domain, MultiDomain):
iden = ScalingOperator(self._domain, 1.)
else:
iden = FieldAdapter(self._domain, 'foo')
x = iden.adjoint(x)
from .adder import Adder
mymap = map(lambda v: self._h(Adder(v) @ iden), self._res_samples)
return utilities.my_sum(mymap).scale(1./len(self._res_samples))(x)
......@@ -71,7 +71,16 @@ class Operator(metaclass=NiftyMeta):
return ContractionOperator(self.target, spaces)(self)
def vdot(self, other):
return (self.conjugate()*other).sum()
from ..field import Field
from ..multi_field import MultiField
from ..sugar import makeOp
if isinstance(other, Operator):
res = self.conjugate()*other
elif isinstance(other, (Field, MultiField)):
res = makeOp(other) @ self.conjugate()
else:
raise TypeError
return res.sum()
@property
def real(self):
......
......@@ -25,11 +25,11 @@ from itertools import product
# hopefully be fixed in the future.
# https://docs.pytest.org/en/latest/proposals/parametrize_with_fixtures.html
SPACES = [
ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)
]
SPACES = [ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)]
for sp in SPACES[:3]:
SPACES.append(ift.MultiDomain.make({'asdf': sp}))
SEEDS = [4, 78, 23]
PARAMS = product(SEEDS, SPACES)
pmp = pytest.mark.parametrize
......@@ -39,16 +39,7 @@ pmp = pytest.mark.parametrize
def field(request):
np.random.seed(request.param[0])
S = ift.ScalingOperator(request.param[1], 1.)
s = S.draw_sample()
return ift.MultiField.from_dict({'s1': s})['s1']
def test_variablecovariancegaussian(field):
dc = {'a': field, 'b': field.exp()}
mf = ift.MultiField.from_dict(dc)
energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b')
ift.extra.check_jacobian_consistency(energy, mf, tol=1e-6)
energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample()
return S.draw_sample()
def test_gaussian(field):
......@@ -56,26 +47,52 @@ def test_gaussian(field):
ift.extra.check_jacobian_consistency(energy, field)
@pmp('icov', [lambda dom:ift.ScalingOperator(dom, 1.),
lambda dom:ift.SandwichOperator.make(ift.GeometryRemover(dom))])
def test_ScaledEnergy(field, icov):
icov = icov(field.domain)
def test_ScaledEnergy(field):
icov = ift.ScalingOperator(field.domain, 1.2)
energy = ift.GaussianEnergy(inverse_covariance=icov)
ift.extra.check_jacobian_consistency(energy.scale(0.3), field)
lin = ift.Linearization.make_var(field, want_metric=True)
met1 = energy(lin).metric
met2 = energy.scale(0.3)(lin).metric
np.testing.assert_allclose(met1(field).val, met2(field).val/0.3, rtol=1e-12)
res1 = met1(field)
res2 = met2(field)/0.3
ift.extra.assert_allclose(res1, res2, 0, 1e-12)
met2.draw_sample()
def test_studentt(field):
if isinstance(field.domain, ift.MultiDomain):
return
energy = ift.StudentTEnergy(domain=field.domain, theta=.5)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-6)
def test_hamiltonian_and_KL(field):
field = field.exp()
space = field.domain
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.StandardHamiltonian(lh)
ift.extra.check_jacobian_consistency(hamiltonian, field)
S = ift.ScalingOperator(space, 1.)
samps = [S.draw_sample() for i in range(3)]
kl = ift.AveragedEnergy(hamiltonian, samps)
ift.extra.check_jacobian_consistency(kl, field)
def test_variablecovariancegaussian(field):
if isinstance(field.domain, ift.MultiDomain):
return
dc = {'a': field, 'b': field.exp()}
mf = ift.MultiField.from_dict(dc)
energy = ift.VariableCovarianceGaussianEnergy(field.domain, 'a', 'b')
ift.extra.check_jacobian_consistency(energy, mf, tol=1e-6)
energy(ift.Linearization.make_var(mf, want_metric=True)).metric.draw_sample()
def test_inverse_gamma(field):
if isinstance(field.domain, ift.MultiDomain):
return
field = field.exp()
space = field.domain
d = np.random.normal(10, size=space.shape)**2
......@@ -85,6 +102,8 @@ def test_inverse_gamma(field):
def testPoissonian(field):
if isinstance(field.domain, ift.MultiDomain):
return
field = field.exp()
space = field.domain
d = np.random.poisson(120, size=space.shape)
......@@ -93,19 +112,9 @@ def testPoissonian(field):
ift.extra.check_jacobian_consistency(energy, field, tol=1e-7)
def test_hamiltonian_and_KL(field):
field = field.exp()
space = field.domain
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.StandardHamiltonian(lh)
ift.extra.check_jacobian_consistency(hamiltonian, field)
S = ift.ScalingOperator(space, 1.)
samps = [S.draw_sample() for i in range(3)]
kl = ift.AveragedEnergy(hamiltonian, samps)
ift.extra.check_jacobian_consistency(kl, field)
def test_bernoulli(field):
if isinstance(field.domain, ift.MultiDomain):
return
field = field.sigmoid()
space = field.domain
d = np.random.binomial(1, 0.1, size=space.shape)
......
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