Commit 988aa85a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_6' into timing_helpers

parents 781be410 8d416d27
Pipeline #64608 passed with stages
in 9 minutes and 51 seconds
...@@ -44,19 +44,19 @@ def _reshaper(x, N): ...@@ -44,19 +44,19 @@ def _reshaper(x, N):
raise TypeError("Shape of parameters cannot be interpreted") raise TypeError("Shape of parameters cannot be interpreted")
def _lognormal_moments(mean, sig, N = 0): def _lognormal_moments(mean, sig, N=0):
if N == 0: if N == 0:
mean, sig = np.asfarray(mean), np.asfarray(sig) mean, sig = np.asfarray(mean), np.asfarray(sig)
else: else:
mean, sig = (_reshaper(param, N) for param in (mean, sig)) mean, sig = (_reshaper(param, N) for param in (mean, sig))
assert np.all(mean > 0 ) assert np.all(mean > 0)
assert np.all(sig > 0) assert np.all(sig > 0)
logsig = np.sqrt(np.log((sig/mean)**2 + 1)) logsig = np.sqrt(np.log((sig/mean)**2 + 1))
logmean = np.log(mean) - logsig**2/2 logmean = np.log(mean) - logsig**2/2
return logmean, logsig return logmean, logsig
def _normal(mean, sig, key, N = 0): def _normal(mean, sig, key, N=0):
if N == 0: if N == 0:
domain = DomainTuple.scalar_domain() domain = DomainTuple.scalar_domain()
mean, sig = np.asfarray(mean), np.asfarray(sig) mean, sig = np.asfarray(mean), np.asfarray(sig)
...@@ -64,7 +64,7 @@ def _normal(mean, sig, key, N = 0): ...@@ -64,7 +64,7 @@ def _normal(mean, sig, key, N = 0):
domain = UnstructuredDomain(N) domain = UnstructuredDomain(N)
mean, sig = (_reshaper(param, N) for param in (mean, sig)) mean, sig = (_reshaper(param, N) for param in (mean, sig))
return Adder(from_global_data(domain, mean)) @ ( return Adder(from_global_data(domain, mean)) @ (
DiagonalOperator(from_global_data(domain,sig)) DiagonalOperator(from_global_data(domain, sig))
@ ducktape(domain, None, key)) @ ducktape(domain, None, key))
...@@ -126,7 +126,7 @@ class _LognormalMomentMatching(Operator): ...@@ -126,7 +126,7 @@ class _LognormalMomentMatching(Operator):
class _SlopeRemover(EndomorphicOperator): class _SlopeRemover(EndomorphicOperator):
def __init__(self, domain, space = 0): def __init__(self, domain, space=0):
self._domain = makeDomain(domain) self._domain = makeDomain(domain)
assert isinstance(self._domain[space], PowerSpace) assert isinstance(self._domain[space], PowerSpace)
logkl = _relative_log_k_lengths(self._domain[space]) logkl = _relative_log_k_lengths(self._domain[space])
...@@ -146,12 +146,12 @@ class _SlopeRemover(EndomorphicOperator): ...@@ -146,12 +146,12 @@ class _SlopeRemover(EndomorphicOperator):
else: else:
res = x.copy() res = x.copy()
res[self._last] -= (x*self._sc[self._extender]).sum( res[self._last] -= (x*self._sc[self._extender]).sum(
axis = self._space, keepdims = True) axis=self._space, keepdims=True)
return from_global_data(self._tgt(mode), res) return from_global_data(self._tgt(mode), res)
class _TwoLogIntegrations(LinearOperator): class _TwoLogIntegrations(LinearOperator):
def __init__(self, target, space = 0): def __init__(self, target, space=0):
self._target = makeDomain(target) self._target = makeDomain(target)
assert isinstance(self.target[space], PowerSpace) assert isinstance(self.target[space], PowerSpace)
dom = list(self._target) dom = list(self._target)
...@@ -164,42 +164,42 @@ class _TwoLogIntegrations(LinearOperator): ...@@ -164,42 +164,42 @@ class _TwoLogIntegrations(LinearOperator):
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
#Maybe make class properties # Maybe make class properties
axis = self._target.axes[self._space][0] axis = self._target.axes[self._space][0]
sl = (slice(None),)*axis sl = (slice(None),)*axis
extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis) extender_sl = (None,)*axis + (slice(None),) + (None,)*(self._target.axes[-1][-1] - axis)
first = sl + (0,) first = sl + (0,)
second = sl + (1,) second = sl + (1,)
from_third = sl + (slice(2,None),) from_third = sl + (slice(2, None),)
no_border = sl + (slice(1,-1),) no_border = sl + (slice(1, -1),)
reverse = sl + (slice(None,None,-1),) reverse = sl + (slice(None, None, -1),)
if mode == self.TIMES: if mode == self.TIMES:
x = x.to_global_data() x = x.to_global_data()
res = np.empty(self._target.shape) res = np.empty(self._target.shape)
res[first] = res[second] = 0 res[first] = res[second] = 0
res[from_third] = np.cumsum(x[second], axis = axis) res[from_third] = np.cumsum(x[second], axis=axis)
res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first] res[from_third] = (res[from_third] + res[no_border])/2*self._log_vol[extender_sl] + x[first]
res[from_third] = np.cumsum(res[from_third], axis = axis) res[from_third] = np.cumsum(res[from_third], axis=axis)
else: else:
x = x.to_global_data_rw() x = x.to_global_data_rw()
res = np.zeros(self._domain.shape) res = np.zeros(self._domain.shape)
x[from_third] = np.cumsum(x[from_third][reverse], axis = axis)[reverse] x[from_third] = np.cumsum(x[from_third][reverse], axis=axis)[reverse]
res[first] += x[from_third] res[first] += x[from_third]
x[from_third] *= (self._log_vol/2.)[extender_sl] x[from_third] *= (self._log_vol/2.)[extender_sl]
x[no_border] += x[from_third] x[no_border] += x[from_third]
res[second] += np.cumsum(x[from_third][reverse], axis = axis)[reverse] res[second] += np.cumsum(x[from_third][reverse], axis=axis)[reverse]
return from_global_data(self._tgt(mode), res) return from_global_data(self._tgt(mode), res)
class _Normalization(Operator): class _Normalization(Operator):
def __init__(self, domain, space = 0): def __init__(self, domain, space=0):
self._domain = self._target = makeDomain(domain) self._domain = self._target = makeDomain(domain)
assert isinstance(self._domain[space], PowerSpace) assert isinstance(self._domain[space], PowerSpace)
hspace = list(self._domain) hspace = list(self._domain)
hspace[space] = hspace[space].harmonic_partner hspace[space] = hspace[space].harmonic_partner
hspace = makeDomain(hspace) hspace = makeDomain(hspace)
pd = PowerDistributor(hspace, power_space=self._domain[space], space = space) pd = PowerDistributor(hspace, power_space=self._domain[space], space=space)
mode_multiplicity = pd.adjoint(full(pd.target, 1.)).to_global_data_rw() mode_multiplicity = pd.adjoint(full(pd.target, 1.)).to_global_data_rw()
zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,) zero_mode = (slice(None),)*self._domain.axes[space][0] + (0,)
mode_multiplicity[zero_mode] = 0 mode_multiplicity[zero_mode] = 0
...@@ -216,7 +216,7 @@ class _Normalization(Operator): ...@@ -216,7 +216,7 @@ class _Normalization(Operator):
class _SpecialSum(EndomorphicOperator): class _SpecialSum(EndomorphicOperator):
def __init__(self, domain, space = 0): def __init__(self, domain, space=0):
self._domain = makeDomain(domain) self._domain = makeDomain(domain)
self._capability = self.TIMES | self.ADJOINT_TIMES self._capability = self.TIMES | self.ADJOINT_TIMES
self._contractor = ContractionOperator(domain, space) self._contractor = ContractionOperator(domain, space)
...@@ -227,7 +227,7 @@ class _SpecialSum(EndomorphicOperator): ...@@ -227,7 +227,7 @@ class _SpecialSum(EndomorphicOperator):
class _Distributor(LinearOperator): class _Distributor(LinearOperator):
def __init__(self, dofdex, domain, target, space = 0): def __init__(self, dofdex, domain, target, space=0):
self._dofdex = dofdex self._dofdex = dofdex
self._target = makeDomain(target) self._target = makeDomain(target)
...@@ -244,7 +244,7 @@ class _Distributor(LinearOperator): ...@@ -244,7 +244,7 @@ class _Distributor(LinearOperator):
res = np.empty(self._tgt(mode).shape) res = np.empty(self._tgt(mode).shape)
res[self._dofdex] = x res[self._dofdex] = x
return from_global_data(self._tgt(mode), res) return from_global_data(self._tgt(mode), res)
class _Amplitude(Operator): class _Amplitude(Operator):
def __init__(self, target, fluctuations, flexibility, asperity, def __init__(self, target, fluctuations, flexibility, asperity,
...@@ -271,16 +271,15 @@ class _Amplitude(Operator): ...@@ -271,16 +271,15 @@ class _Amplitude(Operator):
N_copies = 0 N_copies = 0
space = 0 space = 0
distributed_tgt = target = makeDomain(target) distributed_tgt = target = makeDomain(target)
azm_expander = ContractionOperator(distributed_tgt, spaces = space).adjoint azm_expander = ContractionOperator(distributed_tgt, spaces=space).adjoint
assert isinstance(target[space], PowerSpace) assert isinstance(target[space], PowerSpace)
twolog = _TwoLogIntegrations(target, space) twolog = _TwoLogIntegrations(target, space)
dom = twolog.domain dom = twolog.domain
shp = dom[space].shape shp = dom[space].shape
expander = ContractionOperator(dom, spaces = space).adjoint expander = ContractionOperator(dom, spaces=space).adjoint
ps_expander = ContractionOperator(twolog.target, spaces = space).adjoint ps_expander = ContractionOperator(twolog.target, spaces=space).adjoint
# Prepare constant fields # Prepare constant fields
foo = np.zeros(shp) foo = np.zeros(shp)
...@@ -294,17 +293,18 @@ class _Amplitude(Operator): ...@@ -294,17 +293,18 @@ class _Amplitude(Operator):
foo = np.ones(shp) foo = np.ones(shp)
foo[0] = _log_vol(target[space])**2/12. foo[0] = _log_vol(target[space])**2/12.
shift = DiagonalOperator(from_global_data(dom[space], foo), dom, space) shift = DiagonalOperator(from_global_data(dom[space], foo), dom, space)
vslope = DiagonalOperator( vslope = DiagonalOperator(
from_global_data(target[space], _relative_log_k_lengths(target[space])), from_global_data(target[space],
target, space) _relative_log_k_lengths(target[space])),
target, space)
foo, bar = [np.zeros(target[space].shape) for _ in range(2)] foo, bar = [np.zeros(target[space].shape) for _ in range(2)]
bar[1:] = foo[0] = totvol bar[1:] = foo[0] = totvol
vol0, vol1 = [DiagonalOperator(from_global_data(target[space], aa), vol0, vol1 = [DiagonalOperator(from_global_data(target[space], aa),
target, space) for aa in (foo, bar)] target, space) for aa in (foo, bar)]
#Prepare fields for Adder # Prepare fields for Adder
shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)] shift, vol0 = [op(full(op.domain, 1)) for op in (shift, vol0)]
# End prepare constant fields # End prepare constant fields
...@@ -322,7 +322,7 @@ class _Amplitude(Operator): ...@@ -322,7 +322,7 @@ class _Amplitude(Operator):
op = Distributor @ op op = Distributor @ op
sig_fluc = Distributor @ sig_fluc sig_fluc = Distributor @ sig_fluc
op = Adder(Distributor(vol0)) @ (sig_fluc*(azm_expander @ azm.one_over())*op) op = Adder(Distributor(vol0)) @ (sig_fluc*(azm_expander @ azm.one_over())*op)
self._fluc = (_Distributor(dofdex, fluctuations.target, distributed_tgt[0]) @ self._fluc = (_Distributor(dofdex, fluctuations.target, distributed_tgt[0]) @
fluctuations) fluctuations)
else: else:
op = Adder(vol0) @ (sig_fluc*(azm_expander @ azm.one_over())*op) op = Adder(vol0) @ (sig_fluc*(azm_expander @ azm.one_over())*op)
...@@ -350,8 +350,8 @@ class CorrelatedFieldMaker: ...@@ -350,8 +350,8 @@ class CorrelatedFieldMaker:
@staticmethod @staticmethod
def make(offset_amplitude_mean, offset_amplitude_stddev, prefix, def make(offset_amplitude_mean, offset_amplitude_stddev, prefix,
total_N = 0, total_N=0,
dofdex = None): dofdex=None):
if dofdex is None: if dofdex is None:
dofdex = np.full(total_N, 0) dofdex = np.full(total_N, 0)
else: else:
...@@ -362,7 +362,7 @@ class CorrelatedFieldMaker: ...@@ -362,7 +362,7 @@ class CorrelatedFieldMaker:
prefix + 'zeromode', prefix + 'zeromode',
N) N)
if total_N > 0: if total_N > 0:
zm = _Distributor(dofdex,zm.target,UnstructuredDomain(total_N)) @ zm zm = _Distributor(dofdex, zm.target, UnstructuredDomain(total_N)) @ zm
return CorrelatedFieldMaker(zm, prefix, total_N) return CorrelatedFieldMaker(zm, prefix, total_N)
def add_fluctuations(self, def add_fluctuations(self,
...@@ -375,10 +375,10 @@ class CorrelatedFieldMaker: ...@@ -375,10 +375,10 @@ class CorrelatedFieldMaker:
asperity_stddev, asperity_stddev,
loglogavgslope_mean, loglogavgslope_mean,
loglogavgslope_stddev, loglogavgslope_stddev,
prefix = '', prefix='',
index = None, index=None,
dofdex = None, dofdex=None,
harmonic_partner = None): harmonic_partner=None):
if harmonic_partner is None: if harmonic_partner is None:
harmonic_partner = position_space.get_default_codomain() harmonic_partner = position_space.get_default_codomain()
else: else:
...@@ -399,7 +399,7 @@ class CorrelatedFieldMaker: ...@@ -399,7 +399,7 @@ class CorrelatedFieldMaker:
N = 0 N = 0
position_space = makeDomain(position_space) position_space = makeDomain(position_space)
prefix = str(prefix) prefix = str(prefix)
#assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace) # assert isinstance(position_space[space], (RGSpace, HPSpace, GLSpace)
fluct = _LognormalMomentMatching(fluctuations_mean, fluct = _LognormalMomentMatching(fluctuations_mean,
fluctuations_stddev, fluctuations_stddev,
...@@ -409,12 +409,12 @@ class CorrelatedFieldMaker: ...@@ -409,12 +409,12 @@ class CorrelatedFieldMaker:
self._prefix + prefix + 'flexibility', self._prefix + prefix + 'flexibility',
N) N)
asp = _LognormalMomentMatching(asperity_mean, asperity_stddev, asp = _LognormalMomentMatching(asperity_mean, asperity_stddev,
self._prefix + prefix + 'asperity', self._prefix + prefix + 'asperity',
N) N)
avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev, avgsl = _normal(loglogavgslope_mean, loglogavgslope_stddev,
self._prefix + prefix + 'loglogavgslope', N) self._prefix + prefix + 'loglogavgslope', N)
amp = _Amplitude(PowerSpace(harmonic_partner), amp = _Amplitude(PowerSpace(harmonic_partner),
fluct, flex, asp, avgsl, self._azm, fluct, flex, asp, avgsl, self._azm,
position_space[-1].total_volume, position_space[-1].total_volume,
self._prefix + prefix + 'spectrum', dofdex) self._prefix + prefix + 'spectrum', dofdex)
...@@ -431,7 +431,8 @@ class CorrelatedFieldMaker: ...@@ -431,7 +431,8 @@ class CorrelatedFieldMaker:
n_amplitudes = len(self._a) n_amplitudes = len(self._a)
if self._total_N > 0: if self._total_N > 0:
hspace = makeDomain([UnstructuredDomain(self._total_N)] + hspace = makeDomain([UnstructuredDomain(self._total_N)] +
[dd.target[-1].harmonic_partner for dd in self._a]) [dd.target[-1].harmonic_partner
for dd in self._a])
spaces = list(1 + np.arange(n_amplitudes)) spaces = list(1 + np.arange(n_amplitudes))
else: else:
hspace = makeDomain( hspace = makeDomain(
...@@ -439,28 +440,28 @@ class CorrelatedFieldMaker: ...@@ -439,28 +440,28 @@ class CorrelatedFieldMaker:
spaces = tuple(range(n_amplitudes)) spaces = tuple(range(n_amplitudes))
spaces = list(np.arange(n_amplitudes)) spaces = list(np.arange(n_amplitudes))
expander = ContractionOperator(hspace, spaces = spaces).adjoint expander = ContractionOperator(hspace, spaces=spaces).adjoint
azm = expander @ self._azm azm = expander @ self._azm
#spaces = np.array(range(n_amplitudes)) + 1 - 1//self._total_N # spaces = np.array(range(n_amplitudes)) + 1 - 1//self._total_N
ht = HarmonicTransformOperator(hspace, ht = HarmonicTransformOperator(hspace,
self._position_spaces[0][self._spaces[0]], self._position_spaces[0][self._spaces[0]],
space=spaces[0]) space=spaces[0])
for i in range(1, n_amplitudes): for i in range(1, n_amplitudes):
ht = (HarmonicTransformOperator(ht.target, ht = (HarmonicTransformOperator(ht.target,
self._position_spaces[i][self._spaces[i]], self._position_spaces[i][self._spaces[i]],
space=spaces[i]) @ ht) space=spaces[i]) @ ht)
pd = PowerDistributor(hspace, self._a[0].target[self._spaces[0]], self._spaces[0]) pd = PowerDistributor(hspace, self._a[0].target[self._spaces[0]], self._spaces[0])
for i in range(1, n_amplitudes): for i in range(1, n_amplitudes):
pd = (pd @ PowerDistributor(pd.domain, pd = (pd @ PowerDistributor(pd.domain,
self._a[i].target[self._spaces[i]], self._a[i].target[self._spaces[i]],
space=spaces[i])) space=spaces[i]))
a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0] a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ self._a[0]
for i in range(1, n_amplitudes): for i in range(1, n_amplitudes):
co = ContractionOperator(pd.domain, co = ContractionOperator(pd.domain,
spaces[:i] + spaces[i+1:]) spaces[:i] + spaces[i+1:])
a = a*(co.adjoint @ self._a[i]) a = a*(co.adjoint @ self._a[i])
return ht(azm*(pd @ a)*ducktape(hspace, None, self._prefix + 'xi')) return ht(azm*(pd @ a)*ducktape(hspace, None, self._prefix + 'xi'))
...@@ -486,7 +487,6 @@ class CorrelatedFieldMaker: ...@@ -486,7 +487,6 @@ class CorrelatedFieldMaker:
lst = [('Offset amplitude', self.amplitude_total_offset), lst = [('Offset amplitude', self.amplitude_total_offset),
('Total fluctuation amplitude', self.total_fluctuation)] ('Total fluctuation amplitude', self.total_fluctuation)]
namps = len(self._a) namps = len(self._a)
if namps > 1: if namps > 1:
for ii in range(namps): for ii in range(namps):
...@@ -507,7 +507,7 @@ class CorrelatedFieldMaker: ...@@ -507,7 +507,7 @@ class CorrelatedFieldMaker:
scm = 1. scm = 1.
for a in self._a: for a in self._a:
op = a.fluctuation_amplitude*self._azm.one_over() op = a.fluctuation_amplitude*self._azm.one_over()
res= np.array([op(from_random('normal',op.domain)).to_global_data() res = np.array([op(from_random('normal', op.domain)).to_global_data()
for _ in range(nsamples)]) for _ in range(nsamples)])
scm *= res**2 + 1. scm *= res**2 + 1.
return fluctuations_slice_mean/np.mean(np.sqrt(scm)) return fluctuations_slice_mean/np.mean(np.sqrt(scm))
......
...@@ -158,8 +158,8 @@ class Operator(metaclass=NiftyMeta): ...@@ -158,8 +158,8 @@ class Operator(metaclass=NiftyMeta):
return None, self return None, self
for f in ["sqrt", "exp", "log", "tanh", "sigmoid", 'sin', 'cos', 'tan', for f in ["sqrt", "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "tanh",
'sinh', 'cosh', 'absolute', 'sinc', 'one_over', 'log10']: "sinc", "sigmoid", "absolute", "one_over", "log10", "log1p", "expm1"]:
def func(f): def func(f):
def func2(self): def func2(self):
fa = _FunctionApplier(self.target, f) fa = _FunctionApplier(self.target, f)
......
...@@ -182,7 +182,7 @@ def _makeplot(name, block=True, dpi=None): ...@@ -182,7 +182,7 @@ def _makeplot(name, block=True, dpi=None):
return return
extension = os.path.splitext(name)[1] extension = os.path.splitext(name)[1]
if extension in (".pdf", ".png", ".svg"): if extension in (".pdf", ".png", ".svg"):
args= {} args = {}
if dpi is not None: if dpi is not None:
args['dpi'] = float(dpi) args['dpi'] = float(dpi)
plt.savefig(name, **args) plt.savefig(name, **args)
......
...@@ -38,7 +38,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator', ...@@ -38,7 +38,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator',
'full', 'from_global_data', 'from_local_data', 'full', 'from_global_data', 'from_local_data',
'makeDomain', 'sqrt', 'exp', 'log', 'tanh', 'sigmoid', 'makeDomain', 'sqrt', 'exp', 'log', 'tanh', 'sigmoid',
'sin', 'cos', 'tan', 'sinh', 'cosh', 'log10', 'sin', 'cos', 'tan', 'sinh', 'cosh', 'log10',
'absolute', 'one_over', 'clip', 'sinc', 'absolute', 'one_over', 'clip', 'sinc', "log1p", "expm1",
'conjugate', 'get_signal_variance', 'makeOp', 'domain_union', 'conjugate', 'get_signal_variance', 'makeOp', 'domain_union',
'get_default_codomain', 'single_plot', 'exec_time', 'get_default_codomain', 'single_plot', 'exec_time',
'calculate_position'] 'calculate_position']
......
...@@ -335,7 +335,7 @@ def test_emptydomain(): ...@@ -335,7 +335,7 @@ def test_emptydomain():
@pmp('dom', [ift.RGSpace((8,), harmonic=True), ()]) @pmp('dom', [ift.RGSpace((8,), harmonic=True), ()])
@pmp('func', [ @pmp('func', [
"exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc", "absolute", "exp", "log", "sin", "cos", "tan", "sinh", "cosh", "sinc", "absolute",
"sign", "log10" "sign", "log10", "log1p", "expm1"
]) ])
def test_funcs(num, dom, func): def test_funcs(num, dom, func):
num = 5 num = 5
......
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