Merge branch 'NIFTy_5' into domain_tuple_ops

demos/bernoulli_demo.py

@@ -73,14 +73,15 @@ if __name__ == '__main__':
 
     # Minimize the Hamiltonian
     H = ift.Hamiltonian(likelihood, ic_sampling)
-    H = ift.EnergyAdapter(position, H)
+    H = ift.EnergyAdapter(position, H, want_metric=True)
     # minimizer = ift.L_BFGS(ic_newton)
     H, convergence = minimizer(H)
 
     reconstruction = sky(H.position)
 
-    ift.plot(reconstruction, title='reconstruction')
-    ift.plot(GR.adjoint_times(data), title='data')
-    ift.plot(sky(mock_position), title='truth')
-    ift.plot_finish(nx=3, xsize=16, ysize=5, title="results",
-                    name="bernoulli.png")
+    plot = ift.Plot()
+    plot.add(reconstruction, title='reconstruction')
+    plot.add(GR.adjoint_times(data), title='data')
+    plot.add(sky(mock_position), title='truth')
+    plot.output(nx=3, xsize=16, ysize=5, title="results",
+                name="bernoulli.png")

demos/getting_started_1.py

@@ -103,18 +103,17 @@ if __name__ == '__main__':
 
     # PLOTTING
     rg = isinstance(position_space, ift.RGSpace)
+    plot = ift.Plot()
     if rg and len(position_space.shape) == 1:
