Merge branch 'NIFTy_5' into energy_cleanup

.gitlab-ci.yml

@@ -15,6 +15,8 @@ build_docker_from_scratch:
     - schedules
   image: docker:stable
   stage: build_docker
+  before_script:
+    - ls
   script:
     - docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
     - docker build -t $CONTAINER_TEST_IMAGE --no-cache .
@@ -25,6 +27,8 @@ build_docker_from_cache:
     - schedules
   image: docker:stable
   stage: build_docker
+  before_script:
+    - ls
   script:
     - docker login -u gitlab-ci-token -p $CI_BUILD_TOKEN gitlab-registry.mpcdf.mpg.de
     - docker build -t $CONTAINER_TEST_IMAGE .
@@ -33,7 +37,6 @@ build_docker_from_cache:
 test_python2_with_coverage:
   stage: test
   script:
-    - python setup.py install --user -f
     - mpiexec -n 2 --bind-to none nosetests -q 2> /dev/null
     - nosetests -q --with-coverage --cover-package=nifty5 --cover-erase
     - >
@@ -44,12 +47,13 @@ test_python2_with_coverage:
 test_python3:
   stage: test
   script:
-    - python3 setup.py install --user -f
     - nosetests3 -q
     - mpiexec -n 2 --bind-to none nosetests3 -q 2> /dev/null
 
 pages:
   stage: release
+  before_script:
+    - ls
   script:
     - python setup.py install --user -f
     - sh docs/generate.sh
@@ -62,13 +66,31 @@ pages:
 
 before_script:
   - export MPLBACKEND="agg"
+  - python setup.py install --user -f
+  - python3 setup.py install --user -f
 
 run_ipynb:
   stage: demo_runs
   script:
-    - python setup.py install --user -f
-    - python3 setup.py install --user -f
     - jupyter nbconvert --execute --ExecutePreprocessor.timeout=None demos/Wiener_Filter.ipynb
   artifacts:
     paths:
       - '*.png'
+
+run_getting_started_1:
+  stage: demo_runs
+  script:
+    - python demos/getting_started_1.py
+    - python3 demos/getting_started_1.py
+
+run_getting_started_2:
+  stage: demo_runs
+  script:
+    - python demos/getting_started_2.py
+    - python3 demos/getting_started_2.py
+
+run_getting_started_3:
+  stage: demo_runs
+  script:
+    - python demos/getting_started_3.py
+    - python3 demos/getting_started_3.py

demos/getting_started_1.py

 import nifty5 as ift
 import numpy as np
