Commit 92638ed6 authored by Philipp Arras's avatar Philipp Arras
Browse files

Remove old files

parent 346cb2a2
Pipeline #63019 failed with stages
in 4 minutes and 31 seconds
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import nifty5 as ift
if __name__ == '__main__':
import numpy as np
np.random.seed(42)
nspecs = 2
a = []
spaces = []
for _ in range(nspecs):
ndim = 2
sspace = ift.RGSpace(
np.linspace(16, 20, num=ndim).astype(np.int),
np.linspace(2.3, 7.99, num=ndim))
hspace = sspace.get_default_codomain()
spaces.append(sspace)
target = ift.PowerSpace(hspace)
a.append(ift.NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1))
tgt = ift.makeDomain(tuple(spaces))
op = ift.CorrelatedFieldNormAmplitude(tgt, a, 0., 1.)
fld = ift.from_random('normal', op.domain)
ift.extra.check_jacobian_consistency(op, fld)
import numpy as np
import nifty5 as ift
if __name__ == '__main__':
np.random.seed(42)
ndim = 1
sspace = ift.RGSpace(
np.linspace(16, 20, num=ndim).astype(np.int),
np.linspace(2.3, 7.99, num=ndim))
hspace = sspace.get_default_codomain()
target = ift.PowerSpace(hspace)
op = ift.NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1)
if not isinstance(op.target[0], ift.PowerSpace):
raise RuntimeError
fld = ift.from_random('normal', op.domain)
ift.extra.check_jacobian_consistency(op, fld)
pd = ift.PowerDistributor(hspace, power_space=target)
ht = ift.HartleyOperator(pd.target)
if np.testing.assert_allclose((pd @ op)(fld).integrate(), 1):
raise RuntimeError
if np.testing.assert_allclose(
(ht @ pd @ op)(fld).to_global_data()[ndim*(0,)], 1):
raise RuntimeError
import nifty5 as ift
import numpy as np
def mfplot(fld, A1, A2, name):
p = ift.Plot()
dom = fld.domain
mi, ma = np.min(fld.val), np.max(fld.val)
for ii in range(dom.shape[-1]):
extr = ift.DomainTupleFieldInserter(op.target, 1, (ii,)).adjoint
p.add(extr(fld), zmin=mi, zmax=ma)
p.add(A1)
p.add(A2)
p.add(fld)
p.output(name=name, xsize=14, ysize=14)
if __name__ == '__main__':
np.random.seed(42)
sspace = ift.RGSpace((512, 512))
hspace = sspace.get_default_codomain()
target = ift.PowerSpace(hspace,
ift.PowerSpace.useful_binbounds(hspace, True))
# A = Amplitude(target, [5, -3, 1, 0], [1, 1, 1, 0.5],
# ['rest', 'smooth', 'wienersigma'])
# fld = ift.from_random('normal', A.domain)
# ift.extra.check_jacobian_consistency(A, fld)
# op = CorrelatedFieldNormAmplitude(sspace, A, 3, 2)
# plts = []
# p = ift.Plot()
# for _ in range(25):
# fld = ift.from_random('normal', op.domain)
# plts.append(A.force(fld))
# p.add(op(fld))
# p.output(name='debug1.png', xsize=20, ysize=20)
# ift.single_plot(plts, name='debug.png')
# Multifrequency
sspace1 = ift.RGSpace((512, 512))
hspace1 = sspace1.get_default_codomain()
target1 = ift.PowerSpace(hspace1,
ift.PowerSpace.useful_binbounds(hspace1, True))
A1 = ift.WPAmplitude(target1, [0, -2, 0], [1E-5, 1, 1], 1, 0.99,
['rest1', 'smooth1', 'wienersigma1'])
sspace2 = ift.RGSpace((20,), distances=(1e7))
hspace2 = sspace2.get_default_codomain()
target2 = ift.PowerSpace(hspace2,
ift.PowerSpace.useful_binbounds(hspace2, True))
A2 = ift.WPAmplitude(target2, [0, -2, 0], [1E-5, 1, 1], 1, 0.99,
['rest2', 'smooth2', 'wienersigma2'])
op = ift.CorrelatedFieldNormAmplitude((sspace1, sspace2), (A1, A2), 3, 2)
for jj in range(10):
fld = ift.from_random('normal', op.domain)
mfplot(op(fld), A1.force(fld), A2.force(fld), 'debug{}.png'.format(jj))
...@@ -76,9 +76,6 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL ...@@ -76,9 +76,6 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL
from .sugar import * from .sugar import *
from .plot import Plot from .plot import Plot
from .library.smooth_linear_amplitude import (
SLAmplitude, LinearSLAmplitude, CepstrumOperator, WPAmplitude)
from .library.normalized_amplitude import NormalizedAmplitude
from .library.inverse_gamma_operator import InverseGammaOperator from .library.inverse_gamma_operator import InverseGammaOperator
from .library.los_response import LOSResponse from .library.los_response import LOSResponse
from .library.dynamic_operator import (dynamic_operator, from .library.dynamic_operator import (dynamic_operator,
...@@ -86,11 +83,10 @@ from .library.dynamic_operator import (dynamic_operator, ...@@ -86,11 +83,10 @@ from .library.dynamic_operator import (dynamic_operator,
from .library.light_cone_operator import LightConeOperator from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import CorrelatedField, MfCorrelatedField, CorrelatedFieldNormAmplitude
from .library.adjust_variances import (make_adjust_variances_hamiltonian, from .library.adjust_variances import (make_adjust_variances_hamiltonian,
do_adjust_variances) do_adjust_variances)
from .library.gridder import GridderMaker from .library.gridder import GridderMaker
from .library.final_amplitude import FinalAmplitude from .library.correlated_fields import CorrelatedFieldMaker
from . import extra from . import extra
......
...@@ -12,222 +12,281 @@ ...@@ -12,222 +12,281 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# #
# Copyright(C) 2013-2019 Max-Planck-Society # Copyright(C) 2013-2019 Max-Planck-Society
# Authors: Philipp Frank, Philipp Arras
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce import numpy as np
from operator import mul from numpy.testing import assert_allclose
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..extra import check_jacobian_consistency
from ..field import Field from ..field import Field
from ..multi_domain import MultiDomain
from ..operators.adder import Adder from ..operators.adder import Adder
from ..operators.contraction_operator import ContractionOperator from ..operators.contraction_operator import ContractionOperator
from ..operators.distributors import PowerDistributor from ..operators.distributors import PowerDistributor
from ..operators.endomorphic_operator import EndomorphicOperator
from ..operators.harmonic_operators import HarmonicTransformOperator from ..operators.harmonic_operators import HarmonicTransformOperator
from ..operators.linear_operator import LinearOperator
from ..operators.operator import Operator from ..operators.operator import Operator
from ..operators.simple_linear_operators import VdotOperator, ducktape from ..operators.simple_linear_operators import VdotOperator, ducktape
from ..sugar import full, makeOp from ..operators.value_inserter import ValueInserter
from ..sugar import from_global_data, from_random, full, makeDomain
def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
"""Constructs an operator which turns a white Gaussian excitation field def _lognormal_moment_matching(mean, sig, key):
into a correlated field. mean, sig, key = float(mean), float(sig), str(key)
assert sig > 0
This function returns an operator which implements: logsig = np.sqrt(np.log((sig/mean)**2 + 1))
logmean = np.log(mean) - logsig**2/2
ht @ (vol * A * xi), return _normal(logmean, logsig, key).exp()
where `ht` is a harmonic transform operator, `A` is the square root of the
prior covariance and `xi` is the excitation field. def _normal(mean, sig, key):
return Adder(Field.scalar(mean)) @ (
Parameters sig*ducktape(DomainTuple.scalar_domain(), None, key))
----------
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly one space. class _SlopeOperator(Operator):
amplitude_operator: Operator def __init__(self, smooth, loglogavgslope):
name : string self._domain = MultiDomain.union(
:class:`MultiField` key for the xi-field. [smooth.domain, loglogavgslope.domain])
codomain : Domain self._target = smooth.target
The codomain for target[0]. If not supplied, it is inferred. self._smooth = smooth
self._llas = loglogavgslope
Returns logkl = _log_k_lengths(self._target[0])
------- assert logkl.shape[0] == self._target[0].shape[0] - 1
Operator logkl -= logkl[0]
Correlated field logkl = np.insert(logkl, 0, 0)
self._t = VdotOperator(from_global_data(self._target, logkl)).adjoint
Notes self._T = float(logkl[-1])
----- ind = (smooth.target.shape[0] - 1,)
In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore self._extr_op = ValueInserter(smooth.target, ind).adjoint
the operator constructed by this method will output a correlated field
with *periodic* boundary conditions. If a non-periodic field is needed, def apply(self, x):
one needs to combine this operator with a :class:`FieldZeroPadder`. self._check_input(x)
""" smooth = self._smooth(x.extract(self._smooth.domain))
tgt = DomainTuple.make(target) res0 = self._llas(x.extract(self._llas.domain))
if len(tgt) > 1: res1 = self._extr_op(smooth)/self._T
raise ValueError return self._t(res0 - res1) + smooth
if codomain is None:
codomain = tgt[0].get_default_codomain()
h_space = codomain def _log_k_lengths(pspace):
ht = HarmonicTransformOperator(h_space, target=tgt[0]) return np.log(pspace.k_lengths[1:])
if isinstance(amplitude_operator, Operator):
p_space = amplitude_operator.target[0]
elif isinstance(amplitude_operator, Field): class _TwoLogIntegrations(LinearOperator):
p_space = amplitude_operator.domain[0] def __init__(self, target):
else: self._target = makeDomain(target)
raise TypeError self._domain = makeDomain(
power_distributor = PowerDistributor(h_space, p_space) UnstructuredDomain((2, self.target.shape[0] - 2)))
A = power_distributor(amplitude_operator) self._capability = self.TIMES | self.ADJOINT_TIMES
vol = h_space.scalar_dvol**-0.5 if not isinstance(self._target[0], PowerSpace):
xi = ducktape(h_space, None, name) raise TypeError
# When doubling the resolution of `tgt` the value of the highest k-mode logk_lengths = _log_k_lengths(self._target[0])
# will scale with a square root. `vol` cancels this effect such that the self._logvol = logk_lengths[1:] - logk_lengths[:-1]
# same power spectrum can be used for the spaces with the same volume,
# different resolutions and the same object in them. def apply(self, x, mode):
if isinstance(amplitude_operator, Field): self._check_input(x, mode)
return ht(makeOp(A) @ xi).scale(vol) if mode == self.TIMES:
return ht(A*xi).scale(vol) x = x.to_global_data()
res = np.empty(self._target.shape)
res[0] = 0
def MfCorrelatedField(target, amplitudes, name='xi'): res[1] = 0
"""Constructs an operator which turns white Gaussian excitation fields res[2:] = np.cumsum(x[1])
into a correlated field defined on a DomainTuple with two entries and two res[2:] = (res[2:] + res[1:-1])/2*self._logvol + x[0]
separate correlation structures. res[2:] = np.cumsum(res[2:])
return from_global_data(self._target, res)
This operator may be used as a model for multi-frequency reconstructions else:
with a correlation structure in both spatial and energy direction. x = x.to_global_data_rw()
res = np.zeros(self._domain.shape)
Parameters x[2:] = np.cumsum(x[2:][::-1])[::-1]
---------- res[0] += x[2:]
target : Domain, DomainTuple or tuple of Domain x[2:] *= self._logvol/2.
Target of the operator. Must contain exactly two spaces. res[1:-1] += res[2:]
amplitudes: iterable of Operator x[1] += np.cumsum(res[2:][::-1])[::-1]
List of two amplitude operators. return from_global_data(self._domain, res)
name : string
:class:`MultiField` key for xi-field.
class _Normalization(Operator):
Returns def __init__(self, domain):
------- self._domain = self._target = makeDomain(domain)
Operator hspace = self._domain[0].harmonic_partner
Correlated field pd = PowerDistributor(hspace, power_space=self._domain[0])
# TODO Does not work on sphere yet
Notes cst = pd.adjoint(full(pd.target, 1.)).to_global_data_rw()
----- cst[0] = 0
In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore self._cst = from_global_data(self._domain, cst)
the operator constructed by this method will output a correlated field self._specsum = _SpecialSum(self._domain)
with *periodic* boundary conditions. If a non-periodic field is needed,
one needs to combine this operator with a :class:`FieldZeroPadder` or even def apply(self, x):
two (one for the energy and one for the spatial subdomain) self._check_input(x)
""" amp = x.exp()
tgt = DomainTuple.make(target) spec = (2*x).exp()
if len(tgt) != 2: # FIXME This normalizes also the zeromode which is supposed to be left
raise ValueError # untouched by this operator
if len(amplitudes) != 2: return self._specsum(self._cst*spec)**(-0.5)*amp
raise ValueError
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt]) class _SpecialSum(EndomorphicOperator):
ht1 = HarmonicTransformOperator(hsp, target=tgt[0], space=0) def __init__(self, domain):
ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1) self._domain = makeDomain(domain)
ht = ht2 @ ht1 self._capability = self.TIMES | self.ADJOINT_TIMES
psp = [aa.target[0] for aa in amplitudes] def apply(self, x, mode):
pd0 = PowerDistributor(hsp, psp[0], 0) self._check_input(x, mode)
pd1 = PowerDistributor(pd0.domain, psp[1], 1) return full(self._tgt(mode), x.sum())
pd = pd0 @ pd1
dd0 = ContractionOperator(pd.domain, 1).adjoint class CorrelatedFieldMaker:
dd1 = ContractionOperator(pd.domain, 0).adjoint def __init__(self):
d = [dd0, dd1] self._amplitudes = []
a = [dd @ amplitudes[ii] for ii, dd in enumerate(d)] def add_fluctuations_from_ops(self, target, fluctuations, flexibility,
a = reduce(mul, a) asperity, loglogavgslope, key):
A = pd @ a """
# For `vol` see comment in `CorrelatedField` fluctuations > 0
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp]) flexibility > 0
return ht(vol*A*ducktape(hsp, None, name)) asperity > 0
loglogavgslope probably negative
"""
def CorrelatedFieldNormAmplitude(target, assert isinstance(fluctuations, Operator)
amplitudes, assert isinstance(flexibility, Operator)
stdmean, assert isinstance(asperity, Operator)
stdstd, assert isinstance(loglogavgslope, Operator)
names=['xi', 'std']): target = makeDomain(target)
"""Constructs an operator which turns white Gaussian excitation fields assert len(target) == 1
into a correlated field defined on a DomainTuple with n entries and n assert isinstance(target[0], PowerSpace)
separate correlation structures.
twolog = _TwoLogIntegrations(target)
This operator may be used as a model for multi-frequency reconstructions dt = twolog._logvol
with a correlation structure in both spatial and energy direction. sc = np.zeros(twolog.domain.shape)
sc[0] = sc[1] = np.sqrt(dt)
Parameters sc = from_global_data(twolog.domain, sc)
---------- expander = VdotOperator(sc).adjoint
target : Domain, DomainTuple or tuple of Domain sigmasq = expander @ flexibility
Target of the operator. Must contain exactly n spaces.
amplitudes: Opertor, iterable of Operator dist = np.zeros(twolog.domain.shape)
Amplitude operator if n = 1 or list of n amplitude operators. dist[0] += 1.
stdmean : float dist = from_global_data(twolog.domain, dist)
Prior mean of the overall standart deviation. scale = VdotOperator(dist).adjoint @ asperity
stdstd : float
Prior standart deviation of the overall standart deviation. shift = np.ones(scale.target.shape)
names : iterable of string shift[0] = dt**2/12.
:class:`MultiField` keys for xi-field and std-field. shift = from_global_data(scale.target, shift)
scale = sigmasq*(Adder(shift) @ scale).sqrt()
Returns
------- smooth = twolog @ (scale*ducktape(scale.target, None, key))
Operator smoothslope = _SlopeOperator(smooth, loglogavgslope)
Correlated field
# move to tests
Notes assert_allclose(
----- smooth(from_random('normal', smooth.domain)).val[0:2], 0)
In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore check_jacobian_consistency(smooth, from_random('normal',
the operator constructed by this method will output a correlated field smooth.domain))
with *periodic* boundary conditions. If a non-periodic field is needed, # end move to tests
one needs to combine this operator with a :class:`FieldZeroPadder` or even
two (one for the energy and one for the spatial subdomain) normal_ampl = _Normalization(target) @ smoothslope
""" vol = target[0].harmonic_partner.get_default_codomain().total_volume
arr = np.zeros(target.shape)
amps = [ arr[1:] = vol
amplitudes, expander = VdotOperator(from_global_data(target, arr)).adjoint
] if isinstance(amplitudes, (Operator, Field)) else amplitudes mask = np.zeros(target.shape)
mask[0] = vol
cls = Operator if isinstance(amps[0], Operator) else Field adder = Adder(from_global_data(target, mask))
assert all([isinstance(aa, cls) for aa in amps]) ampl = adder @ ((expander @ fluctuations)*normal_ampl)
tgt = DomainTuple.make(target) # Move to tests
if len(tgt) != len(amps): # FIXME This test fails but it is not relevant for the final result
raise ValueError # assert_allclose(
stdmean, stdstd = float(stdmean), float(stdstd) # normal_ampl(from_random('normal', normal_ampl.domain)).val[0], 1)
if stdstd <= 0: assert_allclose(ampl(from_random('normal', ampl.domain)).val[0], vol)
raise ValueError # End move to tests
if cls == Field: self._amplitudes.append(ampl)
psp = [aa.domain[0] for aa in amps]
else: def add_fluctuations(self, target, fluctuations_mean, fluctuations_stddev,
psp = [aa.target[0] for aa in amps] flexibility_mean, flexibility_stddev, asperity_mean,
hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt]) asperity_stddev, loglogavgslope_mean,
loglogavgslope_stddev, prefix):
ht = HarmonicTransformOperator(hsp, target=tgt[0], space=0) fluctuations_mean = float(fluctuations_mean)
pd = PowerDistributor(hsp, psp[0], 0) fluctuations_stddev = float(fluctuations_stddev)
flexibility_mean = float(flexibility_mean)
for i in range(1, len(amps)): flexibility_stddev = float(flexibility_stddev)
ht = HarmonicTransformOperator(ht.target, target=tgt[i], space=i) @ ht asperity_mean = float(asperity_mean)
pd = pd @ PowerDistributor(pd.domain, psp[i], space=i) asperity_stddev = float(asperity_stddev)
loglogavgslope_mean = float(loglogavgslope_mean)
spaces = tuple(range(len(amps))) loglogavgslope_stddev = float(loglogavgslope_stddev)
prefix = str(prefix)
a = ContractionOperator(pd.domain, spaces[1:]).adjoint(amps[0]) assert fluctuations_stddev > 0
for i in range(1, len(amps)): assert fluctuations_mean > 0
a = a*(ContractionOperator( assert flexibility_stddev > 0
pd.domain, spaces[:i] + spaces[(i + 1):]).adjoint(amps[i])) assert flexibility_mean > 0
assert asperity_stddev > 0
expander = VdotOperator(full(pd.domain, 1.)).adjoint assert asperity_mean > 0
assert loglogavgslope_stddev > 0
Std = stdstd*ducktape(expander.domain, None, names[1])
Std = expander @ (Adder(full(expander.domain, stdmean)) @ Std).exp() fluct = _lognormal_moment_matching(fluctuations_mean,
fluctuations_stddev,
if cls == Field: prefix + 'fluctuations')
A = pd @ (makeOp(a) @ Std) flex = _lognormal_moment_matching(flexibility_mean, flexibility_stddev,
else: prefix + 'flexibility')
A = pd @ (Std*a) asp = _lognormal_moment_matching(asperity_mean, asperity_stddev,
# For `vol` see comment in `CorrelatedField` prefix + 'asperity')
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp]) avgsl = _normal(loglogavgslope_mean, loglogavgslope_mean,
return ht(vol*A*ducktape(hsp, None, names[0])) prefix + 'loglogavgslope')
self.add_fluctuations_from_ops(target, fluct, flex, asp, avgsl,
prefix + 'spectrum')
def finalize_from_op(self, zeromode):
raise NotImplementedError
def finalize(self,
offset_amplitude_mean,
offset_amplitude_stddev,
prefix,
offset=None):
"""
offset vs zeromode: volume factor
"""
offset_amplitude_stddev = float(offset_amplitude_stddev)
offset_amplitude_mean = float(offset_amplitude_mean)
assert offset_amplitude_stddev > 0
assert offset_amplitude_mean > 0
if offset is not None:
offset = float(offset)
hspace = makeDomain(
[dd.target