Commit 549bfa3d authored by Martin Reinecke's avatar Martin Reinecke

add Jakob's changes

parent c2697357
Pipeline #22361 failed with stage
in 4 minutes and 7 seconds
import nifty2go as ift
from nifty2go.library import NonlinearWienerFilterEnergy, NonlinearPowerEnergy
from nifty2go.library.nonlinearities import *
import numpy as np
np.random.seed(42)
def adjust_zero_mode(m0, t0):
zero_position = len(m0.shape)*(0,)
zero_mode = m0.val[zero_position]
m0.val[zero_position] = zero_mode / abs(zero_mode)
t0.val[0] += 2 * np.log(abs(zero_mode))
return m0, t0
if __name__ == "__main__":
noise_level = 1.
p_spec = (lambda k: (1. / (k + 1) ** 2))
# nonlinearity = Linear()
nonlinearity = Exponential()
# Set up position space
# s_space = ift.RGSpace([1024])
s_space = ift.HPSpace(32)
# Define harmonic transformation and associated harmonic space
FFT = ift.FFTOperator(s_space)
h_space = FFT.target[0]
# Setting up power space
p_space = ift.PowerSpace(h_space,
binbounds=ift.PowerSpace.useful_binbounds(
h_space, logarithmic=True))
s_spec = ift.Field(p_space,val=1.)
# Choosing the prior correlation structure and defining correlation operator
p = ift.Field(p_space,val = p_spec(p_space.k_lengths))
log_p = ift.log(p)
S = ift.create_power_operator(h_space, power_spectrum=s_spec)
# Drawing a sample sh from the prior distribution in harmonic space
sp = ift.Field(p_space, val=s_spec)
sh = ift.power_synthesize(sp)
# Choosing the measurement instrument
# Instrument = SmoothingOperator(s_space, sigma=0.01)
mask = np.ones(s_space.shape)
mask[6000:8000] = 0.
mask = ift.Field(s_space,val=mask)
MaskOperator = ift.DiagonalOperator(mask)
InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
MeasurementOperator = ift.ComposedOperator([MaskOperator, InstrumentResponse])
d_space = MeasurementOperator.target
noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
N = ift.DiagonalOperator(noise_covariance)
n = ift.Field.from_random(domain=d_space,
random_type='normal',
std=noise_level,
mean=0)
Projection = ift.PowerProjectionOperator(domain= h_space, power_space=p_space)
power = Projection.adjoint_times(ift.exp(0.5 * log_p))
# Creating the mock data
true_sky = nonlinearity(FFT.adjoint_times(power * sh))
d = MeasurementOperator(true_sky) + n
flat_power = ift.Field(p_space,val=10e-8)
m0 = ift.power_synthesize(flat_power)
t0 = ift.Field(p_space,val=-4.)
power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
IC1 = ift.GradientNormController(verbose=True, iteration_limit=100,
tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(IC1)
ICI = ift.GradientNormController(verbose=False, name="ICI",
iteration_limit=500,
tol_abs_gradnorm=1e-3)
inverter = ift.ConjugateGradient(controller=ICI)
for i in range(20):
power0 = Projection.adjoint_times(ift.exp(0.5*t0))
map0_energy = NonlinearWienerFilterEnergy(m0, d, MeasurementOperator,
nonlinearity, FFT, power0, N, S, inverter=inverter)
# Minimization with chosen minimizer
map0_energy, convergence = minimizer(map0_energy)
m0 = map0_energy.position
# Updating parameters for correlation structure reconstruction
D0 = map0_energy.curvature
# Initializing the power energy with updated parameters
power0_energy = NonlinearPowerEnergy(position=t0, d=d, N=N, m=m0,
D=D0, FFT=FFT, Instrument=MeasurementOperator,
nonlinearity=nonlinearity,Projection=Projection,
sigma=1., samples=2, inverter=inverter)
(power0_energy, convergence) = minimizer(power0_energy)
# Setting new power spectrum
t0_ = t0.copy()
t0 = power0_energy.position
# break degeneracy between amplitude and excitation by setting excitation monopole to 1
m0, t0 = adjust_zero_mode(m0,t0)
ift.plotting.plot(true_sky)
ift.plotting.plot(nonlinearity(FFT.adjoint_times(power0*m0)),title='reconstructed_sky')
ift.plotting.plot(MeasurementOperator.adjoint_times(d))
......@@ -137,14 +137,15 @@ class Field(object):
raise TypeError("could not infer domain from value")
return DomainTuple.make(domain)
# MR: this needs some rethinking ... do we need to have at least float64?
@staticmethod
def _infer_dtype(dtype, val):
if val is None or dtype is not None:
return np.result_type(dtype, np.float64)
if dtype is not None:
return dtype
if val is None:
raise ValueError("could not infer dtype")
if isinstance(val, Field):
return val.dtype
return np.result_type(val, np.float64)
return np.result_type(val)
@staticmethod
def from_random(random_type, domain, dtype=np.float64, **kwargs):
......
from __future__ import division
from .linear_operator import LinearOperator
from .diagonal_operator import DiagonalOperator
from .endomorphic_operator import EndomorphicOperator
from .fft_smoothing_operator import FFTSmoothingOperator
from .fft_operator import FFTOperator
from .inversion_enabler import InversionEnabler
from .composed_operator import ComposedOperator
from .response_operator import ResponseOperator
from .laplace_operator import LaplaceOperator
from .smoothness_operator import SmoothnessOperator
from .power_projection_operator import PowerProjectionOperator
from .dof_projection_operator import DOFProjectionOperator
from .linear_operator import LinearOperator
from .. import Field
from ..spaces import DOFSpace
class DOFProjectionOperator(LinearOperator):
def __init__(self, domain, dofdex, space=None):
super(DOFProjectionOperator, self).__init__()
self._domain = DomainTuple.make(domain)
if space is None and len(self._domain) == 1:
space = 0
space = int(space)
if space < 0 or space >= len(self.domain):
raise ValueError("space index out of range")
partner = self._domain[space]
if not isinstance(dofdex, Field):
raise TypeError("dofdex must be a Field")
if not isinstance(dofdex.dtype, np.integer):
raise TypeError("dofdex must contain integer numbers")
if partner != dofdex.domain:
raise ValueError("incorrect dofdex domain")
nbin = dofdex.max()
if partner.scalar_dvol() is not None:
wgt = np.bincount(dobj.local_data(dofdex.val).ravel(),
minlength=nbin)
wgt *= partner.scalar_dvol()
else:
dvol = dobj.local_data(partner.dvol())
wgt = np.bincount(dobj.local_data(dofdex.val).ravel(),
minlength=nbin, weights=dvol)
# The explicit conversion to float64 is necessary because bincount
# sometimes returns its result as an integer array, even when
# floating-point weights are present ...
wgt = wgt.astype(np.float64, copy=False)
wgt = dobj.np_allreduce_sum(wgt)
if (wgt == 0).any():
raise ValueError("empty bins detected")
self._space = space
tgt = list(self._domain)
tgt[self._space] = DOFSpace(wgt)
self._target = DomainTuple.make(tgt)
if dobj.default_distaxis() in self.domain.axes[self._space]:
dofdex = dobj.local_data(dofdex)
else: # dofdex must be available fully on every task
dofdex = dobj.to_global_data(dofdex)
self._dofdex = dofdex.ravel()
firstaxis = self._domain.axes[self._space][0]
lastaxis = self._domain.axes[self._space][-1]
arrshape = dobj.local_shape(self._domain.shape, 0)
presize = np.prod(arrshape[0:firstaxis], dtype=np.int)
postsize = np.prod(arrshape[lastaxis+1:], dtype=np.int)
self._hshape = (presize, self._target[self._space].shape[0], postsize)
self._pshape = (presize, self._dofdex.size, postsize)
def _times(self, x):
arr = dobj.local_data(x.weight(1).val)
arr = arr.reshape(self._pshape)
oarr = np.zeros(self._hshape, dtype=x.dtype)
np.add.at(oarr, (slice(None), self._dofdex, slice(None)), arr)
if dobj.distaxis(x.val) in x.domain.axes[self._space]:
oarr = dobj.np_allreduce_sum(oarr).reshape(self._target.shape)
res = Field(self._target, dobj.from_global_data(oarr))
else:
oarr = oarr.reshape(dobj.local_shape(self._target.shape,
dobj.distaxis(x.val)))
res = Field(self._target,
dobj.from_local_data(self._target.shape, oarr,
dobj.default_distaxis()))
return res.weight(-1, spaces=self._space)
def _adjoint_times(self, x):
res = Field.empty(self._domain, dtype=x.dtype)
if dobj.distaxis(x.val) in x.domain.axes[self._space]:
arr = dobj.to_global_data(x.val)
else:
arr = dobj.local_data(x.val)
arr = arr.reshape(self._hshape)
oarr = dobj.local_data(res.val).reshape(self._pshape)
oarr[()] = arr[(slice(None), self._dofdex, slice(None))]
return res
@property
def domain(self):
return self._domain
@property
def target(self):
return self._target
@property
def unitary(self):
return False
......@@ -2,6 +2,7 @@ from .gl_space import GLSpace
from .hp_space import HPSpace
from .lm_space import LMSpace
from .power_space import PowerSpace
from .dof_space import DOFSpace
from .rg_space import RGSpace
from .space import Space
from .field_array import FieldArray
from .space import Space
import numpy as np
class DOFSpace(Space):
def __init__(self, dof_weights):
super(DOFSpace, self).__init__()
self._dvol = tuple(dof_weights)
self._needed_for_hash +=['_dvol']
@property
def harmonic(self):
return False
@property
def shape(self):
return (len(self._dvol),)
@property
def dim(self):
return len(self._dvol)
def scalar_dvol(self):
return None
def dvol(self):
return self._dvol
def __repr__(self):
return 'this is a dof space'
......@@ -155,7 +155,7 @@ class PowerSpace(Space):
# floating-point weights are present ...
temp_k_lengths = np.bincount(dobj.local_data(temp_pindex).ravel(),
weights=dobj.local_data(k_length_array.val).ravel(),
minlength=nbin).astype(np.float64)
minlength=nbin).astype(np.float64, copy=False)
temp_k_lengths = dobj.np_allreduce_sum(temp_k_lengths) / temp_rho
temp_dvol = temp_rho*pdvol
self._powerIndexCache[key] = (binbounds,
......
......@@ -37,7 +37,7 @@ class Test_Interface(unittest.TestCase):
def test_return_types(self, domain, attribute_desired_type):
attribute = attribute_desired_type[0]
desired_type = attribute_desired_type[1]
f = ift.Field(domain=domain)
f = ift.Field(domain=domain, val=1.)
assert_equal(isinstance(getattr(f, attribute), desired_type), True)
......
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