-from global_newton.models_other.apply_data import ApplyData
-from global_newton.models_energy.hamiltonian import Hamiltonian
-from nifty5 import GaussianEnergy
 
 
-if __name__ == '__main__':
-    # s_space = ift.RGSpace([1024])
-    s_space = ift.RGSpace([128,128])
-    # s_space = ift.HPSpace(64)
+def make_chess_mask():
+    mask = np.ones(position_space.shape)
+    for i in range(4):
+        for j in range(4):
+            if (i+j) % 2 == 0:
+                mask[i*128//4:(i+1)*128//4, j*128//4:(j+1)*128//4] = 0
+    return mask
 
-    h_space = s_space.get_default_codomain()
-    total_domain = ift.MultiDomain.make({'xi': h_space})
-    HT = ift.HarmonicTransformOperator(h_space, s_space)
 
-    def sqrtpspec(k):
-        return 16. / (20.+k**2)
+def make_random_mask():
+    mask = ift.from_random('pm1', position_space)
+    mask = (mask+1)/2
+    return mask.to_global_data()
 
-    GR = ift.GeometryRemover(s_space)
 
-    d_space = GR.target
-    B = ift.FFTSmoothingOperator(s_space,0.1)
-    mask = np.ones(s_space.shape)
-    mask[64:89,76:100] = 0.
-    mask = ift.Field(s_space,val=mask)
-    Mask = ift.DiagonalOperator(mask)
-    R = GR * Mask * B
-    noise = 1.
-    N = ift.ScalingOperator(noise, d_space)
+if __name__ == '__main__':
+    # # description of the tutorial ###
 
-    p_space = ift.PowerSpace(h_space)
-    pd = ift.PowerDistributor(h_space, p_space)
-    position = ift.from_random('normal', total_domain)
-    xi = ift.Variable(position)['xi']
-    a = ift.Constant(position, ift.PS_field(p_space, sqrtpspec))
-    A = pd(a)
-    s_h = A * xi
-    s = HT(s_h)
-    Rs = R(s)
+    # Choose problem geometry and masking
 
+    # One dimensional regular grid
+    position_space = ift.RGSpace([1024])
+    mask = np.ones(position_space.shape)
 
+    # # Two dimensional regular grid with chess mask
+    # position_space = ift.RGSpace([128,128])
+    # mask = make_chess_mask()
 
-    MOCK_POSITION = ift.from_random('normal',total_domain)
-    data = Rs.at(MOCK_POSITION).value + N.draw_sample()
+    # # Sphere with half of its locations randomly masked
+    # position_space = ift.HPSpace(128)
+    # mask = make_random_mask()
 
-    NWR = ApplyData(data, ift.Field(d_space,val=noise), Rs)
+    harmonic_space = position_space.get_default_codomain()
+    HT = ift.HarmonicTransformOperator(harmonic_space, target=position_space)
 
-    INITIAL_POSITION = ift.from_random('normal',total_domain)
-    likelihood = GaussianEnergy(INITIAL_POSITION, NWR)
+    # set correlation structure with a power spectrum and build
+    # prior correlation covariance
+    def power_spectrum(k):
+        return 100. / (20.+k**3)
+    power_space = ift.PowerSpace(harmonic_space)
+    PD = ift.PowerDistributor(harmonic_space, power_space)
+    prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
 
-    IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
-    inverter = ift.ConjugateGradient(controller=IC)
-    IC2 = ift.GradientNormController(name='Newton', iteration_limit=15)
-    minimizer = ift.RelaxedNewton(IC2)
+    S = ift.DiagonalOperator(prior_correlation_structure)
 
+    # build instrument response consisting of a discretization, mask
+    # and harmonic transformaion
+    GR = ift.GeometryRemover(position_space)
+    mask = ift.Field.from_global_data(position_space, mask)
+    Mask = ift.DiagonalOperator(mask)
+    R = GR * Mask * HT
 
-    H = Hamiltonian(likelihood, inverter)
-    H, convergence = minimizer(H)
-    result = s.at(H.position).value
+    data_space = GR.target
+
+    # setting the noise covariance
+    noise = 5.
+    N = ift.ScalingOperator(noise, data_space)
+
+    # creating mock data
+    MOCK_SIGNAL = S.draw_sample()
+    MOCK_NOISE = N.draw_sample()
+    data = R(MOCK_SIGNAL) + MOCK_NOISE
+
+    # building propagator D and information source j
+    j = R.adjoint_times(N.inverse_times(data))
+    D_inv = R.adjoint * N.inverse * R + S.inverse
+    # make it invertible
+    IC = ift.GradientNormController(iteration_limit=500, tol_abs_gradnorm=1e-3)
+    D = ift.InversionEnabler(D_inv, IC, approximation=S.inverse).inverse
 
+    # WIENER FILTER
+    m = D(j)
 
+    # PLOTTING
+    # Truth, data, reconstruction, residuals

demos/getting_started_2.py

+import nifty5 as ift
+import numpy as np
+
+
+def get_2D_exposure():
+    x_shape, y_shape = position_space.shape
+
+    exposure = np.ones(position_space.shape)
+    exposure[x_shape//3:x_shape//2, :] *= 2.
+    exposure[x_shape*4//5:x_shape, :] *= .1
+    exposure[x_shape//2:x_shape*3//2, :] *= 3.
+    exposure[:, x_shape//3:x_shape//2] *= 2.
+    exposure[:, x_shape*4//5:x_shape] *= .1
+    exposure[:, x_shape//2:x_shape*3//2] *= 3.
+
+    exposure = ift.Field.from_global_data(position_space, exposure)
+    return exposure
+
+
+if __name__ == '__main__':
+    # ABOUT THIS CODE
+    np.random.seed(41)
+
+    # Set up the position space of the signal
+    #
+    # # One dimensional regular grid with uniform exposure
+    # position_space = ift.RGSpace(1024)
+    # exposure = np.ones(position_space.shape)
+
+    # Two-dimensional regular grid with inhomogeneous exposure
+    position_space = ift.RGSpace([512, 512])
+    exposure = get_2D_exposure()
+
+    # # Sphere with with uniform exposure
+    # position_space = ift.HPSpace(128)
+    # exposure = ift.Field.full(position_space, 1.)
+
+    # Defining harmonic space and transform
+    harmonic_space = position_space.get_default_codomain()
+    HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
+
+    domain = ift.MultiDomain.make({'xi': harmonic_space})
+    position = ift.from_random('normal', domain)
+
+    # Define power spectrum and amplitudes
+    def sqrtpspec(k):
+        return 1. / (20. + k**2)
+
+    p_space = ift.PowerSpace(harmonic_space)
+    pd = ift.PowerDistributor(harmonic_space, p_space)
+    a = ift.PS_field(p_space, sqrtpspec)
+    A = pd(a)
+
+    # Set up a sky model
+    xi = ift.Variable(position)['xi']
+    logsky_h = xi * A
+    logsky = HT(logsky_h)
+    sky = ift.PointwiseExponential(logsky)
+
+    M = ift.DiagonalOperator(exposure)
+    GR = ift.GeometryRemover(position_space)
+    # Set up instrumental response
+    R = GR * M
+
+    # Generate mock data
+    d_space = R.target[0]
+    lamb = R(sky)
+    mock_position = ift.from_random('normal', lamb.position.domain)
+    data = lamb.at(mock_position).value
+    data = np.random.poisson(data.to_global_data().astype(np.float64))
+    data = ift.Field.from_global_data(d_space, data)
+
+    # Compute likelihood and Hamiltonian
+    position = ift.from_random('normal', lamb.position.domain)
+    likelihood = ift.PoissonianEnergy(lamb, data)
+    ic_cg = ift.GradientNormController(iteration_limit=50)
+    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
+                                           tol_abs_gradnorm=1e-3)
+    minimizer = ift.RelaxedNewton(ic_newton)
+
+    # Minimize the Hamiltonian
+    H = ift.Hamiltonian(likelihood, ic_cg)
+    H, convergence = minimizer(H)
+
+    # Plot results
+    result_sky = sky.at(H.position).value
+    # PLOTTING

demos/getting_started_3.py

+import nifty5 as ift
+from nifty5.library.los_response import LOSResponse
+from nifty5.library.amplitude_model import make_amplitude_model
+from nifty5.library.smooth_sky import make_correlated_field
+import numpy as np
+from scipy.io import loadmat
+
+
+def get_random_LOS(n_los):
+    starts = list(np.random.uniform(0, 1, (n_los, 2)).T)
+    ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
+
+    return starts, ends
+
+
+if __name__ == '__main__':
+    # ## ABOUT THIS TUTORIAL
+    np.random.seed(42)
+    position_space = ift.RGSpace([128, 128])
+
+    # Setting up an amplitude model
+    A, amplitude_internals = make_amplitude_model(
+        position_space, 16, 1, 10, -4., 1, 0., 1.)
+
+    # Building the model for a correlated signal
+    harmonic_space = position_space.get_default_codomain()
+    ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
+    power_space = A.value.domain[0]
+    power_distributor = ift.PowerDistributor(harmonic_space, power_space)
+    position = {}
+    position['xi'] = ift.Field.from_random('normal', harmonic_space)
+    position = ift.MultiField(position)
+
+    xi = ift.Variable(position)['xi']
+    Amp = power_distributor(A)
+    correlated_field_h = Amp * xi
+    correlated_field = ht(correlated_field_h)
+    # # alternatively to the block above one can do:
+    # correlated_field, _ = make_correlated_field(position_space, A)
+
+    # apply some nonlinearity
+    signal = ift.PointwisePositiveTanh(correlated_field)
+    # Building the Line of Sight response
+    LOS_starts, LOS_ends = get_random_LOS(100)
+    R = LOSResponse(position_space, starts=LOS_starts, ends=LOS_ends)
+    # build signal response model and model likelihood
+    signal_response = R(signal)
+    # specify noise
+    data_space = R.target
+    noise = .001
+    N = ift.ScalingOperator(noise, data_space)
+
+    # generate mock data
+    MOCK_POSITION = ift.from_random('normal', signal.position.domain)
+    data = signal_response.at(MOCK_POSITION).value + N.draw_sample()
+
+    # set up model likelihood
+    likelihood = ift.GaussianEnergy(signal_response, mean=data, covariance=N)
+
+    # set up minimization and inversion schemes
+    ic_cg = ift.GradientNormController(iteration_limit=10)
+    ic_sampling = ift.GradientNormController(iteration_limit=100)
+    ic_newton = ift.GradientNormController(name='Newton', iteration_limit=100)
+    minimizer = ift.RelaxedNewton(ic_newton)
+
+    # build model Hamiltonian
+    H = ift.Hamiltonian(likelihood, ic_cg,
+                        iteration_controller_sampling=ic_sampling)
+
+    INITIAL_POSITION = ift.from_random('normal', H.position.domain)
+    position = INITIAL_POSITION
+
+    ift.plot(signal.at(MOCK_POSITION).value, name='truth.pdf')
+    ift.plot(R.adjoint_times(data), name='data.pdf')
+    ift.plot([A.at(MOCK_POSITION).value], name='power.pdf')
+
+    # number of samples used to estimate the KL
+    N_samples = 20
+    for i in range(5):
+        H = H.at(position)
+        samples = [H.curvature.draw_sample(from_inverse=True)
+                   for _ in range(N_samples)]
+
+        KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg)
+        KL, convergence = minimizer(KL)
+        position = KL.position
+
+        ift.plot(signal.at(position).value, name='reconstruction.pdf')
+
+        ift.plot([A.at(position).value, A.at(MOCK_POSITION).value],
+                 name='power.pdf')
+
+    avrg = 0.
+    va = 0.
+    powers = []
+    for sample in samples:
+        sam = signal.at(sample + position).value
+        powers.append(A.at(sample+position).value)
+        avrg += sam
+        va += sam**2
+
+    avrg /= len(samples)
+    va /= len(samples)
+    va -= avrg**2
+    std = ift.sqrt(va)
+    ift.plot(avrg, name='avrg.pdf')
+    ift.plot(std, name='std.pdf')
+    ift.plot([A.at(position).value, A.at(MOCK_POSITION).value]+powers,
+             name='power.pdf')

nifty5/domains/dof_space.py

@@ -21,7 +21,17 @@ from .structured_domain import StructuredDomain
 
 
 class DOFSpace(StructuredDomain):
-    """Generic degree-of-freedom space."""
+    """Generic degree-of-freedom space. It is defined as the domain of some
+    DOFDistributor.
+    Its entries represent the underlying degrees of freedom of some other
+    space, according to the dofdex.
+
+    Parameters
+    ----------
+    dof_weights: 1-D numpy array
+        A numpy array containing the multiplicity of each individual degree of
+        freedom.
+    """
 
     _needed_for_hash = ["_dvol"]
 

nifty5/energies/kl.py

@@ -19,7 +19,6 @@ class SampledKullbachLeiblerDivergence(Energy):
     def at(self, position):
         return self.__class__(self._h.at(position), self._res_samples)
 
-
     @property
     @memo
     def value(self):

nifty5/library/__init__.py

@@ -4,6 +4,5 @@ from .gaussian_energy import GaussianEnergy
 from .los_response import LOSResponse
 from .point_sources import PointSources
 from .poissonian_energy import PoissonianEnergy
-from .smooth_sky import make_smooth_mf_sky_model, make_smooth_sky_model
 from .wiener_filter_curvature import WienerFilterCurvature
 from .wiener_filter_energy import WienerFilterEnergy

nifty5/library/amplitude_model.py

@@ -49,7 +49,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
     phi_sig = np.array([sv, iv])
 
     slope = SlopeOperator(param_space, logk_space, phi_sig)
-    norm_phi_mean = Field(param_space, val=phi_mean/phi_sig)
+    norm_phi_mean = Field.from_global_data(param_space, phi_mean/phi_sig)
 
     fields = {keys[0]: Field.from_random('normal', dof_space),
               keys[1]: Field.from_random('normal', param_space)}
@@ -94,7 +94,7 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
     # Prepare q_array
     q_array = np.zeros((dim,) + shape)
     if dim == 1:
-        ks = domain.get_k_length_array().val
+        ks = domain.get_k_length_array().to_global_data()
         q_array = np.array([ks])
     else:
         for i in range(dim):
@@ -121,4 +121,4 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
         # Do summation
         cepstrum_field[sl2] = np.sum(cepstrum_field[sl], axis=i)
 
-    return Field(domain, val=cepstrum_field)
+    return Field.from_global_data(domain, cepstrum_field)

nifty5/library/point_sources.py

@@ -21,40 +21,40 @@ class PointSources(Model):
     @property
     @memo
     def value(self):
-        points = self.position['points'].to_global_data()
+        points = self.position['points'].local_data
         points = np.clip(points, None, 8.2)
-        points = Field(self.position['points'].domain, points)
+        points = Field.from_local_data(self.position['points'].domain, points)
         return self.IG(points, self._alpha, self._q)
 
     @property
     @memo
     def gradient(self):
-        u = self.position['points']
-        inner = norm.pdf(u.val)
-        outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u.val),
+        u = self.position['points'].local_data
+        inner = norm.pdf(u)
+        outer_inv = invgamma.pdf(invgamma.ppf(norm.cdf(u),
                                               self._alpha,
                                               scale=self._q),
                                  self._alpha, scale=self._q)
         # FIXME
         outer_inv = np.clip(outer_inv, 1e-20, None)
         outer = 1/outer_inv
-        grad = Field(u.domain, val=inner*outer)
+        grad = Field.from_local_data(u.domain, inner*outer)
         grad = makeOp(MultiField({'points': grad}))
         return SelectionOperator(grad.target, 'points')*grad
 
     @staticmethod
     def IG(field, alpha, q):
-        foo = invgamma.ppf(norm.cdf(field.val), alpha, scale=q)
-        return Field(field.domain, val=foo)
+        foo = invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q)
+        return Field.from_local_data(field.domain, foo)
 
     @staticmethod
     def IG_prime(field, alpha, q):
-        inner = norm.pdf(field.val)
-        outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.val), alpha, scale=q), alpha, scale=q)
+        inner = norm.pdf(field.local_data)
+        outer = invgamma.pdf(invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q), alpha, scale=q)
         # # FIXME
         # outer = np.clip(outer, 1e-20, None)
         outer = 1/outer
