Commit d698a084 authored by Jait Dixit's avatar Jait Dixit

WIP: Refactor spaces

- Remove smooth from abstract Space class
- Remove smooth from RGSpace
- Remove smooth and transform from HPSpace
- Move all HPSpace specific attributes to paradict
parent 30623565
......@@ -42,8 +42,7 @@ from nifty.spaces.lm_space import LMSpace
from nifty.spaces.space import Space
from nifty.config import about,\
nifty_configuration as gc,\
from nifty.config import about, nifty_configuration as gc,\
dependency_injector as gdi
from hp_space_paradict import HPSpaceParadict
from nifty.nifty_random import random
......@@ -127,33 +126,24 @@ class HPSpace(Space):
"""
# check imports
if not gc['use_healpy']:
raise ImportError(about._errors.cstring(
"ERROR: healpy not available."))
raise ImportError("ERROR: healpy not available.")
self._cache_dict = {'check_codomain': {}}
# check parameters
self.paradict = HPSpaceParadict(nside=nside)
# let the abstract Space class initialize dtype
super(HPSpace, self).__init__(dtype=np.dtype('float64'))
self.dtype = np.dtype('float64')
self.harmonic = False
self.distances = (np.float(4*np.pi / (12*self.paradict['nside']**2)),)
@property
def para(self):
temp = np.array([self.paradict['nside']], dtype=int)
return temp
@para.setter
def para(self, x):
self.paradict['nside'] = x[0]
self.paradict = HPSpaceParadict(nside=nside,
distances=np.float(
4 * np.pi / (12 * nside ** 2))
)
# HPSpace is never harmonic
self._harmonic = False
def copy(self):
return HPSpace(nside=self.paradict['nside'])
@property
def shape(self):
return (np.int(12 * self.paradict['nside']**2),)
return (np.int(12 * self.paradict['nside'] ** 2),)
@property
def meta_volume(self):
......@@ -184,8 +174,7 @@ class HPSpace(Space):
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.float)
return self.calc_weight(mol, power=1)
pass
# TODO: Extend to binning/log
def enforce_power(self, spec, size=None, kindex=None):
......@@ -197,243 +186,6 @@ class HPSpace(Space):
size=size,
kindex=kindex)
def _check_codomain(self, codomain):
"""
Checks whether a given codomain is compatible to the space or not.
Parameters
----------
codomain : nifty.space
Space to be checked for compatibility.
Returns
-------
check : bool
Whether or not the given codomain is compatible to the space.
Notes
-----
Compatible codomains are instances of :py:class:`hp_space` and
:py:class:`lm_space`.
"""
if codomain is None:
return False
if not isinstance(codomain, Space):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if isinstance(codomain, LMSpace):
nside = self.paradict['nside']
lmax = codomain.paradict['lmax']
mmax = codomain.paradict['mmax']
# 3*nside-1==lmax
# lmax==mmax
if (3*nside-1 == lmax) and (lmax == mmax):
return True
return False
def get_codomain(self, **kwargs):
"""
Generates a compatible codomain to which transformations are
reasonable, i.e.\ an instance of the :py:class:`lm_space` class.
Returns
-------
codomain : nifty.lm_space
A compatible codomain.
"""
lmax = 3*self.paradict['nside'] - 1
mmax = lmax
return LMSpace(lmax=lmax, mmax=mmax, dtype=np.dtype('complex128'))
def get_random_values(self, **kwargs):
"""
Generates random field values according to the specifications given
by the parameters.
Returns
-------
x : numpy.ndarray
Valid field values.
Other parameters
----------------
random : string, *optional*
Specifies the probability distribution from which the random
numbers are to be drawn.
Supported distributions are:
- "pm1" (uniform distribution over {+1,-1}
- "gau" (normal distribution with zero-mean and a given
standard
deviation or variance)
- "syn" (synthesizes from a given power spectrum)
- "uni" (uniform distribution over [vmin,vmax[)
(default: None).
dev : float, *optional*
Standard deviation (default: 1).
var : float, *optional*
Variance, overriding `dev` if both are specified
(default: 1).
spec : {scalar, list, numpy.array, nifty.field, function},
*optional*
Power spectrum (default: 1).
codomain : nifty.lm_space, *optional*
A compatible codomain for power indexing (default: None).
vmin : float, *optional*
Lower limit for a uniform distribution (default: 0).
vmax : float, *optional*
Upper limit for a uniform distribution (default: 1).
"""
arg = random.parse_arguments(self, **kwargs)
# if arg is None:
# x = np.zeros(self.shape, dtype=self.dtype)
#
# elif arg['random'] == "pm1":
# x = random.pm1(dtype=self.dtype, shape=self.shape)
#
# elif arg['random'] == "gau":
# x = random.gau(dtype=self.dtype, shape=self.shape,
# mean=arg['mean'],
# std=arg['std'])
if arg['random'] == "syn":
nside = self.paradict['nside']
lmax = 3*nside-1
sample = hp.synfast(arg['spec'],
nside,
lmax=lmax, mmax=lmax,
alm=False, pol=True, pixwin=False,
fwhm=0.0, sigma=None)
else:
sample = super(HPSpace, self).get_random_values(**arg)
# elif arg['random'] == "uni":
# x = random.uni(dtype=self.dtype, shape=self.shape,
# vmin=arg['vmin'],
# vmax=arg['vmax'])
#
# else:
# raise KeyError(about._errors.cstring(
# "ERROR: unsupported random key '" + str(arg['random']) + "'."))
sample = self.cast(sample)
return sample
def calc_transform(self, x, codomain=None, niter=0, **kwargs):
"""
Computes the transform of a given array of field values.
Parameters
----------
x : numpy.ndarray
Array to be transformed.
codomain : nifty.space, *optional*
codomain space to which the transformation shall map
(default: self).
Returns
-------
Tx : numpy.ndarray
Transformed array
Other parameters
----------------
iter : int, *optional*
Number of iterations performed in the HEALPix basis
transformation.
Notes
-----
Only instances of the :py:class:`lm_space` or :py:class:`hp_space`
classes are allowed as `codomain`.
"""
x = self.cast(x)
# Check if the given codomain is suitable for the transformation
if not self.check_codomain(codomain):
raise ValueError(about._errors.cstring(
"ERROR: unsupported codomain."))
# TODO look at these kinds of checks maybe need replacement
# if self.datamodel != 'not':
# about.warnings.cprint(
# "WARNING: Field data is consolidated to all nodes for "
# "external map2alm method!")
np_x = x.get_full_data()
if isinstance(codomain, LMSpace):
# transform
np_Tx = hp.map2alm(np_x.astype(np.float64, copy=False),
lmax=codomain.paradict['lmax'],
mmax=codomain.paradict['mmax'],
iter=niter, pol=True, use_weights=False,
datapath=None)
else:
raise ValueError(about._errors.cstring(
"ERROR: unsupported transformation."))
return codomain.cast(np_Tx)
def calc_smooth(self, x, sigma=0, niter=0, **kwargs):
"""
Smoothes an array of field values by convolution with a Gaussian
kernel.
Parameters
----------
x : numpy.ndarray
Array of field values to be smoothed.
sigma : float, *optional*
Standard deviation of the Gaussian kernel, specified in units
of length in position space; for testing: a sigma of -1 will be
reset to a reasonable value (default: 0).
Returns
-------
Gx : numpy.ndarray
Smoothed array.
Other parameters
----------------
iter : int, *optional*
Number of iterations performed in the HEALPix basis
transformation.
"""
nside = self.paradict['nside']
x = self.cast(x)
# check sigma
if sigma == 0:
return self.unary_operation(x, op='copy')
elif sigma == -1:
about.infos.cprint("INFO: invalid sigma reset.")
# sqrt(2)*pi/(lmax+1)
sigma = np.sqrt(2) * np.pi / (3. * nside)
elif sigma < 0:
raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
# smooth
# if self.datamodel != 'not':
# about.warnings.cprint(
# "WARNING: Field data is consolidated to all nodes for "
# "external smoothalm method!")
np_x = x.get_full_data()
lmax = 3*nside-1
mmax = lmax
result = hp.smoothing(np_x, fwhm=0.0, sigma=sigma, pol=True,
iter=niter, lmax=lmax, mmax=mmax,
use_weights=False, datapath=None)
return self.cast(result)
def calc_power(self, x, niter=0, **kwargs):
"""
Computes the power of an array of field values.
......@@ -458,7 +210,7 @@ class HPSpace(Space):
x = self.cast(x)
nside = self.paradict['nside']
lmax = 3*nside-1
lmax = 3 * nside - 1
mmax = lmax
# if self.datamodel != 'not':
......@@ -529,49 +281,71 @@ class HPSpace(Space):
except AttributeError:
pass
if(not pl.isinteractive())and(not bool(kwargs.get("save", False))):
if (not pl.isinteractive()) and (not bool(kwargs.get("save", False))):
about.warnings.cprint("WARNING: interactive mode off.")
if(power):
if (power):
x = self.calc_power(x, **kwargs)
fig = pl.figure(num=None, figsize=(6.4, 4.8), dpi=None, facecolor="none",
edgecolor="none", frameon=False, FigureClass=pl.Figure)
fig = pl.figure(num=None, figsize=(6.4, 4.8), dpi=None,
facecolor="none",
edgecolor="none", frameon=False,
FigureClass=pl.Figure)
ax0 = fig.add_axes([0.12, 0.12, 0.82, 0.76])
xaxes = np.arange(3 * self.para[0], dtype=np.int)
if(vmin is None):
if (vmin is None):
vmin = np.min(x[:mono].tolist(
) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None, out=None)
if(vmax is None):
) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None,
out=None)
if (vmax is None):
vmax = np.max(x[:mono].tolist(
) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None, out=None)
) + (xaxes * (2 * xaxes + 1) * x)[1:].tolist(), axis=None,
out=None)
ax0.loglog(xaxes[1:], (xaxes * (2 * xaxes + 1) * x)[1:], color=[0.0,
0.5, 0.0], label="graph 0", linestyle='-', linewidth=2.0, zorder=1)
if(mono):
ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), x[0], s=20, color=[0.0, 0.5, 0.0], marker='o',
cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, zorder=1)
if(other is not None):
if(isinstance(other, tuple)):
0.5,
0.0],
label="graph 0", linestyle='-', linewidth=2.0, zorder=1)
if (mono):
ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), x[0], s=20,
color=[0.0, 0.5, 0.0], marker='o',
cmap=None, norm=None, vmin=None, vmax=None,
alpha=None, linewidths=None, verts=None, zorder=1)
if (other is not None):
if (isinstance(other, tuple)):
other = list(other)
for ii in xrange(len(other)):
if(isinstance(other[ii], Field)):
if (isinstance(other[ii], Field)):
other[ii] = other[ii].power(**kwargs)
else:
other[ii] = self.enforce_power(other[ii])
elif(isinstance(other, Field)):
elif (isinstance(other, Field)):
other = [other.power(**kwargs)]
else:
other = [self.enforce_power(other)]
imax = max(1, len(other) - 1)
for ii in xrange(len(other)):
ax0.loglog(xaxes[1:], (xaxes * (2 * xaxes + 1) * other[ii])[1:], color=[max(0.0, 1.0 - (2 * ii / imax)**2), 0.5 * ((2 * ii - imax) / imax)
** 2, max(0.0, 1.0 - (2 * (ii - imax) / imax)**2)], label="graph " + str(ii + 1), linestyle='-', linewidth=1.0, zorder=-ii)
if(mono):
ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), other[ii][0], s=20, color=[max(0.0, 1.0 - (2 * ii / imax)**2), 0.5 * ((2 * ii - imax) / imax)**2, max(
0.0, 1.0 - (2 * (ii - imax) / imax)**2)], marker='o', cmap=None, norm=None, vmin=None, vmax=None, alpha=None, linewidths=None, verts=None, zorder=-ii)
if(legend):
ax0.loglog(xaxes[1:],
(xaxes * (2 * xaxes + 1) * other[ii])[1:],
color=[max(0.0, 1.0 - (2 * ii / imax) ** 2),
0.5 * ((2 * ii - imax) / imax)
** 2, max(0.0, 1.0 - (
2 * (ii - imax) / imax) ** 2)],
label="graph " + str(ii + 1), linestyle='-',
linewidth=1.0, zorder=-ii)
if (mono):
ax0.scatter(0.5 * (xaxes[1] + xaxes[2]), other[ii][0],
s=20,
color=[max(0.0, 1.0 - (2 * ii / imax) ** 2),
0.5 * ((2 * ii - imax) / imax) ** 2,
max(
0.0, 1.0 - (
2 * (ii - imax) / imax) ** 2)],
marker='o', cmap=None, norm=None, vmin=None,
vmax=None, alpha=None, linewidths=None,
verts=None, zorder=-ii)
if (legend):
ax0.legend()
ax0.set_xlim(xaxes[1], xaxes[-1])
......@@ -581,24 +355,30 @@ class HPSpace(Space):
ax0.set_title(title)
else:
if(norm == "log"):
if(vmin is not None):
if(vmin <= 0):
if (norm == "log"):
if (vmin is not None):
if (vmin <= 0):
raise ValueError(about._errors.cstring(
"ERROR: nonpositive value(s)."))
elif(np.min(x, axis=None, out=None) <= 0):
elif (np.min(x, axis=None, out=None) <= 0):
raise ValueError(about._errors.cstring(
"ERROR: nonpositive value(s)."))
if(cmap is None):
if (cmap is None):
cmap = pl.cm.jet # default
cmap.set_under(color='k', alpha=0.0) # transparent box
hp.mollview(x, fig=None, rot=None, coord=None, unit=unit, xsize=800, title=title, nest=False, min=vmin, max=vmax, flip="astro", remove_dip=False,
remove_mono=False, gal_cut=0, format="%g", format2="%g", cbar=cbar, cmap=cmap, notext=False, norm=norm, hold=False, margins=None, sub=None)
hp.mollview(x, fig=None, rot=None, coord=None, unit=unit, xsize=800,
title=title, nest=False, min=vmin, max=vmax,
flip="astro", remove_dip=False,
remove_mono=False, gal_cut=0, format="%g", format2="%g",
cbar=cbar, cmap=cmap, notext=False, norm=norm,
hold=False, margins=None, sub=None)
fig = pl.gcf()
if(bool(kwargs.get("save", False))):
fig.savefig(str(kwargs.get("save")), dpi=None, facecolor="none", edgecolor="none", orientation="portrait",
papertype=None, format=None, transparent=False, bbox_inches=None, pad_inches=0.1)
if (bool(kwargs.get("save", False))):
fig.savefig(str(kwargs.get("save")), dpi=None, facecolor="none",
edgecolor="none", orientation="portrait",
papertype=None, format=None, transparent=False,
bbox_inches=None, pad_inches=0.1)
pl.close(fig)
else:
fig.canvas.draw()
......@@ -6,17 +6,22 @@ from nifty.spaces.space import SpaceParadict
class HPSpaceParadict(SpaceParadict):
def __init__(self, nside):
SpaceParadict.__init__(self, nside=nside)
def __init__(self, nside, distances):
if not hasattr(self, 'parameters'):
self.parameters = {}
SpaceParadict.__init__(self, nside=nside, distances=distances)
def __setitem__(self, key, arg):
if key not in ['nside']:
if key not in ['nside', 'distances']:
raise ValueError(about._errors.cstring(
"ERROR: Unsupported hp_space parameter"))
temp = int(arg)
# if(not hp.isnsideok(nside)):
if ((temp & (temp - 1)) != 0) or (temp < 2):
raise ValueError(about._errors.cstring(
"ERROR: invalid parameter ( nside <> 2**n )."))
if key == 'nside':
temp = int(arg)
# if(not hp.isnsideok(nside)):
if ((temp & (temp - 1)) != 0) or (temp < 2):
raise ValueError("ERROR: invalid parameter ( nside <> 2**n ).")
elif key == 'distances':
temp = (arg,)
self.parameters.__setitem__(key, temp)
......@@ -149,11 +149,12 @@ class RGSpace(Space):
None
"""
super(RGSpace, self).__init__(dtype=np.dtype('float64'))
self.paradict = RGSpaceParadict(shape=shape,
zerocenter=zerocenter,
distances=distances,
harmonic=harmonic)
self.dtype = np.dtype(dtype)
@property
def harmonic(self):
......@@ -242,98 +243,9 @@ class RGSpace(Space):
dists = np.sqrt(dists)
return dists
def smooth(self, x, sigma=0, codomain=None, axes=None):
"""
Smoothes an array of field values by convolution with a Gaussian
kernel.
Parameters
----------
x : numpy.ndarray
Array of field values to be smoothed.
sigma : float, *optional*
Standard deviation of the Gaussian kernel, specified in units
of length in position space; for testing: a sigma of -1 will be
reset to a reasonable value (default: 0).
axes: None, tuple
Axes which should be smoothed
Returns
-------
Gx : numpy.ndarray
Smoothed array.
"""
# Check sigma
if sigma == 0:
return x.copy()
elif sigma == -1:
about.infos.cprint(
"INFO: Resetting sigma to sqrt(2)*max(dist).")
sigma = np.sqrt(2) * np.max(self.distances)
elif(sigma < 0):
raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
# if a codomain was given...
if codomain is not None:
# ...check if it was suitable
if not self.check_codomain(codomain):
raise ValueError(about._errors.cstring(
"ERROR: the given codomain is not a compatible!"))
else:
codomain = self.get_codomain()
# TODO: Use the Fourier Transformation Operator for the switch into
# hormonic space.
raise NotImplementedError
x = self.calc_transform(x, codomain=codomain, axes=axes)
x = codomain._calc_smooth_helper(x, sigma, axes=axes)
x = codomain.calc_transform(x, codomain=self, axes=axes)
return x
def _calc_smooth_helper(self, x, sigma, axes=None):
# multiply the gaussian kernel, etc...
if axes is None:
axes = range(len(x.shape))
# if x is hermitian it remains hermitian during smoothing
# TODO look at this later
# if self.datamodel in RG_DISTRIBUTION_STRATEGIES:
remember_hermitianQ = x.hermitian
# Define the Gaussian kernel function
gaussian = lambda x: np.exp(-2. * np.pi**2 * x**2 * sigma**2)
# Define the variables in the dialect of the legacy smoothing.py
nx = np.array(self.shape)
dx = 1 / nx / self.distances
# Multiply the data along each axis with suitable the gaussian kernel
for i in range(len(nx)):
# Prepare the exponent
dk = 1. / nx[i] / dx[i]
nk = nx[i]
k = -0.5 * nk * dk + np.arange(nk) * dk
if self.paradict['zerocenter'][i] == False:
k = np.fft.fftshift(k)
# compute the actual kernel vector
gaussian_kernel_vector = gaussian(k)
# blow up the vector to an array of shape (1,.,1,len(nk),1,.,1)
blown_up_shape = [1, ] * len(x.shape)
blown_up_shape[axes[i]] = len(gaussian_kernel_vector)
gaussian_kernel_vector =\
gaussian_kernel_vector.reshape(blown_up_shape)
# apply the blown-up gaussian_kernel_vector
x = x*gaussian_kernel_vector
try:
x.hermitian = remember_hermitianQ
except AttributeError:
pass
return x
def get_plot(self,x,title="",vmin=None,vmax=None,power=None,unit="",
norm=None,cmap=None,cbar=True,other=None,legend=False,mono=True,**kwargs):
norm=None,cmap=None,cbar=True,other=None,legend=False,
mono=True,**kwargs):
"""
Creates a plot of field values according to the specifications
given by the parameters.
......@@ -396,7 +308,8 @@ class RGSpace(Space):