-        ift.plot([HT(MOCK_SIGNAL), GR.adjoint(data), HT(m)],
+        plot.add([HT(MOCK_SIGNAL), GR.adjoint(data), HT(m)],
                  label=['Mock signal', 'Data', 'Reconstruction'],
                  alpha=[1, .3, 1])
-        ift.plot(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals')
-        ift.plot_finish(nx=2, ny=1, xsize=10, ysize=4,
-                        title="getting_started_1")
+        plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals')
+        plot.output(nx=2, ny=1, xsize=10, ysize=4, title="getting_started_1")
     else:
-        ift.plot(HT(MOCK_SIGNAL), title='Mock Signal')
-        ift.plot(mask_to_nan(mask, (GR(Mask)).adjoint(data)),
+        plot.add(HT(MOCK_SIGNAL), title='Mock Signal')
+        plot.add(mask_to_nan(mask, (GR(Mask)).adjoint(data)),
                  title='Data')
-        ift.plot(HT(m), title='Reconstruction')
-        ift.plot(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals')
-        ift.plot_finish(nx=2, ny=2, xsize=10, ysize=10,
-                        title="getting_started_1")
+        plot.add(HT(m), title='Reconstruction')
+        plot.add(mask_to_nan(mask, HT(m-MOCK_SIGNAL)), title='Residuals')
+        plot.output(nx=2, ny=2, xsize=10, ysize=10, title="getting_started_1")

demos/getting_started_2.py

@@ -93,14 +93,15 @@ if __name__ == '__main__':
 
     # Minimize the Hamiltonian
     H = ift.Hamiltonian(likelihood)
-    H = ift.EnergyAdapter(position, H)
+    H = ift.EnergyAdapter(position, H, want_metric=True)
     H, convergence = minimizer(H)
 
     # Plot results
     signal = sky(mock_position)
     reconst = sky(H.position)
-    ift.plot(signal, title='Signal')
-    ift.plot(GR.adjoint(data), title='Data')
-    ift.plot(reconst, title='Reconstruction')
-    ift.plot(reconst - signal, title='Residuals')
-    ift.plot_finish(name='getting_started_2.png', xsize=16, ysize=16)
+    plot = ift.Plot()
+    plot.add(signal, title='Signal')
+    plot.add(GR.adjoint(data), title='Data')
+    plot.add(reconst, title='Reconstruction')
+    plot.add(reconst - signal, title='Residuals')
+    plot.output(name='getting_started_2.png', xsize=16, ysize=16)

demos/getting_started_3.py

@@ -41,11 +41,12 @@ if __name__ == '__main__':
     power_space = A.target[0]
     power_distributor = ift.PowerDistributor(harmonic_space, power_space)
     dummy = ift.Field.from_random('normal', harmonic_space)
-    domain = ift.MultiDomain.union(
-        (A.domain, ift.MultiDomain.make({'xi': harmonic_space})))
+    domain = ift.MultiDomain.union((A.domain,
+                                    ift.MultiDomain.make({
+                                        'xi': harmonic_space
+                                    })))
 
-    correlated_field = ht(
-        power_distributor(A)*ift.FieldAdapter(domain, "xi"))
+    correlated_field = ht(power_distributor(A)*ift.FieldAdapter(domain, "xi"))
     # alternatively to the block above one can do:
     # correlated_field = ift.CorrelatedField(position_space, A)
 
@@ -54,8 +55,7 @@ if __name__ == '__main__':
 
     # Building the Line of Sight response
     LOS_starts, LOS_ends = get_random_LOS(100)
-    R = ift.LOSResponse(position_space, starts=LOS_starts,
-                        ends=LOS_ends)
+    R = ift.LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
     # build signal response model and model likelihood
     signal_response = R(signal)
     # specify noise
@@ -68,8 +68,7 @@ if __name__ == '__main__':
     data = signal_response(MOCK_POSITION) + N.draw_sample()
 
     # set up model likelihood
-    likelihood = ift.GaussianEnergy(
-        mean=data, covariance=N)(signal_response)
+    likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
 
     # set up minimization and inversion schemes
     ic_sampling = ift.GradientNormController(iteration_limit=100)
@@ -83,34 +82,33 @@ if __name__ == '__main__':
     INITIAL_POSITION = ift.from_random('normal', domain)
     position = INITIAL_POSITION
 
-    ift.plot(signal(MOCK_POSITION), title='ground truth')
-    ift.plot(R.adjoint_times(data), title='data')
-    ift.plot([A(MOCK_POSITION)], title='power')
-    ift.plot_finish(nx=3, xsize=16, ysize=5, title="setup", name="setup.png")
+    plot = ift.Plot()
+    plot.add(signal(MOCK_POSITION), title='Ground Truth')
+    plot.add(R.adjoint_times(data), title='Data')
+    plot.add([A(MOCK_POSITION)], title='Power Spectrum')
+    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
 
     # number of samples used to estimate the KL
     N_samples = 20
     for i in range(2):
-        metric = H(ift.Linearization.make_var(position)).metric
-        samples = [metric.draw_sample(from_inverse=True)
-                   for _ in range(N_samples)]
-
-        KL = ift.SampledKullbachLeiblerDivergence(H, samples)
-        KL = ift.EnergyAdapter(position, KL)
+        KL = ift.KL_Energy(position, H, N_samples, want_metric=True)
         KL, convergence = minimizer(KL)
         position = KL.position
 
-        ift.plot(signal(position), title="reconstruction")
-        ift.plot([A(position), A(MOCK_POSITION)], title="power")
-        ift.plot_finish(nx=2, xsize=12, ysize=6, title="loop", name="loop.png")
+        plot = ift.Plot()
+        plot.add(signal(KL.position), title="reconstruction")
+        plot.add([A(KL.position), A(MOCK_POSITION)], title="power")
+        plot.output(ny=1, ysize=6, xsize=16, name="loop.png")
 
+    plot = ift.Plot()
     sc = ift.StatCalculator()
-    for sample in samples:
-        sc.add(signal(sample+position))
-    ift.plot(sc.mean, title="mean")
-    ift.plot(ift.sqrt(sc.var), title="std deviation")
-
-    powers = [A(s+position) for s in samples]
-    ift.plot([A(position), A(MOCK_POSITION)]+powers, title="power")
-    ift.plot_finish(nx=3, xsize=16, ysize=5, title="results",
-                    name="results.png")
+    for sample in KL.samples:
+        sc.add(signal(sample+KL.position))
+    plot.add(sc.mean, title="Posterior Mean")
+    plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
+
+    powers = [A(s+KL.position) for s in KL.samples]
+    plot.add(
+        [A(KL.position), A(MOCK_POSITION)]+powers,
+        title="Sampled Posterior Power Spectrum")
+    plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")

demos/plot_test.py

@@ -20,21 +20,24 @@ def plot_test():
 
     # Start various plotting tests
 
-    ift.plot(field_rg1_1, title='Single plot')
-    ift.plot_finish()
-
-    ift.plot(field_rg2, title='2d rg')
-    ift.plot([field_rg1_1, field_rg1_2], title='list 1d rg', label=['1', '2'])
-    ift.plot(field_rg1_2, title='1d rg, xmin, ymin', xmin=0.5, ymin=0.,
+    plot = ift.Plot()
+    plot.add(field_rg1_1, title='Single plot')
+    plot.output()
+
+    plot = ift.Plot()
+    plot.add(field_rg2, title='2d rg')
+    plot.add([field_rg1_1, field_rg1_2], title='list 1d rg', label=['1', '2'])
+    plot.add(field_rg1_2, title='1d rg, xmin, ymin', xmin=0.5, ymin=0.,
              xlabel='xmin=0.5', ylabel='ymin=0')
-    ift.plot_finish(title='Three plots')
-
-    ift.plot(field_hp, title='HP planck-color', colormap='Planck-like')
-    ift.plot(field_rg1_2, title='1d rg')
-    ift.plot(field_ps)
-    ift.plot(field_gl, title='GL')
-    ift.plot(field_rg2, title='2d rg')
-    ift.plot_finish(nx=2, ny=3, title='Five plots')
+    plot.output(title='Three plots')
+
+    plot = ift.Plot()
+    plot.add(field_hp, title='HP planck-color', colormap='Planck-like')
+    plot.add(field_rg1_2, title='1d rg')
+    plot.add(field_ps)
+    plot.add(field_gl, title='GL')
+    plot.add(field_rg2, title='2d rg')
+    plot.output(nx=2, ny=3, title='Five plots')
 
 
 if __name__ == '__main__':

demos/polynomial_fit.py

@@ -86,15 +86,16 @@ N = ift.DiagonalOperator(ift.from_global_data(d_space, var))
 
 IC = ift.GradientNormController(tol_abs_gradnorm=1e-8)
 likelihood = ift.GaussianEnergy(d, N)(R)
-H = ift.Hamiltonian(likelihood, IC)
-H = ift.EnergyAdapter(params, H, IC)
+Ham = ift.Hamiltonian(likelihood, IC)
+H = ift.EnergyAdapter(params, Ham, want_metric=True)
 
 # Minimize
 minimizer = ift.NewtonCG(IC)
 H, _ = minimizer(H)
 
 # Draw posterior samples
-samples = [H.metric.draw_sample(from_inverse=True) + H.position
+metric = Ham(ift.Linearization.make_var(H.position, want_metric=True)).metric
+samples = [metric.draw_sample(from_inverse=True) + H.position
            for _ in range(N_samples)]
 
 # Plotting

nifty5/__init__.py

@@ -22,7 +22,8 @@ from .operators.operator import Operator
 from .operators.central_zero_padder import CentralZeroPadder
 from .operators.diagonal_operator import DiagonalOperator
 from .operators.distributors import DOFDistributor, PowerDistributor
-from .operators.domain_tuple_operators import DomainDistributor, DomainTupleFieldInserter
+from .operators.domain_tuple_operators import DomainTupleFieldInserter
+from .operators.contraction_operator import ContractionOperator
 from .operators.endomorphic_operator import EndomorphicOperator
 from .operators.exp_transform import ExpTransform
 from .operators.harmonic_operators import (
@@ -67,9 +68,10 @@ from .minimization.energy import Energy
 from .minimization.quadratic_energy import QuadraticEnergy
 from .minimization.line_energy import LineEnergy
 from .minimization.energy_adapter import EnergyAdapter
+from .minimization.kl_energy import KL_Energy
 
 from .sugar import *
-from .plotting.plot import plot, plot_finish
+from .plotting.plot import Plot
 
 from .library.amplitude_model import AmplitudeModel
 from .library.inverse_gamma_model import InverseGammaModel

nifty5/extra/energy_and_model_tests.py

@@ -41,7 +41,7 @@ def _get_acceptable_location(op, loc, lin):
     for i in range(50):
         try:
             loc2 = loc+dir
-            lin2 = op(Linearization.make_var(loc2))
+            lin2 = op(Linearization.make_var(loc2, lin.want_metric))
             if np.isfinite(lin2.val.sum()) and abs(lin2.val.sum()) < 1e20:
                 break
         except FloatingPointError:
@@ -54,14 +54,14 @@ def _get_acceptable_location(op, loc, lin):
 
 def _check_consistency(op, loc, tol, ntries, do_metric):
     for _ in range(ntries):
-        lin = op(Linearization.make_var(loc))
+        lin = op(Linearization.make_var(loc, do_metric))
         loc2, lin2 = _get_acceptable_location(op, loc, lin)
         dir = loc2-loc
         locnext = loc2
         dirnorm = dir.norm()
         for i in range(50):
             locmid = loc + 0.5*dir
-            linmid = op(Linearization.make_var(locmid))
+            linmid = op(Linearization.make_var(locmid, do_metric))
             dirder = linmid.jac(dir)
             numgrad = (lin2.val-lin.val)
             xtol = tol * dirder.norm() / np.sqrt(dirder.size)

nifty5/library/correlated_fields.py

@@ -21,8 +21,8 @@ from __future__ import absolute_import, division, print_function
 from ..compat import *
 from ..domain_tuple import DomainTuple
 from ..multi_domain import MultiDomain
+from ..operators.contraction_operator import ContractionOperator
 from ..operators.distributors import PowerDistributor
-from ..operators.domain_tuple_operators import DomainDistributor
 from ..operators.harmonic_operators import HarmonicTransformOperator
 from ..operators.simple_linear_operators import FieldAdapter
 from ..sugar import exp
@@ -65,8 +65,8 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
     pd_energy = PowerDistributor(pd_spatial.domain, p_space_energy, 1)
     pd = pd_spatial(pd_energy)
 
-    dom_distr_spatial = DomainDistributor(pd.domain, 0)
-    dom_distr_energy = DomainDistributor(pd.domain, 1)
+    dom_distr_spatial = ContractionOperator(pd.domain, 0).adjoint
+    dom_distr_energy = ContractionOperator(pd.domain, 1).adjoint
 
     a_spatial = dom_distr_spatial(amplitude_model_spatial)
     a_energy = dom_distr_energy(amplitude_model_energy)

nifty5/library/inverse_gamma_model.py

@@ -53,7 +53,7 @@ class InverseGammaModel(Operator):
         outer = 1/outer_inv
         jac = makeOp(Field.from_local_data(self._domain, inner*outer))
         jac = jac(x.jac)
-        return Linearization(points, jac)
+        return x.new(points, jac)
 
     @staticmethod
     def IG(field, alpha, q):

nifty5/linearization.py

@@ -9,13 +9,17 @@ from .sugar import makeOp
 
 
 class Linearization(object):
-    def __init__(self, val, jac, metric=None):
+    def __init__(self, val, jac, metric=None, want_metric=False):
         self._val = val
         self._jac = jac
         if self._val.domain != self._jac.target:
             raise ValueError("domain mismatch")
+        self._want_metric = want_metric
         self._metric = metric
 
+    def new(self, val, jac, metric=None):
+        return Linearization(val, jac, metric, self._want_metric)
+
     @property
     def domain(self):
         return self._jac.domain
@@ -37,6 +41,10 @@ class Linearization(object):
         """Only available if target is a scalar"""
         return self._jac.adjoint_times(Field.scalar(1.))
 
+    @property
+    def want_metric(self):
+        return self._want_metric
+
     @property
     def metric(self):
         """Only available if target is a scalar"""
@@ -44,35 +52,34 @@ class Linearization(object):
 
     def __getitem__(self, name):
         from .operators.simple_linear_operators import FieldAdapter
-        return Linearization(self._val[name], FieldAdapter(self.domain, name))
+        return self.new(self._val[name], FieldAdapter(self.domain, name))
 
     def __neg__(self):
-        return Linearization(
-            -self._val, -self._jac,
-            None if self._metric is None else -self._metric)
+        return self.new(-self._val, -self._jac,
+                        None if self._metric is None else -self._metric)
 
     def conjugate(self):
-        return Linearization(
+        return self.new(
             self._val.conjugate(), self._jac.conjugate(),
             None if self._metric is None else self._metric.conjugate())
 
     @property
     def real(self):
-        return Linearization(self._val.real, self._jac.real)
+        return self.new(self._val.real, self._jac.real)
 
     def _myadd(self, other, neg):
         if isinstance(other, Linearization):
             met = None
             if self._metric is not None and other._metric is not None:
                 met = self._metric._myadd(other._metric, neg)
-            return Linearization(
+            return self.new(
                 self._val.flexible_addsub(other._val, neg),
                 self._jac._myadd(other._jac, neg), met)
         if isinstance(other, (int, float, complex, Field, MultiField)):
             if neg:
-                return Linearization(self._val-other, self._jac, self._metric)
+                return self.new(self._val-other, self._jac, self._metric)
             else:
-                return Linearization(self._val+other, self._jac, self._metric)
+                return self.new(self._val+other, self._jac, self._metric)
 
     def __add__(self, other):
         return self._myadd(other, False)
@@ -91,7 +98,7 @@ class Linearization(object):
         if isinstance(other, Linearization):
             if self.target != other.target:
                 raise ValueError("domain mismatch")
-            return Linearization(
+            return self.new(
                 self._val*other._val,
                 (makeOp(other._val)(self._jac))._myadd(
                  makeOp(self._val)(other._jac), False))
@@ -99,11 +106,11 @@ class Linearization(object):
             if other == 1:
                 return self
             met = None if self._metric is None else self._metric.scale(other)
-            return Linearization(self._val*other, self._jac.scale(other), met)
+            return self.new(self._val*other, self._jac.scale(other), met)
         if isinstance(other, (Field, MultiField)):
             if self.target != other.domain:
                 raise ValueError("domain mismatch")
-            return Linearization(self._val*other, makeOp(other)(self._jac))
+            return self.new(self._val*other, makeOp(other)(self._jac))
 
     def __rmul__(self, other):
         return self.__mul__(other)
@@ -111,46 +118,48 @@ class Linearization(object):
     def vdot(self, other):
         from .operators.simple_linear_operators import VdotOperator
         if isinstance(other, (Field, MultiField)):
-            return Linearization(
+            return self.new(
                 Field.scalar(self._val.vdot(other)),
                 VdotOperator(other)(self._jac))
-        return Linearization(
+        return self.new(
             Field.scalar(self._val.vdot(other._val)),
             VdotOperator(self._val)(other._jac) +
             VdotOperator(other._val)(self._jac))
 
     def sum(self):
         from .operators.simple_linear_operators import SumReductionOperator
-        return Linearization(
+        return self.new(
             Field.scalar(self._val.sum()),
             SumReductionOperator(self._jac.target)(self._jac))
 
     def exp(self):
         tmp = self._val.exp()
-        return Linearization(tmp, makeOp(tmp)(self._jac))
+        return self.new(tmp, makeOp(tmp)(self._jac))
 
     def log(self):
         tmp = self._val.log()
-        return Linearization(tmp, makeOp(1./self._val)(self._jac))
+        return self.new(tmp, makeOp(1./self._val)(self._jac))
 
     def tanh(self):
         tmp = self._val.tanh()
-        return Linearization(tmp, makeOp(1.-tmp**2)(self._jac))
+        return self.new(tmp, makeOp(1.-tmp**2)(self._jac))
 
     def positive_tanh(self):
         tmp = self._val.tanh()
         tmp2 = 0.5*(1.+tmp)
-        return Linearization(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
+        return self.new(tmp2, makeOp(0.5*(1.-tmp**2))(self._jac))
 
     def add_metric(self, metric):
-        return Linearization(self._val, self._jac, metric)
+        return self.new(self._val, self._jac, metric)
 
     @staticmethod
-    def make_var(field):
+    def make_var(field, want_metric=False):
         from .operators.scaling_operator import ScalingOperator
-        return Linearization(field, ScalingOperator(1., field.domain))
+        return Linearization(field, ScalingOperator(1., field.domain),
+                             want_metric=want_metric)
 
     @staticmethod
-    def make_const(field):
+    def make_const(field, want_metric=False):
         from .operators.simple_linear_operators import NullOperator
-        return Linearization(field, NullOperator(field.domain, field.domain))
+        return Linearization(field, NullOperator(field.domain, field.domain),
+                             want_metric=want_metric)

nifty5/minimization/conjugate_gradient.py

@@ -75,7 +75,7 @@ class ConjugateGradient(Minimizer):
             return energy, controller.CONVERGED
 
         while True:
-            q = energy.metric(d)
+            q = energy.apply_metric(d)
             ddotq = d.vdot(q).real
             if ddotq == 0.:
                 logger.error("Error: ConjugateGradient: ddotq==0.")

nifty5/minimization/descent_minimizers.py

@@ -180,7 +180,7 @@ class NewtonCG(DescentMinimizer):
         while True:
             if abs(ri).sum() <= termcond:
                 return xsupi
-            Ap = energy.metric(psupi)
+            Ap = energy.apply_metric(psupi)
             # check curvature
             curv = psupi.vdot(Ap)
             if 0 <= curv <= 3*float64eps:

nifty5/minimization/energy.py

@@ -109,6 +109,20 @@ class Energy(NiftyMetaBase()):
         """
         raise NotImplementedError
 
+    def apply_metric(self, x):
+        """
+        Parameters
+        ----------
+        x: Field/MultiField
+            Argument for the metric operator
+
+        Returns
+        -------
+        Field/MultiField:
+            Output of the metric operator
+        """
+        raise NotImplementedError
+
     def longest_step(self, dir):
         """Returns the longest allowed step size along `dir`
 

nifty5/minimization/energy_adapter.py

@@ -8,58 +8,38 @@ from ..operators.scaling_operator import ScalingOperator
 
 
 class EnergyAdapter(Energy):
-    def __init__(self, position, op, controller=None, preconditioner=None,
-                 constants=[]):
+    def __init__(self, position, op, constants=[], want_metric=False):
         super(EnergyAdapter, self).__init__(position)
         self._op = op
-        self._val = self._grad = self._metric = None
-        self._controller = controller
-        self._preconditioner = preconditioner
         self._constants = constants
-
-    def at(self, position):
-        return EnergyAdapter(position, self._op, self._controller,
-                             self._preconditioner, self._constants)
-
-    def _fill_all(self):
+        self._want_metric = want_metric
         if len(self._constants) == 0:
-            tmp = self._op(Linearization.make_var(self._position))
+            tmp = self._op(Linearization.make_var(self._position, want_metric))
         else:
             ops = [ScalingOperator(0. if key in self._constants else 1., dom)
                    for key, dom in self._position.domain.items()]
             bdop = BlockDiagonalOperator(self._position.domain, tuple(ops))
-            tmp = self._op(Linearization(self._position, bdop))
+            tmp = self._op(Linearization(self._position, bdop,
+                                         want_metric=want_metric))
         self._val = tmp.val.local_data[()]
         self._grad = tmp.gradient
-        if self._controller is not None:
-            from ..operators.linear_operator import LinearOperator
-            from ..operators.inversion_enabler import InversionEnabler
+        self._metric = tmp._metric
 
-            if self._preconditioner is None:
-                precond = None
-            elif isinstance(self._preconditioner, LinearOperator):
-                precond = self._preconditioner
-            elif isinstance(self._preconditioner, Energy):
-                precond = self._preconditioner.at(self._position).metric
-            self._metric = InversionEnabler(tmp._metric, self._controller,
-                                            precond)
-        else:
-            self._metric = tmp._metric
+    def at(self, position):
+        return EnergyAdapter(position, self._op, self._constants,
+                             self._want_metric)
 
     @property
     def value(self):
-        if self._val is None:
-            self._val = self._op(self._position).local_data[()]
         return self._val
 
     @property
     def gradient(self):
-        if self._grad is None:
-            self._fill_all()
         return self._grad
 
     @property
     def metric(self):
-        if self._metric is None:
-            self._fill_all()
         return self._metric
+
+    def apply_metric(self, x):
+        return self._metric(x)

nifty5/minimization/kl_energy.py

+from __future__ import absolute_import, division, print_function
+
+from ..compat import *
+from .energy import Energy
+from ..linearization import Linearization
+from ..operators.scaling_operator import ScalingOperator
+from ..operators.block_diagonal_operator import BlockDiagonalOperator
+from .. import utilities