-        return Field(field.domain, val=inner*outer)
+        return Field.from_local_data(field.domain, inner*outer)
 
     @staticmethod
     def inverseIG(u, alpha, q):

nifty5/library/smooth_sky.py

-def make_smooth_sky_model(s_space, amplitude_model):
+def make_correlated_field(s_space, amplitude_model):
     '''
     Method for construction of correlated sky model
 
@@ -22,12 +22,12 @@ def make_smooth_sky_model(s_space, amplitude_model):
 
     xi = Variable(position)['xi']
     A = power_distributor(amplitude_model)
-    logsky_h = A * xi
-    logsky = ht(logsky_h)
-    internals = {'logsky_h': logsky_h,
+    correlated_field_h = A * xi
+    correlated_field = ht(correlated_field_h)
+    internals = {'correlated_field_h': correlated_field_h,
                  'power_distributor': power_distributor,
                  'ht': ht}
-    return PointwiseExponential(logsky), internals
+    return correlated_field, internals
 
 
 def make_smooth_mf_sky_model(s_space_spatial, s_space_energy,

nifty5/models/binary_helpers.py

@@ -21,9 +21,9 @@ from ..sugar import makeOp
 from .model import Model
 
 
-def _joint_position(op1, op2):
-    a = op1.position._val
-    b = op2.position._val
+def _joint_position(model1, model2):
+    a = model1.position._val
+    b = model2.position._val
     # Note: In python >3.5 one could do {**a, **b}
     ab = a.copy()
     ab.update(b)
@@ -31,57 +31,78 @@ def _joint_position(op1, op2):
 
 
 class ScalarMul(Model):
-    def __init__(self, factor, op):
-        # TODO op -> model
-        super(ScalarMul, self).__init__(op.position)
+    """Class representing a model multiplied by a scalar factor."""
+    def __init__(self, factor, model):
+        super(ScalarMul, self).__init__(model.position)
         # TODO -> floating
         if not isinstance(factor, (float, int)):
             raise TypeError
 
-        self._op = op
+        self._model = model
         self._factor = factor
 
-        self._value = self._factor * self._op.value
-        self._gradient = self._factor * self._op.gradient
+        self._value = self._factor * self._model.value
+        self._gradient = self._factor * self._model.gradient
 
     def at(self, position):
-        return self.__class__(self._factor, self._op.at(position))
+        return self.__class__(self._factor, self._model.at(position))
 
 
 class Add(Model):
-    def __init__(self, position, op1, op2):
+    """Class representing the sum of two models."""
+    def __init__(self, position, model1, model2):
         super(Add, self).__init__(position)
 
-        self._op1 = op1.at(position)
-        self._op2 = op2.at(position)
+        self._model1 = model1.at(position)
+        self._model2 = model2.at(position)
 
-        self._value = self._op1.value + self._op2.value
-        self._gradient = self._op1.gradient + self._op2.gradient
+        self._value = self._model1.value + self._model2.value
+        self._gradient = self._model1.gradient + self._model2.gradient
 
     @staticmethod
-    def make(op1, op2):
-        position = _joint_position(op1, op2)
-        return Add(position, op1, op2)
+    def make(model1, model2):
+        """Build the sum of two models.
+
+        Parameters
+        ----------
+        model1: Model
+            First model.
+        model2: Model
+            Second model
+        """