diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py
index acc270d5b1002b83d6d929bbb04b1ccce8a1f713..edf136e3bb990faa50e36c6ee1594cb98f21ec36 100644
--- a/demos/critical_filtering.py
+++ b/demos/critical_filtering.py
@@ -97,10 +97,8 @@ if __name__ == "__main__":
# The information source
j = R.adjoint_times(N.inverse_times(d))
- realized_power = log(sh.power_analyze(logarithmic=p_space.config["logarithmic"],
- nbin=p_space.config["nbin"]))
- data_power = log(fft(d).power_analyze(logarithmic=p_space.config["logarithmic"],
- nbin=p_space.config["nbin"]))
+ realized_power = log(sh.power_analyze(binbounds=p_space.binbounds))
+ data_power = log(fft(d).power_analyze(binbounds=p_space.binbounds))
d_data = d.val.get_full_data().real
if rank == 0:
pl.plot([go.Heatmap(z=d_data)], filename='data.html')
diff --git a/nifty/field.py b/nifty/field.py
index 60710495aecd7e0efde606612286de8d0b91619f..ee12368b7421bef1747cca1a6a08401e9d0253da 100644
--- a/nifty/field.py
+++ b/nifty/field.py
@@ -274,7 +274,7 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
- def power_analyze(self, spaces=None, logarithmic=False, nbin=None,
+ def power_analyze(self, spaces=None, logarithmic=None, nbin=None,
binbounds=None, keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`.
@@ -291,14 +291,15 @@ class Field(Loggable, Versionable, object):
(default : None).
logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning.
- {default : False}
+ {default : None}
nbin : int *optional*
The number of bins the resulting PowerSpace shall have
(default : None).
if nbin==None : maximum number of bins is used
binbounds : array-like *optional*
Inner bounds of the bins (default : None).
- if binbounds==None : bins are inferred. Overwrites nbins and log
+ Overrides nbin and logarithmic.
+ if binbounds==None : bins are inferred.
keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum
of the input Field.
@@ -401,14 +402,9 @@ class Field(Loggable, Versionable, object):
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
- # extract pindex and rho from power_domain
- pindex = power_domain.pindex
- rho = power_domain.rho
-
power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val,
- pindex=pindex,
- rho=rho,
+ pdomain=power_domain,
axes=work_field.domain_axes[space_index])
# create the result field and put power_spectrum into it
@@ -425,8 +421,11 @@ class Field(Loggable, Versionable, object):
return result_field
@classmethod
- def _calculate_power_spectrum(cls, field_val, pindex, rho, axes=None):
+ def _calculate_power_spectrum(cls, field_val, pdomain, axes=None):
+ pindex = pdomain.pindex
+ # MR FIXME: how about iterating over slices, instead of replicating
+ # pindex? Would save memory and probably isn't slower.
if axes is not None:
pindex = cls._shape_up_pindex(
pindex=pindex,
@@ -435,6 +434,7 @@ class Field(Loggable, Versionable, object):
axes=axes)
power_spectrum = pindex.bincount(weights=field_val,
axis=axes)
+ rho = pdomain.rho
if axes is not None:
new_rho_shape = [1, ] * len(power_spectrum.shape)
new_rho_shape[axes[0]] = len(rho)
@@ -759,7 +759,7 @@ class Field(Loggable, Versionable, object):
Returns
-------
out : tuple
- The output object. The tuple contains the dimansions of the spaces
+ The output object. The tuple contains the dimensions of the spaces
in domain.
See Also
diff --git a/nifty/library/critical_filter/critical_power_energy.py b/nifty/library/critical_filter/critical_power_energy.py
index 7848aee5901b71af033225fb1404a329bfe73411..e94e2d35c97be477d11f96792300897da6e1dd55 100644
--- a/nifty/library/critical_filter/critical_power_energy.py
+++ b/nifty/library/critical_filter/critical_power_energy.py
@@ -103,14 +103,12 @@ class CriticalPowerEnergy(Energy):
posterior_sample = generate_posterior_sample(
self.m, self.D)
projected_sample = posterior_sample.power_analyze(
- logarithmic=self.position.domain[0].config["logarithmic"],
- nbin=self.position.domain[0].config["nbin"])
+ binbounds=self.position.domain[0].binbounds)
w += (projected_sample) * self.rho
w /= float(self.samples)
else:
w = self.m.power_analyze(
- logarithmic=self.position.domain[0].config["logarithmic"],
- nbin=self.position.domain[0].config["nbin"])
+ binbounds=self.position.domain[0].binbounds)
w *= self.rho
self._w = w
return self._w
diff --git a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
index 4517ff24d12d5aa00e21a21fe19d30dcfd680950..10b52362331aabfd87b6fef49e3ad07774c7a7b8 100644
--- a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
+++ b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
@@ -39,7 +39,7 @@ class InvertibleOperatorMixin(object):
(default: ConjugateGradient)
preconditioner : LinearOperator
- Preconditioner that is used by ConjugateGraduent if no minimizer was
+ Preconditioner that is used by ConjugateGradient if no minimizer was
given.
Attributes
diff --git a/nifty/operators/response_operator/response_operator.py b/nifty/operators/response_operator/response_operator.py
index b23ab13dbf2a43771327d7d57d45ebf019a1bea8..2c86031721b440e7ceb2ef941c12686ef1557a80 100644
--- a/nifty/operators/response_operator/response_operator.py
+++ b/nifty/operators/response_operator/response_operator.py
@@ -34,8 +34,7 @@ class ResponseOperator(LinearOperator):
domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
The domain on which the Operator's input Field lives.
target : tuple of DomainObjects, i.e. Spaces and FieldTypes
- The domain in which the outcome of the operator lives. As the Operator
- is endomorphic this is the same as its domain.
+ The domain in which the outcome of the operator lives.
unitary : boolean
Indicates whether the Operator is unitary or not.
diff --git a/nifty/spaces/lm_space/lm_space.py b/nifty/spaces/lm_space/lm_space.py
index 7fa584220093e264f528a0bfca94a0e5405a61bb..66c9f08608dff9fff48cfbf4a2e6ea6f7c132cb0 100644
--- a/nifty/spaces/lm_space/lm_space.py
+++ b/nifty/spaces/lm_space/lm_space.py
@@ -127,6 +127,22 @@ class LMSpace(Space):
return x.copy()
def get_distance_array(self, distribution_strategy):
+ if distribution_strategy == 'not': # short cut
+ lmax = self.lmax
+ ldist = np.empty((self.dim,), dtype=np.float64)
+ ldist[0:lmax+1] = np.arange(lmax+1, dtype=np.float64)
+ tmp = np.empty((2*lmax+2), dtype=np.float64)
+ tmp[0::2] = np.arange(lmax+1)
+ tmp[1::2] = np.arange(lmax+1)
+ idx = lmax+1
+ for l in range(1, lmax+1):
+ ldist[idx:idx+2*(lmax+1-l)] = tmp[2*l:]
+ idx += 2*(lmax+1-l)
+ dists = arange(start=0, stop=self.shape[0],
+ distribution_strategy=distribution_strategy)
+ dists.set_local_data(ldist)
+ return dists
+
dists = arange(start=0, stop=self.shape[0],
distribution_strategy=distribution_strategy)
@@ -136,6 +152,12 @@ class LMSpace(Space):
return dists
+ def get_unique_distances(self):
+ return np.arange(self.lmax+1, dtype=np.float64)
+
+ def get_natural_binbounds(self):
+ return np.arange(self.lmax, dtype=np.float64) + 0.5
+
@staticmethod
def _distance_array_helper(index_array, lmax):
u = 2*lmax + 1
diff --git a/nifty/spaces/power_space/__init__.py b/nifty/spaces/power_space/__init__.py
index 8f160171c3147480f2848074c3a8b61e3c134a6e..36c265fcfa9ad642f552148a589e7eb5bd440568 100644
--- a/nifty/spaces/power_space/__init__.py
+++ b/nifty/spaces/power_space/__init__.py
@@ -17,4 +17,3 @@
# and financially supported by the Studienstiftung des deutschen Volkes.
from .power_space import PowerSpace
-from .power_index_factory import PowerIndexFactory
diff --git a/nifty/spaces/power_space/power_index_factory.py b/nifty/spaces/power_space/power_index_factory.py
deleted file mode 100644
index 545e47744a908b35440edc4a9fe1446d3008180d..0000000000000000000000000000000000000000
--- a/nifty/spaces/power_space/power_index_factory.py
+++ /dev/null
@@ -1,44 +0,0 @@
-# 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-2017 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
-# and financially supported by the Studienstiftung des deutschen Volkes.
-
-from builtins import object
-from .power_indices import PowerIndices
-
-
-class _PowerIndexFactory(object):
- def __init__(self):
- self.power_indices_storage = {}
-
- def get_power_index(self, domain, distribution_strategy,
- logarithmic=False, nbin=None, binbounds=None):
- key = (domain, distribution_strategy)
-
- if key not in self.power_indices_storage:
- self.power_indices_storage[key] = \
- PowerIndices(domain, distribution_strategy,
- logarithmic=logarithmic,
- nbin=nbin,
- binbounds=binbounds)
- power_indices = self.power_indices_storage[key]
- power_index = power_indices.get_index_dict(logarithmic=logarithmic,
- nbin=nbin,
- binbounds=binbounds)
- return power_index
-
-
-PowerIndexFactory = _PowerIndexFactory()
diff --git a/nifty/spaces/power_space/power_indices.py b/nifty/spaces/power_space/power_indices.py
deleted file mode 100644
index bb210be857d71b018ea320383ef4794925feff47..0000000000000000000000000000000000000000
--- a/nifty/spaces/power_space/power_indices.py
+++ /dev/null
@@ -1,422 +0,0 @@
-# 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-2017 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
-# and financially supported by the Studienstiftung des deutschen Volkes.
-
-from __future__ import division
-from builtins import range
-from builtins import object
-import numpy as np
-from d2o import distributed_data_object,\
- STRATEGIES as DISTRIBUTION_STRATEGIES
-
-from d2o.config import dependency_injector as d2o_di
-from d2o.config import configuration as d2o_config
-
-
-class PowerIndices(object):
- """Computes helpful quantities to deal with power spectra.
-
- Given the shape and the density of a underlying rectangular grid this
- class provides the user
- with the pindex, kindex, rho and pundex. The indices are binned
- according to the supplied parameter scheme. If wanted, computed
- results are stored for future reuse.
-
- Parameters
- ----------
- domain : NIFTy harmonic space
- The space for which the power indices get computed
- distribution_strategy : str
- The distribution_strategy that will be used for the k_array and pindex
- distributed_data_object.
- logarithmic : bool *optional*
- Flag specifying if the binning of the default indices is
- performed on logarithmic scale.
- nbin : integer *optional*
- Number of used bins for the binning of the default indices.
- binbounds : {list, array}
- Array-like inner boundaries of the used bins of the default
- indices.
- """
- def __init__(self, domain, distribution_strategy,
- logarithmic=False, nbin=None, binbounds=None):
- self.domain = domain
- self.distribution_strategy = distribution_strategy
-
- # Compute the global k_array
- self.k_array = self.domain.get_distance_array(distribution_strategy)
- # Initialize the dictionary which stores all individual index-dicts
- self.global_dict = {}
- # Set self.default_parameters
- self.set_default(config_dict={'logarithmic': logarithmic,
- 'nbin': nbin,
- 'binbounds': binbounds})
-
- # Redirect the direct calls approaching a power_index instance to the
- # default_indices dict
- @property
- def default_indices(self):
- return self.get_index_dict(**self.default_parameters)
-
- def __getitem__(self, x):
- return self.default_indices.get(x)
-
- def __contains__(self, x):
- return self.default_indices.__contains__(x)
-
- def __iter__(self):
- return self.default_indices.__iter__()
-
- def __getattr__(self, x):
- return self.default_indices.__getattribute__(x)
-
- def set_default(self, **kwargs):
- """
- Sets the index-set which is specified by the parameters as the
- default for the power_index instance.
-
- Parameters
- ----------
- logarithmic : bool
- Flag specifying if the binning is performed on logarithmic
- scale.
- nbin : integer
- Number of used bins.
- binbounds : {list, array}
- Array-like inner boundaries of the used bins.
-
- Returns
- -------
- None
- """
- parsed_kwargs = self._cast_config(**kwargs)
- self.default_parameters = parsed_kwargs
-
- def _cast_config(self, **kwargs):
- """
- internal helper function which casts the various combinations of
- possible parameters into a properly defaulted dictionary
- """
- temp_config_dict = kwargs.get('config_dict', None)
- if temp_config_dict is not None:
- return self._cast_config_helper(**temp_config_dict)
- else:
- defaults = self.default_parameters
- temp_logarithmic = kwargs.get("logarithmic",
- defaults['logarithmic'])
- temp_nbin = kwargs.get("nbin", defaults['nbin'])
- temp_binbounds = kwargs.get("binbounds", defaults['binbounds'])
-
- return self._cast_config_helper(logarithmic=temp_logarithmic,
- nbin=temp_nbin,
- binbounds=temp_binbounds)
-
- def _cast_config_helper(self, logarithmic, nbin, binbounds):
- """
- internal helper function which sets the defaults for the
- _cast_config function
- """
-
- try:
- temp_logarithmic = bool(logarithmic)
- except(TypeError):
- temp_logarithmic = False
-
- try:
- temp_nbin = int(nbin)
- except(TypeError):
- temp_nbin = None
-
- try:
- temp_binbounds = tuple(np.array(binbounds))
- except(TypeError):
- temp_binbounds = None
-
- temp_dict = {"logarithmic": temp_logarithmic,
- "nbin": temp_nbin,
- "binbounds": temp_binbounds}
- return temp_dict
-
- def get_index_dict(self, **kwargs):
- """
- Returns a dictionary containing the pindex, kindex, rho and pundex
- binned according to the supplied parameter scheme and a
- configuration dict containing this scheme.
-
- Parameters
- ----------
- store : bool
- Flag specifying if the calculated index dictionary should be
- stored in the global_dict for future use.
- logarithmic : bool
- Flag specifying if the binning is performed on logarithmic
- scale.
- nbin : integer
- Number of used bins.
- binbounds : {list, array}
- Array-like inner boundaries of the used bins.
-
- Returns
- -------
- index_dict : dict
- Contains the keys: 'config', 'pindex', 'kindex', 'rho' and
- 'pundex'
- """
- # Cast the input arguments
- temp_config_dict = self._cast_config(**kwargs)
- # Compute a hashable identifier from the config which will be used
- # as dict key
- temp_key = self._freeze_config(temp_config_dict)
- # Check if the result should be stored for future use.
- storeQ = kwargs.get("store", True)
- # Try to find the requested index dict in the global_dict
- try:
- return self.global_dict[temp_key]
- except(KeyError):
- # If it is not found, calculate it.
- temp_index_dict = self._compute_index_dict(temp_config_dict)
- # Store it, if required
- if storeQ:
- self.global_dict[temp_key] = temp_index_dict
- # Important: If the result is stored, return a reference to
- # the dictionary entry, not anly a plain copy. Otherwise,
- # set_default breaks!
- return self.global_dict[temp_key]
- else:
- # Return the plain result.
- return temp_index_dict
-
- def _freeze_config(self, config_dict):
- """
- a helper function which forms a hashable identifying object from
- a config dictionary which can be used as key of a dict
- """
- return frozenset(list(config_dict.items()))
-
- def _compute_index_dict(self, config_dict):
- """
- Internal helper function which takes a config_dict, asks for the
- pindex/kindex/rho/pundex set, and bins them according to the config
- """
- # if no binning is requested, compute the indices, build the dict,
- # and return it straight.
- if not config_dict["logarithmic"] and config_dict["nbin"] is None and \
- config_dict["binbounds"] is None:
- (temp_pindex, temp_kindex, temp_rho, temp_pundex) =\
- self._compute_indices(self.k_array)
- temp_k_array = self.k_array
-
- # if binning is required, make a recursive call to get the unbinned
- # indices, bin them, compute the pundex and then return everything.
- else:
- # Get the unbinned indices
- temp_unbinned_indices = self.get_index_dict(nbin=None,
- binbounds=None,
- logarithmic=False,
- store=False)
- # Bin them
- (temp_pindex, temp_kindex, temp_rho, temp_pundex) = \
- self._bin_power_indices(
- temp_unbinned_indices, **config_dict)
- # Make a binned version of k_array
- temp_k_array = self._compute_k_array_from_pindex_kindex(
- temp_pindex, temp_kindex)
-
- temp_index_dict = {"config": config_dict,
- "pindex": temp_pindex,
- "kindex": temp_kindex,
- "rho": temp_rho,
- "pundex": temp_pundex,
- "k_array": temp_k_array}
- return temp_index_dict
-
- def _compute_k_array_from_pindex_kindex(self, pindex, kindex):
- tempindex = pindex.copy(dtype=kindex.dtype)
- result = tempindex.apply_scalar_function(
- lambda x: kindex[x.astype(np.dtype('int'))])
- return result
-
- def _compute_indices(self, k_array):
- """
- Internal helper function which computes pindex, kindex, rho and pundex
- from a given k_array
- """
- ##########
- # kindex #
- ##########
- global_kindex = k_array.unique()
-
- ##########
- # pindex #
- ##########
- # compute the local pindex slice on basis of the local k_array data
- local_pindex = np.searchsorted(global_kindex, k_array.get_local_data())
- # prepare the distributed_data_object
- global_pindex = distributed_data_object(
- global_shape=k_array.shape,
- dtype=local_pindex.dtype,
- distribution_strategy=self.distribution_strategy)
- # store the local pindex data in the global_pindex d2o
- global_pindex.set_local_data(local_pindex)
-
- #######
- # rho #
- #######
- global_rho = global_pindex.bincount().get_full_data()
-
- ##########
- # pundex #
- ##########
- global_pundex = self._compute_pundex(global_pindex,
- global_kindex)
-
- return global_pindex, global_kindex, global_rho, global_pundex
-
- def _compute_pundex(self, global_pindex, global_kindex):
- """
- Internal helper function which computes the pundex array from a
- pindex and a kindex array. This function is separated from the
- _compute_indices function as it is needed in _bin_power_indices,
- too.
- """
- if self.distribution_strategy in DISTRIBUTION_STRATEGIES['slicing']:
- ##########
- # pundex #
- ##########
- # Prepare the local data
- local_pindex = global_pindex.get_local_data()
- # Compute the local pundices for the local pindices
- (temp_uniqued_pindex, local_temp_pundex) = np.unique(
- local_pindex,
- return_index=True)
- # cure a bug in numpy
- # https://github.com/numpy/numpy/issues/9405
- local_temp_pundex = np.asarray(local_temp_pundex, dtype=np.int)
- # Shift the local pundices by the nodes' local_dim_offset
- local_temp_pundex += global_pindex.distributor.local_dim_offset
-
- # Prepare the pundex arrays used for the Allreduce operation
- # pundex has the same length as the kindex array
- local_pundex = np.zeros(shape=global_kindex.shape, dtype=np.int)
- # Set the default value higher than the maximal possible pundex
- # value so that MPI.MIN can sort out the default
- local_pundex += np.prod(global_pindex.shape) + 1
- # Set the default value higher than the length
- global_pundex = np.empty_like(local_pundex)
- # Store the individual pundices in the local_pundex array
- local_pundex[temp_uniqued_pindex] = local_temp_pundex
- # Extract the MPI module from the global_pindex d2o
- MPI = d2o_di[d2o_config['mpi_module']]
- # Use Allreduce to find the first occurences/smallest pundices
- global_pindex.comm.Allreduce(local_pundex,
- global_pundex,
- op=MPI.MIN)
- return global_pundex
-
- elif self.distribution_strategy in DISTRIBUTION_STRATEGIES['not']:
- ##########
- # pundex #
- ##########
- pundex = np.unique(global_pindex.get_local_data(),
- return_index=True)[1]
- return pundex
- else:
- raise NotImplementedError(
- "_compute_pundex_d2o not implemented for given "
- "distribution_strategy")
-
- def _bin_power_indices(self, index_dict, **kwargs):
- """
- Returns the binned power indices associated with the Fourier grid.
-
- Parameters
- ----------
- pindex : distributed_data_object
- Index of the Fourier grid points in a distributed_data_object.
- kindex : ndarray
- Array of all k-vector lengths.
- rho : ndarray
- Degeneracy factor of the individual k-vectors.
- logarithmic : bool
- Flag specifying if the binning is performed on logarithmic
- scale.
- nbin : integer
- Number of used bins.
- binbounds : {list, array}
- Array-like inner boundaries of the used bins.
-
- Returns
- -------
- pindex : distributed_data_object
- kindex, rho, pundex : ndarrays
- The (re)binned power indices.
-
- """
- # Cast the given config
- temp_config_dict = self._cast_config(**kwargs)
- logarithmic = temp_config_dict['logarithmic']
- nbin = temp_config_dict['nbin']
- binbounds = temp_config_dict['binbounds']
-
- # Extract the necessary indices from the supplied index dict
- pindex = index_dict["pindex"]
- kindex = index_dict["kindex"]
- rho = index_dict["rho"]
-
- # boundaries
- if(binbounds is not None):
- binbounds = np.sort(binbounds)
- # equal binning
- else:
- if(logarithmic is None):
- logarithmic = False
- if(logarithmic):
- k = np.r_[0, np.log(kindex[1:])]
- else:
- k = kindex
- dk = np.max(k[2:] - k[1:-1]) # minimal dk
- if(nbin is None):
- nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
- dk - 0.5) # maximal nbin
- else:
- nbin = min(int(nbin), int(
- (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
- dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
- binbounds = np.r_[0.5 * (3 * k[1] - k[2]),
- 0.5 * (k[1] + k[2]) + dk * np.arange(nbin - 2)]
- if(logarithmic):
- binbounds = np.exp(binbounds)
- # reordering
- reorder = np.searchsorted(binbounds, kindex)
- rho_ = np.zeros(len(binbounds) + 1, dtype=rho.dtype)
- kindex_ = np.empty(len(binbounds) + 1, dtype=kindex.dtype)
- for ii in range(len(reorder)):
- if(rho_[reorder[ii]] == 0):
- kindex_[reorder[ii]] = kindex[ii]
- rho_[reorder[ii]] += rho[ii]
- else:
- kindex_[reorder[ii]] = ((kindex_[reorder[ii]] *
- rho_[reorder[ii]] +
- kindex[ii] * rho[ii]) /
- (rho_[reorder[ii]] + rho[ii]))
- rho_[reorder[ii]] += rho[ii]
-
- pindex_ = pindex.copy_empty()
- pindex_.set_local_data(reorder[pindex.get_local_data()])
- pundex_ = self._compute_pundex(pindex_, kindex_)
-
- return pindex_, kindex_, rho_, pundex_
diff --git a/nifty/spaces/power_space/power_space.py b/nifty/spaces/power_space/power_space.py
index d96ff42cd09dc5dd4e4d622c4c8dafc607fef136..84ebcf186ff183887fbc31e4c85472c3963e34f4 100644
--- a/nifty/spaces/power_space/power_space.py
+++ b/nifty/spaces/power_space/power_space.py
@@ -19,9 +19,7 @@
# from builtins import str
import numpy as np
-import d2o
-
-from .power_index_factory import PowerIndexFactory
+from d2o import distributed_data_object
from ..space import Space
from functools import reduce
@@ -39,18 +37,24 @@ class PowerSpace(Space):
derived from this PowerSpace, e.g. the pindex.
(default : 'not')
logarithmic : bool *optional*
- True if logarithmic binning should be used (default : False).
+ True if logarithmic binning should be used (default : None).
nbin : {int, None} *optional*
The number of bins that should be used for power spectrum binning
(default : None).
if nbin == None, then nbin is set to the length of kindex.
binbounds : {list, array-like} *optional*
- Array-like inner boundaries of the used bins of the default
- indices.
+ Boundaries between the power spectrum bins.
+ (If binbounds has n entries, there will be n+1 bins, the first bin
+ starting at -inf, the last bin ending at +inf.)
(default : None)
- if binbounds == None :
+ if binbounds == None:
Calculates the bounds from the kindex while applying the
logarithmic and nbin keywords.
+ Note: if "bindounds" is not None, both "logarithmic" and "nbin" must be
+ None!
+ Note: if "binbounds", "logarithmic", and "nbin" are all None, then
+ "natural" binning is performed, i.e. there will be one bin for every
+ distinct k-vector length.
Attributes
----------
@@ -59,9 +63,6 @@ class PowerSpace(Space):
mapped to which power bin
kindex : numpy.ndarray
Sorted array of all k-modes.
- pundex : numpy.ndarray
- Flat index of the first occurence of a k-vector with length==kindex[n]
- in the k_array.
rho : numpy.ndarray
The amount of k-modes that get mapped to one power bin is given by
rho.
@@ -73,9 +74,8 @@ class PowerSpace(Space):
The total volume of the space.
shape : tuple of np.ints
The shape of the space's data array.
- config : {logarithmic, nbin, binbounds}
- Dictionary storing the values for `logarithmic`, `nbin`, and
- `binbounds` that were used during initialization.
+ binbounds : tuple or None
+ Boundaries between the power spectrum bins
Notes
-----
@@ -84,42 +84,107 @@ class PowerSpace(Space):
"""
+ _powerIndexCache = {}
+
# ---Overwritten properties and methods---
- def __init__(self, harmonic_partner,
- distribution_strategy='not',
- logarithmic=False, nbin=None, binbounds=None):
+ def __init__(self, harmonic_partner, distribution_strategy='not',
+ logarithmic=None, nbin=None, binbounds=None):
super(PowerSpace, self).__init__()
- self._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
- '_k_array']
+ self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
- if not isinstance(harmonic_partner, Space):
- raise ValueError(
- "harmonic_partner must be a Space.")
- if not harmonic_partner.harmonic:
- raise ValueError(
- "harmonic_partner must be a harmonic space.")
+ if not (isinstance(harmonic_partner, Space) and
+ harmonic_partner.harmonic):
+ raise ValueError("harmonic_partner must be a harmonic space.")
self._harmonic_partner = harmonic_partner
- power_index = PowerIndexFactory.get_power_index(
- domain=self.harmonic_partner,
- distribution_strategy=distribution_strategy,
- logarithmic=logarithmic,
- nbin=nbin,
- binbounds=binbounds)
-
- self._config = power_index['config']
+ # sanity check
+ if binbounds is not None and not(nbin is None and logarithmic is None):
+ raise ValueError(
+ "if binbounds is defined, nbin and logarithmic must be None")
+
+ if binbounds is not None:
+ binbounds = tuple(binbounds)
+
+ key = (harmonic_partner, distribution_strategy, logarithmic, nbin,
+ binbounds)
+ if self._powerIndexCache.get(key) is None:
+ distance_array = \
+ self.harmonic_partner.get_distance_array(distribution_strategy)
+ temp_binbounds = self._compute_binbounds(
+ harmonic_partner=self.harmonic_partner,
+ distribution_strategy=distribution_strategy,
+ logarithmic=logarithmic,
+ nbin=nbin,
+ binbounds=binbounds)
+ temp_pindex = self._compute_pindex(
+ harmonic_partner=self.harmonic_partner,
+ distance_array=distance_array,
+ binbounds=temp_binbounds,
+ distribution_strategy=distribution_strategy)
+ temp_rho = temp_pindex.bincount().get_full_data()
+ temp_kindex = \
+ (temp_pindex.bincount(weights=distance_array).get_full_data() /
+ temp_rho)
+ self._powerIndexCache[key] = (temp_binbounds,
+ temp_pindex,
+ temp_kindex,
+ temp_rho)
+
+ (self._binbounds, self._pindex, self._kindex, self._rho) = \
+ self._powerIndexCache[key]
+
+ @staticmethod
+ def _compute_binbounds(harmonic_partner, distribution_strategy,
+ logarithmic, nbin, binbounds):
+
+ if logarithmic is None and nbin is None and binbounds is None:
+ result = None
+ else:
+ if binbounds is not None:
+ bb = np.sort(np.array(binbounds))
+ else:
+ if logarithmic is not None:
+ logarithmic = bool(logarithmic)
+ if nbin is not None:
+ nbin = int(nbin)
+
+ # equidistant binning (linear or log)
+ # MR FIXME: this needs to improve
+ kindex = harmonic_partner.get_unique_distances()
+ if (logarithmic):
+ k = np.r_[0, np.log(kindex[1:])]
+ else:
+ k = kindex
+ dk = np.max(k[2:] - k[1:-1]) # minimum dk to avoid empty bins
+ if(nbin is None):
+ nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
+ dk - 0.5) # maximal nbin
+ else:
+ nbin = min(int(nbin), int(
+ (k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
+ dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
+ bb = np.r_[0.5 * (3 * k[1] - k[2]),
+ 0.5 * (k[1] + k[2]) + dk * np.arange(nbin-2)]
+ if(logarithmic):
+ bb = np.exp(bb)
+ result = tuple(bb)
+ return result
- self._pindex = power_index['pindex']
- self._kindex = power_index['kindex']
- self._rho = power_index['rho']
- self._pundex = power_index['pundex']
- self._k_array = power_index['k_array']
+ @staticmethod
+ def _compute_pindex(harmonic_partner, distance_array, binbounds,
+ distribution_strategy):
- if self.config['nbin'] is not None:
- if self.config['nbin'] > len(self.kindex):
- self.logger.warn("nbin was set to a value being larger than "
- "the length of kindex!")
+ # Compute pindex, kindex and rho according to bb
+ pindex = distributed_data_object(
+ global_shape=distance_array.shape,
+ dtype=np.int,
+ distribution_strategy=distribution_strategy)
+ if binbounds is None:
+ binbounds = harmonic_partner.get_natural_binbounds()
+ pindex.set_local_data(
+ np.searchsorted(binbounds, distance_array.get_local_data()))
+ return pindex
def pre_cast(self, x, axes):
""" Casts power spectrum functions to discretized power spectra.
@@ -144,19 +209,15 @@ class PowerSpace(Space):
"""
- if callable(x):
- return x(self.kindex)
- else:
- return x
+ return x(self.kindex) if callable(x) else x
# ---Mandatory properties and methods---
def __repr__(self):
return ("PowerSpace(harmonic_partner=%r, distribution_strategy=%r, "
- "logarithmic=%r, nbin=%r, binbounds=%r)"
+ "binbounds=%r)"
% (self.harmonic_partner, self.pindex.distribution_strategy,
- self.config['logarithmic'], self.config['nbin'],
- self.config['binbounds']))
+ self._binbounds))
@property
def harmonic(self):
@@ -179,11 +240,9 @@ class PowerSpace(Space):
distribution_strategy = self.pindex.distribution_strategy
return self.__class__(harmonic_partner=self.harmonic_partner,
distribution_strategy=distribution_strategy,
- logarithmic=self.config["logarithmic"],
- nbin=self.config["nbin"],
- binbounds=self.config["binbounds"])
+ binbounds=self._binbounds)
- def weight(self, x, power=1, axes=None, inplace=False):
+ def weight(self, x, power, axes, inplace=False):
reshaper = [1, ] * len(x.shape)
# we know len(axes) is always 1
reshaper[axes[0]] = self.shape[0]
@@ -201,10 +260,9 @@ class PowerSpace(Space):
return result_x
def get_distance_array(self, distribution_strategy):
- result = d2o.distributed_data_object(
+ return distributed_data_object(
self.kindex, dtype=np.float64,
distribution_strategy=distribution_strategy)
- return result
def get_fft_smoothing_kernel_function(self, sigma):
raise NotImplementedError(
@@ -219,11 +277,8 @@ class PowerSpace(Space):
return self._harmonic_partner
@property
- def config(self):
- """ Returns the configuration which was used for `logarithmic`, `nbin`
- and `binbounds` during initialization.
- """
- return self._config
+ def binbounds(self):
+ return self._binbounds
@property
def pindex(self):
@@ -244,67 +299,20 @@ class PowerSpace(Space):
"""
return self._rho
- @property
- def pundex(self):
- """ An array for which the n-th entry gives the flat index of the
- first occurence of a k-vector with length==kindex[n] in the
- k_array.
- """
- return self._pundex
-
- @property
- def k_array(self):
- """ An array containing distances to the grid center (i.e. zero-mode)
- for every k-mode in the grid of the harmonic partner space.
- """
- return self._k_array
-
# ---Serialization---
def _to_hdf5(self, hdf5_group):
- hdf5_group['kindex'] = self.kindex
- hdf5_group['rho'] = self.rho
- hdf5_group['pundex'] = self.pundex
- hdf5_group['logarithmic'] = self.config["logarithmic"]
- # Store nbin as string, since it can be None
- hdf5_group.attrs['nbin'] = str(self.config["nbin"])
- hdf5_group.attrs['binbounds'] = str(self.config["binbounds"])
-
- #MR FIXME: why not "return None" as happens everywhere else?
+ hdf5_group.attrs['binbounds'] = str(self._binbounds)
+ hdf5_group.attrs['distribution_strategy'] = \
+ self._pindex.distribution_strategy
+
return {
'harmonic_partner': self.harmonic_partner,
- 'pindex': self.pindex,
- 'k_array': self.k_array
}
@classmethod
def _from_hdf5(cls, hdf5_group, repository):
- # make an empty PowerSpace object
- new_ps = EmptyPowerSpace()
- # reset class
- new_ps.__class__ = cls
- # call instructor so that classes are properly setup
- super(PowerSpace, new_ps).__init__()
- # set all values
- new_ps._harmonic_partner = repository.get('harmonic_partner',
- hdf5_group)
-
- new_ps._config = {}
- new_ps._config['logarithmic'] = hdf5_group['logarithmic'][()]
- exec("new_ps._config['nbin'] = " + hdf5_group.attrs['nbin'])
- exec("new_ps._config['binbounds'] = " + hdf5_group.attrs['binbounds'])
-
- new_ps._pindex = repository.get('pindex', hdf5_group)
- new_ps._kindex = hdf5_group['kindex'][:]
- new_ps._rho = hdf5_group['rho'][:]
- new_ps._pundex = hdf5_group['pundex'][:]
- new_ps._k_array = repository.get('k_array', hdf5_group)
- new_ps._ignore_for_hash += ['_pindex', '_kindex', '_rho', '_pundex',
- '_k_array']
-
- return new_ps
-
-
-class EmptyPowerSpace(PowerSpace):
- def __init__(self):
- pass
+ hp = repository.get('harmonic_partner', hdf5_group)
+ exec("bb = " + hdf5_group.attrs['binbounds'])
+ ds = hdf5_group.attrs['distribution_strategy']
+ return PowerSpace(hp, ds, binbounds=bb)
diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py
index d9109607645fce4094e767bbf3e03c73ebe9d1e7..2c2ae67d5cff890a719165e440487d7699ff49d4 100644
--- a/nifty/spaces/rg_space/rg_space.py
+++ b/nifty/spaces/rg_space/rg_space.py
@@ -249,6 +249,36 @@ class RGSpace(Space):
dists = np.sqrt(dists)
return dists
+ def get_unique_distances(self):
+ dimensions = len(self.shape)
+ if dimensions == 1: # extra easy
+ maxdist = self.shape[0]//2
+ return np.arange(maxdist+1, dtype=np.float64) * self.distances[0]
+ if np.all(self.distances == self.distances[0]): # shortcut
+ maxdist = np.asarray(self.shape)//2
+ tmp = np.sum(maxdist*maxdist)
+ tmp = np.zeros(tmp+1, dtype=np.bool)
+ t2 = np.arange(maxdist[0]+1, dtype=np.int64)
+ t2 *= t2
+ for i in range(1, dimensions):
+ t3 = np.arange(maxdist[i]+1, dtype=np.int64)
+ t3 *= t3
+ t2 = np.add.outer(t2, t3)
+ tmp[t2] = True
+ return np.sqrt(np.nonzero(tmp)[0])*self.distances[0]
+ else: # do it the hard way
+ tmp = self.get_distance_array('not').unique() # expensive!
+ tol = 1e-12*tmp[-1]
+ # remove all points that are closer than tol to their right
+ # neighbors.
+ # I'm appending the last value*2 to the array to treat the
+ # rightmost point correctly.
+ return tmp[np.diff(np.r_[tmp, 2*tmp[-1]]) > tol]
+
+ def get_natural_binbounds(self):
+ tmp = self.get_unique_distances()
+ return 0.5*(tmp[:-1]+tmp[1:])
+
def get_fft_smoothing_kernel_function(self, sigma):
return lambda x: np.exp(-2. * np.pi*np.pi * x*x * sigma*sigma)
diff --git a/nifty/spaces/space/space.py b/nifty/spaces/space/space.py
index 73421ad317d344c024e5dae204e9b35a1baa6cac..da60dde377e7eaee6aace767b8d449a2a76220a3 100644
--- a/nifty/spaces/space/space.py
+++ b/nifty/spaces/space/space.py
@@ -117,6 +117,20 @@ class Space(DomainObject):
raise NotImplementedError(
"There is no generic distance structure for Space base class.")
+ def get_unique_distances(self):
+ raise NotImplementedError
+
+ def get_natural_binbounds(self):
+ """ The boundaries for natural power spectrum binning.
+
+ Returns
+ -------
+ distributed_data_object
+ A numpy array containing the binbounds
+
+ """
+ raise NotImplementedError
+
def get_fft_smoothing_kernel_function(self, sigma):
""" This method returns a smoothing kernel function.
diff --git a/test/test_operators/test_fft_operator.py b/test/test_operators/test_fft_operator.py
index a43bbea45e91cfc05976e81abe7fbf8858280b75..48122b9d95657311824996886d99879180078054 100644
--- a/test/test_operators/test_fft_operator.py
+++ b/test/test_operators/test_fft_operator.py
@@ -78,10 +78,13 @@ class FFTOperatorTests(unittest.TestCase):
b = RGSpace(dim1, zerocenter=zc2, distances=1./(dim1*d), harmonic=True)
fft = FFTOperator(domain=a, target=b, domain_dtype=itp,
target_dtype=_harmonic_type(itp), module=module)
+ np.random.seed(16)
inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3,
dtype=itp)
out = fft.adjoint_times(fft.times(inp))
- assert_allclose(inp.val, out.val, rtol=tol, atol=tol)
+ assert_allclose(inp.val.get_full_data(),
+ out.val.get_full_data(),
+ rtol=tol, atol=tol)
@expand(product(["numpy", "fftw", "fftw_mpi"],
[10, 11], [9, 12], [False, True],
diff --git a/test/test_spaces/test_power_space.py b/test/test_spaces/test_power_space.py
index 237055b9ef62c3c7e53f345e44474ee9190e6746..4d75ad467e6e01000689f1195dad24769466dcc1 100644
--- a/test/test_spaces/test_power_space.py
+++ b/test/test_spaces/test_power_space.py
@@ -47,14 +47,12 @@ HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
#Try all sensible kinds of combinations of spaces, distributuion strategy and
#binning parameters
-_maybe_fftw = ["fftw"] if ('pyfftw' in gdi) else []
-
CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES,
["not", "equal", "fftw"],
[None], [None, 3, 4], [True, False])
CONSISTENCY_CONFIGS_EXPLICIT = product(HARMONIC_SPACES,
["not", "equal", "fftw"],
- [[0., 1.3]], [None], [False])
+ [[0., 1.3]], [None], [None])
CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
CONSISTENCY_CONFIGS_EXPLICIT)
@@ -63,18 +61,16 @@ CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
CONSTRUCTOR_CONFIGS = [
[1, 'not', False, None, None, {'error': ValueError}],
[RGSpace((8,)), 'not', False, None, None, {'error': ValueError}],
- [RGSpace((8,), harmonic=True), 'not', False, None, None, {
+ [RGSpace((8,), harmonic=True), 'not', None, None, None, {
'harmonic': True,
'shape': (5,),
'dim': 5,
'total_volume': 8.0,
'harmonic_partner': RGSpace((8,), harmonic=True),
- 'config': {'logarithmic': False, 'nbin': None, 'binbounds': None},
+ 'binbounds': None,
'pindex': distributed_data_object([0, 1, 2, 3, 4, 3, 2, 1]),
'kindex': np.array([0., 1., 2., 3., 4.]),
'rho': np.array([1, 2, 2, 2, 1]),
- 'pundex': np.array([0, 1, 2, 3, 4]),
- 'k_array': np.array([0., 1., 2., 3., 4., 3., 2., 1.]),
}],
[RGSpace((8,), harmonic=True), 'not', True, None, None, {
'harmonic': True,
@@ -82,13 +78,10 @@ CONSTRUCTOR_CONFIGS = [
'dim': 2,
'total_volume': 8.0,
'harmonic_partner': RGSpace((8,), harmonic=True),
- 'config': {'logarithmic': True, 'nbin': None, 'binbounds': None},
+ 'binbounds': (0.70710678118654757,),
'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]),
'kindex': np.array([0., 2.28571429]),
'rho': np.array([1, 7]),
- 'pundex': np.array([0, 1]),
- 'k_array': np.array([0., 2.28571429, 2.28571429, 2.28571429,
- 2.28571429, 2.28571429, 2.28571429, 2.28571429]),
}],
]
@@ -121,12 +114,10 @@ def get_weight_configs():
class PowerSpaceInterfaceTest(unittest.TestCase):
@expand([
['harmonic_partner', Space],
- ['config', dict],
+ ['binbounds', type(None)],
['pindex', distributed_data_object],
['kindex', np.ndarray],
['rho', np.ndarray],
- ['pundex', np.ndarray],
- ['k_array', distributed_data_object],
])
def test_property_ret_type(self, attribute, expected_type):
r = RGSpace((4, 4), harmonic=True)
@@ -136,23 +127,8 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
class PowerSpaceConsistencyCheck(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS)
- def test_pipundexInversion(self, harmonic_partner, distribution_strategy,
- binbounds, nbin, logarithmic):
- if distribution_strategy == "fftw":
- if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
- raise SkipTest
- p = PowerSpace(harmonic_partner=harmonic_partner,
- distribution_strategy=distribution_strategy,
- logarithmic=logarithmic, nbin=nbin,
- binbounds=binbounds)
- assert_equal(p.pindex.flatten().get_full_data()[p.pundex],
- np.arange(p.dim),
- err_msg='pundex is not right-inverse of pindex!')
-
- @expand(CONSISTENCY_CONFIGS)
- def test_rhopindexConsistency(self, harmonic_partner,
- distribution_strategy, binbounds, nbin,
- logarithmic):
+ def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
+ binbounds, nbin, logarithmic):
if distribution_strategy == "fftw":
if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
raise SkipTest
@@ -160,8 +136,6 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
- assert_equal(p.pindex.flatten()[p.pundex],np.arange(p.dim),
- err_msg='pundex is not right-inverse of pindex!')
assert_equal(p.pindex.flatten().bincount(), p.rho,
err_msg='rho is not equal to pindex degeneracy')
@@ -173,7 +147,6 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
if distribution_strategy == "fftw":
if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
raise SkipTest
- raise SkipTest
if 'error' in expected:
with assert_raises(expected['error']):
PowerSpace(harmonic_partner=harmonic_partner,