There is a maintenance of MPCDF Gitlab on Thursday, April 22st 2020, 9:00 am CEST - Expect some service interruptions during this time

Commit f98653ce authored by Martin Reinecke's avatar Martin Reinecke

stage2

parent eda70a19
...@@ -167,7 +167,7 @@ class RGSpace(StructuredDomain): ...@@ -167,7 +167,7 @@ class RGSpace(StructuredDomain):
raise NotImplementedError raise NotImplementedError
op = HarmonicTransformOperator(self, self.get_default_codomain()) op = HarmonicTransformOperator(self, self.get_default_codomain())
dist = op.target[0]._get_dist_array() dist = op.target[0]._get_dist_array()
kernel = Field.from_local_data(op.target, func(dist.local_data)) kernel = Field(op.target, func(dist.val))
kernel = kernel / kernel.integrate() kernel = kernel / kernel.integrate()
return op.adjoint_times(kernel.weight(1)) return op.adjoint_times(kernel.weight(1))
......
...@@ -32,8 +32,7 @@ __all__ = ["consistency_check", "check_jacobian_consistency", ...@@ -32,8 +32,7 @@ __all__ = ["consistency_check", "check_jacobian_consistency",
def assert_allclose(f1, f2, atol, rtol): def assert_allclose(f1, f2, atol, rtol):
if isinstance(f1, Field): if isinstance(f1, Field):
return np.testing.assert_allclose(f1.local_data, f2.local_data, return np.testing.assert_allclose(f1.val, f2.val, atol=atol, rtol=rtol)
atol=atol, rtol=rtol)
for key, val in f1.items(): for key, val in f1.items():
assert_allclose(val, f2[key], atol=atol, rtol=rtol) assert_allclose(val, f2[key], atol=atol, rtol=rtol)
......
...@@ -67,7 +67,7 @@ class InverseGammaOperator(Operator): ...@@ -67,7 +67,7 @@ class InverseGammaOperator(Operator):
def apply(self, x): def apply(self, x):
self._check_input(x) self._check_input(x)
lin = isinstance(x, Linearization) lin = isinstance(x, Linearization)
val = x.val.local_data if lin else x.local_data val = x.val.val if lin else x.val
val = (np.clip(val, self._xmin, self._xmax) - self._xmin) / self._delta val = (np.clip(val, self._xmin, self._xmax) - self._xmin) / self._delta
...@@ -76,26 +76,26 @@ class InverseGammaOperator(Operator): ...@@ -76,26 +76,26 @@ class InverseGammaOperator(Operator):
w = val - fi w = val - fi
res = np.exp((1 - w)*self._table[fi] + w*self._table[fi + 1]) res = np.exp((1 - w)*self._table[fi] + w*self._table[fi + 1])
points = Field.from_local_data(self._domain, res) points = Field(self._domain, res)
if not lin: if not lin:
return points return points
# Derivative of linear interpolation # Derivative of linear interpolation
der = self._deriv[fi]*res der = self._deriv[fi]*res
jac = makeOp(Field.from_local_data(self._domain, der)) jac = makeOp(Field(self._domain, der))
jac = jac(x.jac) jac = jac(x.jac)
return x.new(points, jac) return x.new(points, jac)
@staticmethod @staticmethod
def IG(field, alpha, q): def IG(field, alpha, q):
foo = invgamma.ppf(norm.cdf(field.local_data), alpha, scale=q) foo = invgamma.ppf(norm.cdf(field.val), alpha, scale=q)
return Field.from_local_data(field.domain, foo) return Field(field.domain, foo)
@staticmethod @staticmethod
def inverseIG(u, alpha, q): def inverseIG(u, alpha, q):
res = norm.ppf(invgamma.cdf(u.local_data, alpha, scale=q)) res = norm.ppf(invgamma.cdf(u.val, alpha, scale=q))
return Field.from_local_data(u.domain, res) return Field(u.domain, res)
@property @property
def alpha(self): def alpha(self):
......
...@@ -227,9 +227,9 @@ class LOSResponse(LinearOperator): ...@@ -227,9 +227,9 @@ class LOSResponse(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
if mode == self.TIMES: if mode == self.TIMES:
result_arr = self._smat.matvec(x.local_data.reshape(-1)) result_arr = self._smat.matvec(x.val.reshape(-1))
return Field.from_global_data(self._target, result_arr, return Field.from_global_data(self._target, result_arr,
sum_up=True) sum_up=True)
local_input_data = x.to_global_data().reshape(-1) local_input_data = x.to_global_data().reshape(-1)
res = self._smat.rmatvec(local_input_data).reshape(self._local_shape) res = self._smat.rmatvec(local_input_data).reshape(self._local_shape)
return Field.from_local_data(self._domain, res) return Field(self._domain, res)
...@@ -320,10 +320,10 @@ class Linearization(object): ...@@ -320,10 +320,10 @@ class Linearization(object):
def sinc(self): def sinc(self):
tmp = self._val.sinc() tmp = self._val.sinc()
tmp2 = ((np.pi*self._val).cos()-tmp)/self._val tmp2 = ((np.pi*self._val).cos()-tmp)/self._val
ind = self._val.local_data == 0 ind = self._val.val == 0
loc = tmp2.local_data.copy() loc = tmp2.val.copy()
loc[ind] = 0 loc[ind] = 0
tmp2 = Field.from_local_data(tmp.domain, loc) tmp2 = Field(tmp.domain, loc)
return self.new(tmp, makeOp(tmp2)(self._jac)) return self.new(tmp, makeOp(tmp2)(self._jac))
def log(self): def log(self):
...@@ -370,10 +370,10 @@ class Linearization(object): ...@@ -370,10 +370,10 @@ class Linearization(object):
tmp = self._val.absolute() tmp = self._val.absolute()
tmp2 = self._val.sign() tmp2 = self._val.sign()
ind = self._val.local_data == 0 ind = self._val.val == 0
loc = tmp2.local_data.copy().astype(float) loc = tmp2.val.copy().astype(float)
loc[ind] = np.nan loc[ind] = np.nan
tmp2 = Field.from_local_data(tmp.domain, loc) tmp2 = Field(tmp.domain, loc)
return self.new(tmp, makeOp(tmp2)(self._jac)) return self.new(tmp, makeOp(tmp2)(self._jac))
......
...@@ -47,7 +47,7 @@ class EnergyAdapter(Energy): ...@@ -47,7 +47,7 @@ class EnergyAdapter(Energy):
self._want_metric = want_metric self._want_metric = want_metric
lin = Linearization.make_partial_var(position, constants, want_metric) lin = Linearization.make_partial_var(position, constants, want_metric)
tmp = self._op(lin) tmp = self._op(lin)
self._val = tmp.val.local_data[()] self._val = tmp.val.val[()]
self._grad = tmp.gradient self._grad = tmp.gradient
self._metric = tmp._metric self._metric = tmp._metric
......
...@@ -111,10 +111,10 @@ class MetricGaussianKL(Energy): ...@@ -111,10 +111,10 @@ class MetricGaussianKL(Energy):
for s in self._samples: for s in self._samples:
tmp = self._hamiltonian(self._lin+s) tmp = self._hamiltonian(self._lin+s)
if v is None: if v is None:
v = tmp.val.local_data[()] v = tmp.val.val[()]
g = tmp.gradient g = tmp.gradient
else: else:
v += tmp.val.local_data[()] v += tmp.val.val[()]
g = g + tmp.gradient g = g + tmp.gradient
self._val = v / len(self._samples) self._val = v / len(self._samples)
self._grad = g * (1./len(self._samples)) self._grad = g * (1./len(self._samples))
......
...@@ -51,10 +51,9 @@ def np_allreduce_sum(arr): ...@@ -51,10 +51,9 @@ def np_allreduce_sum(arr):
def allreduce_sum_field(fld): def allreduce_sum_field(fld):
if isinstance(fld, Field): if isinstance(fld, Field):
return Field.from_local_data(fld.domain, return Field(fld.domain, np_allreduce_sum(fld.val))
np_allreduce_sum(fld.local_data))
res = tuple( res = tuple(
Field.from_local_data(f.domain, np_allreduce_sum(f.local_data)) Field(f.domain, np_allreduce_sum(f.val))
for f in fld.values()) for f in fld.values())
return MultiField(fld.domain, res) return MultiField(fld.domain, res)
...@@ -183,16 +182,16 @@ class MetricGaussianKL_MPI(Energy): ...@@ -183,16 +182,16 @@ class MetricGaussianKL_MPI(Energy):
v, g = None, None v, g = None, None
if len(self._samples) == 0: # hack if there are too many MPI tasks if len(self._samples) == 0: # hack if there are too many MPI tasks
tmp = self._hamiltonian(self._lin) tmp = self._hamiltonian(self._lin)
v = 0. * tmp.val.local_data v = 0. * tmp.val.val
g = 0. * tmp.gradient g = 0. * tmp.gradient
else: else:
for s in self._samples: for s in self._samples:
tmp = self._hamiltonian(self._lin+s) tmp = self._hamiltonian(self._lin+s)
if v is None: if v is None:
v = tmp.val.local_data.copy() v = tmp.val.val.copy()
g = tmp.gradient g = tmp.gradient
else: else:
v += tmp.val.local_data v += tmp.val.val
g = g + tmp.gradient g = g + tmp.gradient
self._val = np_allreduce_sum(v)[()] / self._n_samples self._val = np_allreduce_sum(v)[()] / self._n_samples
self._grad = allreduce_sum_field(g) / self._n_samples self._grad = allreduce_sum_field(g) / self._n_samples
......
...@@ -32,7 +32,7 @@ def _multiToArray(fld): ...@@ -32,7 +32,7 @@ def _multiToArray(fld):
ofs = 0 ofs = 0
for val in fld.values(): for val in fld.values():
sz2 = 2*val.size if iscomplextype(val.dtype) else val.size sz2 = 2*val.size if iscomplextype(val.dtype) else val.size
locdat = val.local_data.reshape(-1) locdat = val.val.reshape(-1)
if iscomplextype(val.dtype): if iscomplextype(val.dtype):
locdat = locdat.view(locdat.real.dtype) locdat = locdat.view(locdat.real.dtype)
res[ofs:ofs+sz2] = locdat res[ofs:ofs+sz2] = locdat
...@@ -42,20 +42,19 @@ def _multiToArray(fld): ...@@ -42,20 +42,19 @@ def _multiToArray(fld):
def _toArray(fld): def _toArray(fld):
if isinstance(fld, Field): if isinstance(fld, Field):
return fld.local_data.reshape(-1) return fld.val.reshape(-1)
return _multiToArray(fld) return _multiToArray(fld)
def _toArray_rw(fld): def _toArray_rw(fld):
if isinstance(fld, Field): if isinstance(fld, Field):
return fld.local_data.copy().reshape(-1) return fld.val.copy().reshape(-1)
return _multiToArray(fld) return _multiToArray(fld)
def _toField(arr, template): def _toField(arr, template):
if isinstance(template, Field): if isinstance(template, Field):
return Field.from_local_data(template.domain, return Field(template.domain, arr.reshape(template.shape).copy())
arr.reshape(template.shape).copy())
ofs = 0 ofs = 0
res = [] res = []
for v in template.values(): for v in template.values():
...@@ -63,7 +62,7 @@ def _toField(arr, template): ...@@ -63,7 +62,7 @@ def _toField(arr, template):
locdat = arr[ofs:ofs+sz2].copy() locdat = arr[ofs:ofs+sz2].copy()
if iscomplextype(v.dtype): if iscomplextype(v.dtype):
locdat = locdat.view(np.complex128) locdat = locdat.view(np.complex128)
res.append(Field.from_local_data(v.domain, locdat.reshape(v.shape))) res.append(Field(v.domain, locdat.reshape(v.shape)))
ofs += sz2 ofs += sz2
return MultiField(template.domain, tuple(res)) return MultiField(template.domain, tuple(res))
......
...@@ -73,9 +73,9 @@ class _DomRemover(LinearOperator): ...@@ -73,9 +73,9 @@ class _DomRemover(LinearOperator):
@staticmethod @staticmethod
def _check_float_dtype(fld): def _check_float_dtype(fld):
if isinstance(fld, MultiField): if isinstance(fld, MultiField):
dts = [ff.local_data.dtype for ff in fld.values()] dts = [ff.dtype for ff in fld.values()]
elif isinstance(fld, Field): elif isinstance(fld, Field):
dts = [fld.local_data.dtype] dts = [fld.dtype]
else: else:
raise TypeError raise TypeError
for dt in dts: for dt in dts:
......
...@@ -54,13 +54,13 @@ class ContractionOperator(LinearOperator): ...@@ -54,13 +54,13 @@ class ContractionOperator(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
if mode == self.ADJOINT_TIMES: if mode == self.ADJOINT_TIMES:
ldat = x.to_global_data() if 0 in self._spaces else x.local_data ldat = x.val
shp = [] shp = []
for i, dom in enumerate(self._domain): for i, dom in enumerate(self._domain):
tmp = dom.shape if i > 0 else dom.local_shape tmp = dom.shape
shp += tmp if i not in self._spaces else (1,)*len(dom.shape) shp += tmp if i not in self._spaces else (1,)*len(dom.shape)
ldat = np.broadcast_to(ldat.reshape(shp), self._domain.local_shape) ldat = np.broadcast_to(ldat.reshape(shp), self._domain.shape)
res = Field.from_local_data(self._domain, ldat) res = Field(self._domain, ldat)
if self._weight != 0: if self._weight != 0:
res = res.weight(self._weight, spaces=self._spaces) res = res.weight(self._weight, spaces=self._spaces)
return res return res
......
...@@ -130,15 +130,15 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -130,15 +130,15 @@ class DiagonalOperator(EndomorphicOperator):
self._check_input(x, mode) self._check_input(x, mode)
# shortcut for most common cases # shortcut for most common cases
if mode == 1 or (not self._complex and mode == 2): if mode == 1 or (not self._complex and mode == 2):
return Field.from_local_data(x.domain, x.local_data*self._ldiag) return Field(x.domain, x.val*self._ldiag)
xdiag = self._ldiag xdiag = self._ldiag
if self._complex and (mode & 10): # adjoint or inverse adjoint if self._complex and (mode & 10): # adjoint or inverse adjoint
xdiag = xdiag.conj() xdiag = xdiag.conj()
if mode & 3: if mode & 3:
return Field.from_local_data(x.domain, x.local_data*xdiag) return Field(x.domain, x.val*xdiag)
return Field.from_local_data(x.domain, x.local_data/xdiag) return Field(x.domain, x.val/xdiag)
def _flip_modes(self, trafo): def _flip_modes(self, trafo):
if trafo == self.ADJOINT_BIT and not self._complex: # shortcut if trafo == self.ADJOINT_BIT and not self._complex: # shortcut
......
...@@ -64,18 +64,18 @@ class DOFDistributor(LinearOperator): ...@@ -64,18 +64,18 @@ class DOFDistributor(LinearOperator):
if partner != dofdex.domain[0]: if partner != dofdex.domain[0]:
raise ValueError("incorrect dofdex domain") raise ValueError("incorrect dofdex domain")
ldat = dofdex.local_data ldat = dofdex.val
if ldat.size == 0: # can happen for weird configurations if ldat.size == 0: # can happen for weird configurations
nbin = 0 nbin = 0
else: else:
nbin = ldat.max() nbin = ldat.max()
nbin = nbin + 1 nbin = nbin + 1
if partner.scalar_dvol is not None: if partner.scalar_dvol is not None:
wgt = np.bincount(dofdex.local_data.ravel(), minlength=nbin) wgt = np.bincount(dofdex.val.ravel(), minlength=nbin)
wgt = wgt*partner.scalar_dvol wgt = wgt*partner.scalar_dvol
else: else:
dvol = Field.from_global_data(partner, partner.dvol).local_data dvol = Field.from_global_data(partner, partner.dvol).val
wgt = np.bincount(dofdex.local_data.ravel(), wgt = np.bincount(dofdex.val.ravel(),
minlength=nbin, weights=dvol) minlength=nbin, weights=dvol)
# The explicit conversion to float64 is necessary because bincount # The explicit conversion to float64 is necessary because bincount
# sometimes returns its result as an integer array, even when # sometimes returns its result as an integer array, even when
......
...@@ -184,7 +184,7 @@ class PoissonianEnergy(EnergyOperator): ...@@ -184,7 +184,7 @@ class PoissonianEnergy(EnergyOperator):
def __init__(self, d): def __init__(self, d):
if not isinstance(d, Field) or not np.issubdtype(d.dtype, np.integer): if not isinstance(d, Field) or not np.issubdtype(d.dtype, np.integer):
raise TypeError raise TypeError
if np.any(d.local_data < 0): if np.any(d.val < 0):
raise ValueError raise ValueError
self._d = d self._d = d
self._domain = DomainTuple.make(d.domain) self._domain = DomainTuple.make(d.domain)
...@@ -192,7 +192,7 @@ class PoissonianEnergy(EnergyOperator): ...@@ -192,7 +192,7 @@ class PoissonianEnergy(EnergyOperator):
def apply(self, x): def apply(self, x):
self._check_input(x) self._check_input(x)
res = x.sum() res = x.sum()
tmp = res.val.local_data if isinstance(res, Linearization) else res tmp = res.val.val if isinstance(res, Linearization) else res
# if we have no infinity here, we can continue with the calculation; # if we have no infinity here, we can continue with the calculation;
# otherwise we know that the result must also be infinity # otherwise we know that the result must also be infinity
if not np.isinf(tmp): if not np.isinf(tmp):
...@@ -231,8 +231,7 @@ class InverseGammaLikelihood(EnergyOperator): ...@@ -231,8 +231,7 @@ class InverseGammaLikelihood(EnergyOperator):
raise TypeError raise TypeError
self._beta = beta self._beta = beta
if np.isscalar(alpha): if np.isscalar(alpha):
alpha = Field.from_local_data( alpha = Field(beta.domain, np.full(beta.shape, alpha))
beta.domain, np.full(beta.local_data.shape, alpha))
elif not isinstance(alpha, Field): elif not isinstance(alpha, Field):
raise TypeError raise TypeError
self._alphap1 = alpha+1 self._alphap1 = alpha+1
...@@ -302,7 +301,7 @@ class BernoulliEnergy(EnergyOperator): ...@@ -302,7 +301,7 @@ class BernoulliEnergy(EnergyOperator):
def __init__(self, d): def __init__(self, d):
if not isinstance(d, Field) or not np.issubdtype(d.dtype, np.integer): if not isinstance(d, Field) or not np.issubdtype(d.dtype, np.integer):
raise TypeError raise TypeError
if not np.all(np.logical_or(d.local_data == 0, d.local_data == 1)): if not np.all(np.logical_or(d.val == 0, d.val == 1)):
raise ValueError raise ValueError
self._d = d self._d = d
self._domain = DomainTuple.make(d.domain) self._domain = DomainTuple.make(d.domain)
......
...@@ -42,7 +42,7 @@ class VdotOperator(LinearOperator): ...@@ -42,7 +42,7 @@ class VdotOperator(LinearOperator):
self._check_mode(mode) self._check_mode(mode)
if mode == self.TIMES: if mode == self.TIMES:
return Field.scalar(self._field.vdot(x)) return Field.scalar(self._field.vdot(x))
return self._field*x.local_data[()] return self._field*x.val[()]
class ConjugationOperator(EndomorphicOperator): class ConjugationOperator(EndomorphicOperator):
......
...@@ -51,10 +51,10 @@ class ValueInserter(LinearOperator): ...@@ -51,10 +51,10 @@ class ValueInserter(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
x = x.to_global_data() x = x.val
if mode == self.TIMES: if mode == self.TIMES:
res = np.zeros(self.target.shape, dtype=x.dtype) res = np.zeros(self.target.shape, dtype=x.dtype)
res[self._index] = x res[self._index] = x
return Field.from_global_data(self._tgt(mode), res) return Field(self._tgt(mode), res)
else: else:
return Field.scalar(x[self._index]) return Field.scalar(x[self._index])
...@@ -56,7 +56,7 @@ def test_inverse_gamma(field): ...@@ -56,7 +56,7 @@ def test_inverse_gamma(field):
field = field.exp() field = field.exp()
space = field.domain space = field.domain
d = np.random.normal(10, size=space.shape)**2 d = np.random.normal(10, size=space.shape)**2
d = ift.Field.from_global_data(space, d) d = ift.Field(space, d)
energy = ift.InverseGammaLikelihood(d) energy = ift.InverseGammaLikelihood(d)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-5) ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
...@@ -65,7 +65,7 @@ def testPoissonian(field): ...@@ -65,7 +65,7 @@ def testPoissonian(field):
field = field.exp() field = field.exp()
space = field.domain space = field.domain
d = np.random.poisson(120, size=space.shape) d = np.random.poisson(120, size=space.shape)
d = ift.Field.from_global_data(space, d) d = ift.Field(space, d)
energy = ift.PoissonianEnergy(d) energy = ift.PoissonianEnergy(d)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-7) ift.extra.check_jacobian_consistency(energy, field, tol=1e-7)
...@@ -86,6 +86,6 @@ def test_bernoulli(field): ...@@ -86,6 +86,6 @@ def test_bernoulli(field):
field = field.sigmoid() field = field.sigmoid()
space = field.domain space = field.domain
d = np.random.binomial(1, 0.1, size=space.shape) d = np.random.binomial(1, 0.1, size=space.shape)
d = ift.Field.from_global_data(space, d) d = ift.Field(space, d)
energy = ift.BernoulliEnergy(d) energy = ift.BernoulliEnergy(d)
ift.extra.check_jacobian_consistency(energy, field, tol=1e-5) ift.extra.check_jacobian_consistency(energy, field, tol=1e-5)
...@@ -68,8 +68,8 @@ def test_quadratic_minimization(minimizer, space): ...@@ -68,8 +68,8 @@ def test_quadratic_minimization(minimizer, space):
assert_equal(convergence, IC.CONVERGED) assert_equal(convergence, IC.CONVERGED)
assert_allclose( assert_allclose(
energy.position.local_data, energy.position.val,
1./covariance_diagonal.local_data, 1./covariance_diagonal.val,
rtol=1e-3, rtol=1e-3,
atol=1e-3) atol=1e-3)
...@@ -90,8 +90,8 @@ def test_WF_curvature(space): ...@@ -90,8 +90,8 @@ def test_WF_curvature(space):
iteration_controller_sampling=IC) iteration_controller_sampling=IC)
m = curv.inverse(required_result) m = curv.inverse(required_result)
assert_allclose( assert_allclose(
m.local_data, m.val,
1./all_diag.local_data, 1./all_diag.val,
rtol=1e-3, rtol=1e-3,
atol=1e-3) atol=1e-3)
curv.draw_sample() curv.draw_sample()
...@@ -107,8 +107,8 @@ def test_WF_curvature(space): ...@@ -107,8 +107,8 @@ def test_WF_curvature(space):
iteration_controller_sampling=IC) iteration_controller_sampling=IC)
m = curv.inverse(required_result) m = curv.inverse(required_result)
assert_allclose( assert_allclose(
m.local_data, m.val,
1./all_diag.local_data, 1./all_diag.val,
rtol=1e-3,