Skip to content

GitLab

Commit d65bd374 authored by theos's avatar theos Committed by csongor
Browse files
Options

Updated LMSpace to new Space interface. 

compute_k_array is still missing.
parent 212238a5
Hide whitespace changes
Inline Side-by-side
Showing with 39 additions and 439 deletions
nifty/spaces/lm_space/__init__.py
View file @ d65bd374
# -*- coding: utf-8 -*-
from lm_space import LMSpace
nifty/spaces/lm_space/lm_space.py
View file @ d65bd374
from __future__ import division
import os
import numpy as np
import pylab as pl
from matplotlib.colors import LogNorm as ln
from matplotlib.ticker import LogFormatter as lf
from d2o import STRATEGIES as DISTRIBUTION_STRATEGIES
from nifty.spaces.space import Space
......@@ -15,13 +9,9 @@ from nifty.config import about,\
nifty_configuration as gc,\
dependency_injector as gdi
from nifty.nifty_random import random
gl = gdi.get('libsharp_wrapper_gl')
hp = gdi.get('healpy')
LM_DISTRIBUTION_STRATEGIES = DISTRIBUTION_STRATEGIES['global']
class LMSpace(Space):
"""
......@@ -120,439 +110,45 @@ class LMSpace(Space):
self._lmax = self._parse_lmax(lmax)
self._mmax = self._parse_mmax(mmax)
def _identifier(self):
# Extract the identifying parts from the vars(self) dict.
temp = [(ii[0],
((lambda x: tuple(x) if
isinstance(x, np.ndarray) else x)(ii[1])))
for ii in vars(self).iteritems()
if ii[0] not in ['_cache_dict', 'power_indices']]
# Return the sorted identifiers as a tuple.
return tuple(sorted(temp))
def copy(self):
return self.__class__(lmax=self.lmax,
mmax=self.mmax,
dtype=self.dtype)
def compute_k_array(self, distribution_strategy):
# return a d2o with the l-value at the individual pixels
# e.g. for 0<=l<=2: [0,1,2, 1,1,2,2, 2,2]
pass
@property
def shape(self):
return (np.int((self.mmax + 1) * (self.lmax + 1)),)
@property
def dim(self):
return np.int((self.mmax + 1) * (self.lmax + 1))
# ---Mandatory properties and methods---
@property
def harmonic(self):
return True
@property
def meta_volume(self):
"""
Calculates the meta volumes.
The meta volumes are the volumes associated with each component of
a field, taking into account field components that are not
explicitly included in the array of field values but are determined
by symmetry conditions.
Parameters
----------
total : bool, *optional*
Whether to return the total meta volume of the space or the
individual ones of each field component (default: False).
Returns
-------
mol : {numpy.ndarray, float}
Meta volume of the field components or the complete space.
Notes
-----
The spherical harmonics components with :math:`m=0` have meta
volume 1, the ones with :math:`m>0` have meta volume 2, sinnce they
each determine another component with negative :math:`m`.
"""
return np.float(self.dof())
def shape(self):
return (self.dim, )
@property
def meta_volume_split(self):
mol = self.cast(1, dtype=np.float)
mol[self.paradict['lmax'] + 1:] = 2 # redundant: (l,m) and (l,-m)
return mol
def complement_cast(self, x, axis=None, **kwargs):
if axis is None:
lmax = self.paradict['lmax']
complexity_mask = x[:lmax+1].iscomplex()
if complexity_mask.any():
about.warnings.cprint("WARNING: Taking the absolute values for " +
"all complex entries where lmax==0")
x[:lmax+1] = abs(x[:lmax+1])
else:
# TODO hermitianize only on specific axis
lmax = self.paradict['lmax']
complexity_mask = x[:lmax+1].iscomplex()
if complexity_mask.any():
about.warnings.cprint("WARNING: Taking the absolute values for " +
"all complex entries where lmax==0")
x[:lmax+1] = abs(x[:lmax+1])
return x
# TODO: Extend to binning/log
def enforce_power(self, spec, size=None, kindex=None):
if kindex is None:
kindex_size = self.paradict['lmax'] + 1
kindex = np.arange(kindex_size,
dtype=np.array(self.distances).dtype)
return self._enforce_power_helper(spec=spec,
size=size,
kindex=kindex)
def _check_codomain(self, codomain):
"""
Checks whether a given codomain is compatible to the
:py:class:`lm_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:`lm_space`,
:py:class:`gl_space`, and :py:class:`hp_space`.
"""
if codomain is None:
return False
from nifty.spaces.hp_space import HPSpace
from nifty.spaces.gl_space import GLSpace
if not isinstance(codomain, Space):
raise TypeError(about._errors.cstring(
"ERROR: The given codomain must be a nifty lm_space."))
elif isinstance(codomain, GLSpace):
# lmax==mmax
# nlat==lmax+1
# nlon==2*lmax+1
if ((self.paradict['lmax'] == self.paradict['mmax']) and
(codomain.paradict['nlat'] == self.paradict['lmax']+1) and
(codomain.paradict['nlon'] == 2*self.paradict['lmax']+1)):
return True
elif isinstance(codomain, HPSpace):
# lmax==mmax
# 3*nside-1==lmax
if ((self.paradict['lmax'] == self.paradict['mmax']) and
(3*codomain.paradict['nside']-1 == self.paradict['lmax'])):
return True
return False
def get_random_values(self, **kwargs):
"""
Generates random field values according to the specifications given
by the parameters, taking into account complex-valuedness and
hermitian symmetry.
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} or {+1,+i,-1,-i}
- "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).
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 = 0
#
# 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":
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
if self.dtype == np.dtype('complex64'):
if gc['use_libsharp']:
sample = gl.synalm_f(arg['spec'], lmax=lmax, mmax=mmax)
else:
sample = hp.synalm(
arg['spec'].astype(np.complex128),
lmax=lmax, mmax=mmax).astype(np.complex64,
copy=False)
else:
if gc['use_healpy']:
sample = hp.synalm(arg['spec'], lmax=lmax, mmax=mmax)
else:
sample = gl.synalm(arg['spec'], lmax=lmax, mmax=mmax)
else:
sample = super(LMSpace, 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 dot_contraction(self, x, axes):
assert len(axes) == 1
axis = axes[0]
lmax = self.paradict['lmax']
# extract the low and high parts of x
extractor = ()
extractor += (slice(None),)*axis
low_extractor = extractor + (slice(None, lmax+1), )
high_extractor = extractor + (slice(lmax+1), )
result = x[low_extractor].sum(axes) + 2 * x[high_extractor].sum(axes)
return result
def get_plot(self, x, title="", vmin=None, vmax=None, power=True,
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.
Parameters
----------
x : numpy.ndarray
Array containing the field values.
def dim(self):
l = self.lmax
m = self.mmax
# the LMSpace consist of the full triangle, minus two little triangles
# if mmax < lmax
# dim = l(l+1)/2 - 2 * (l-m)(l-m+1)/2
return np.int(l*(l+1.)/2. - 2.*(l-m)*(l-m+1.)/2.)
Returns
-------
None
Other parameters
----------------
title : string, *optional*
Title of the plot (default: "").
vmin : float, *optional*
Minimum value to be displayed (default: ``min(x)``).
vmax : float, *optional*
Maximum value to be displayed (default: ``max(x)``).
power : bool, *optional*
Whether to plot the power contained in the field or the field
values themselves (default: True).
unit : string, *optional*
Unit of the field values (default: "").
norm : string, *optional*
Scaling of the field values before plotting (default: None).
cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
Color map to be used for two-dimensional plots (default: None).
cbar : bool, *optional*
Whether to show the color bar or not (default: True).
other : {single object, tuple of objects}, *optional*
Object or tuple of objects to be added, where objects can be
scalars, arrays, or fields (default: None).
legend : bool, *optional*
Whether to show the legend or not (default: False).
mono : bool, *optional*
Whether to plot the monopole or not (default: True).
save : string, *optional*
Valid file name where the figure is to be stored, by default
the figure is not saved (default: False).
@property
def total_volume(self):
# the individual pixels have a fixed volume of 1.
return np.float(self.dim)
"""
from nifty.field import Field
try:
x = x.get_full_data()
except AttributeError:
pass
if(not pl.isinteractive())and(not bool(kwargs.get("save", False))):
about.warnings.cprint("WARNING: interactive mode off.")
if(power):
x = self.calc_power(x)
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(self.para[0] + 1, dtype=np.int)
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):
vmax = np.max(x[:mono].tolist(
) + (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)):
other = list(other)
for ii in xrange(len(other)):
if(isinstance(other[ii], Field)):
other[ii] = other[ii].power(**kwargs)
else:
other[ii] = self.enforce_power(other[ii])
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.legend()
ax0.set_xlim(xaxes[1], xaxes[-1])
ax0.set_xlabel(r"$\ell$")
ax0.set_ylim(vmin, vmax)
ax0.set_ylabel(r"$\ell(2\ell+1) C_\ell$")
ax0.set_title(title)
def copy(self):
return self.__class__(lmax=self.lmax,
mmax=self.mmax,
dtype=self.dtype)
def weight(self, x, power=1, axes=None, inplace=False):
if inplace:
return x
else:
if(np.iscomplexobj(x)):
if(title):
title += " "
if(bool(kwargs.get("save", False))):
save_ = os.path.splitext(
os.path.basename(str(kwargs.get("save"))))
kwargs.update(save=save_[0] + "_absolute" + save_[1])
self.get_plot(np.absolute(x), title=title + "(absolute)", vmin=vmin, vmax=vmax,
power=False, norm=norm, cmap=cmap, cbar=cbar, other=None, legend=False, **kwargs)
# self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
if(cmap is None):
cmap = pl.cm.hsv_r
if(bool(kwargs.get("save", False))):
kwargs.update(save=save_[0] + "_phase" + save_[1])
self.get_plot(np.angle(x, deg=False), title=title + "(phase)", vmin=-3.1416, vmax=3.1416, power=False,
norm=None, cmap=cmap, cbar=cbar, other=None, legend=False, **kwargs) # values in [-pi,pi]
return None # leave method
else:
if(vmin is None):
vmin = np.min(x, axis=None, out=None)
if(vmax is None):
vmax = np.max(x, axis=None, out=None)
if(norm == "log")and(vmin <= 0):
raise ValueError(about._errors.cstring(
"ERROR: nonpositive value(s)."))
# not a number
xmesh = np.nan * \
np.empty(self.para[::-1] + 1, dtype=np.float16, order='C')
xmesh[4, 1] = None
xmesh[1, 4] = None
lm = 0
for mm in xrange(self.para[1] + 1):
xmesh[mm][mm:] = x[lm:lm + self.para[0] + 1 - mm]
lm += self.para[0] + 1 - mm
s_ = np.array([1, self.para[1] / self.para[0]
* (1.0 + 0.159 * bool(cbar))])
fig = pl.figure(num=None, figsize=(
6.4 * s_[0], 6.4 * s_[1]), dpi=None, facecolor="none", edgecolor="none", frameon=False, FigureClass=pl.Figure)
ax0 = fig.add_axes(
[0.06 / s_[0], 0.06 / s_[1], 1.0 - 0.12 / s_[0], 1.0 - 0.12 / s_[1]])
ax0.set_axis_bgcolor([0.0, 0.0, 0.0, 0.0])
xaxes = np.arange(self.para[0] + 2, dtype=np.int) - 0.5
yaxes = np.arange(self.para[1] + 2, dtype=np.int) - 0.5
if(norm == "log"):
n_ = ln(vmin=vmin, vmax=vmax)
else:
n_ = None
sub = ax0.pcolormesh(xaxes, yaxes, np.ma.masked_where(np.isnan(
xmesh), xmesh), cmap=cmap, norm=n_, vmin=vmin, vmax=vmax, clim=(vmin, vmax))
ax0.set_xlim(xaxes[0], xaxes[-1])
ax0.set_xticks([0], minor=False)
ax0.set_xlabel(r"$\ell$")
ax0.set_ylim(yaxes[0], yaxes[-1])
ax0.set_yticks([0], minor=False)
ax0.set_ylabel(r"$m$")
ax0.set_aspect("equal")
if(cbar):
if(norm == "log"):
f_ = lf(10, labelOnlyBase=False)
b_ = sub.norm.inverse(
np.linspace(0, 1, sub.cmap.N + 1))
v_ = np.linspace(
sub.norm.vmin, sub.norm.vmax, sub.cmap.N)
else:
f_ = None
b_ = None
v_ = None
fig.colorbar(sub, ax=ax0, orientation="horizontal", fraction=0.1, pad=0.05, shrink=0.75, aspect=20, ticks=[
vmin, vmax], format=f_, drawedges=False, boundaries=b_, values=v_)
ax0.set_title(title)
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()
def getlm(self): # > compute all (l,m)
index = np.arange(self.dim)
n = 2 * self.paradict['lmax'] + 1
m = np.ceil(
(n - np.sqrt(n**2 - 8 * (index - self.paradict['lmax']))) / 2
).astype(np.int)
l = index - self.paradict['lmax'] * m + m * (m - 1) // 2
return l, m
return x.copy()
# ---Added properties and methods---
......@@ -568,7 +164,7 @@ class LMSpace(Space):
lmax = np.int(lmax)
if lmax < 1:
raise ValueError(about._errors.cstring(
"ERROR: lmax: nonpositive number."))
"ERROR: negative lmax is not allowed."))
# exception lmax == 2 (nside == 1)
if (lmax % 2 == 0) and (lmax > 2):
about.warnings.cprint(
......@@ -576,13 +172,18 @@ class LMSpace(Space):
return lmax
def _parse_mmax(self, mmax):
mmax = int(mmax)
if (mmax < 1) or(mmax > self.lmax):
about.warnings.cprint(
"WARNING: mmax parameter set to default.")
if mmax is None:
mmax = self.lmax
if(mmax != self.lmax):
else:
mmax = int(mmax)
if mmax < 1:
raise ValueError(about._errors.cstring(
"ERROR: mmax < 1 is not allowed."))
if mmax > self.lmax:
raise ValueError(about._errors.cstring(
"ERROR: mmax > lmax is not allowed."))
if mmax != self.lmax:
about.warnings.cprint(
"WARNING: unrecommended parameter set (mmax <> lmax).")
"WARNING: unrecommended parameter combination (mmax <> lmax).")
return mmax
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