Skip to content

GitLab

Commit 8bc3a2d0 authored by Ultima's avatar Ultima
Browse files
Options

More tests for lm_space. 

Some bugfixes.
parent e10b54aa
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 688 additions and 1849 deletions
__init__.py
View file @ 8bc3a2d0
......@@ -38,7 +38,6 @@ from nifty_core import space,\
field
from nifty_mpi_data import distributed_data_object, d2o_librarian
from nifty_power import *
from nifty_random import random
from nifty_simple_math import *
from nifty_utilities import *
......@@ -46,7 +45,6 @@ from nifty_utilities import *
from nifty_paradict import space_paradict,\
point_space_paradict,\
nested_space_paradict
from smoothing import *
from operators import *
## optional submodule `rg`
......
demos/demo_wf1.py
View file @ 8bc3a2d0
......@@ -32,25 +32,33 @@
"""
from __future__ import division
import matplotlib as mpl
mpl.use('Agg')
import imp
nifty = imp.load_module('nifty', None,
'/home/steininger/Downloads/nifty', ('','',5))
from nifty import * # version 0.8.0
about.warnings.off()
# some signal space; e.g., a two-dimensional regular grid
x_space = rg_space([1280, 1280]) # define signal space
#x_space = rg_space(512)
x_space = rg_space([4280, 1280]) # define signal space
#y_space = point_space(1280*1280)
#x_space = hp_space(32)
#x_space = gl_space(96)
#x_space = gl_space(800)
k_space = x_space.get_codomain() # get conjugate space
y_space = point_space(1280*1280)
# some power spectrum
power = (lambda k: 42 / (k + 1) ** 3)
power = (lambda k: 42 / (k + 1) ** 4)
S = power_operator(k_space, spec=power) # define signal covariance
S = power_operator(k_space, codomain=x_space, spec=power) # define signal covariance
s = S.get_random_field(domain=x_space) # generate signal
R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None, target = y_space) # define response
R = response_operator(x_space, sigma=0.0, mask=1.0, assign=None) # define response
d_space = R.target # get data space
# some noise variance; e.g., signal-to-noise ratio of 1
......@@ -58,15 +66,15 @@ N = diagonal_operator(d_space, diag=s.var(), bare=True) # define noise
n = N.get_random_field(domain=d_space) # generate noise
d = R(s) #+ n # compute data
d = R(s) + n # compute data
j = R.adjoint_times(N.inverse_times(d)) # define information source
D = propagator_operator(S=S, N=N, R=R) # define information propagator
m = D(j, W=S, tol=1E-2, note=True) # reconstruct map
s.plot(title="signal", save = 'plot_s.png') # plot signal
d_ = field(x_space, val=d.val, target=k_space)
d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png') # plot data
m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png') # plot map
m = D(j, W=S, tol=1E-8, limii=10, note=True) # reconstruct map
#s.plot(title="signal", save = 'plot_s.png') # plot signal
#d_ = field(x_space, val=d.val, target=k_space)
#d_.plot(title="data", vmin=s.min(), vmax=s.max(), save = 'plot_d.png') # plot data
#m.plot(title="reconstructed map", vmin=s.min(), vmax=s.max(), save = 'plot_m.png') # plot map
#
dummys/MPI_dummy.py
View file @ 8bc3a2d0
......@@ -21,7 +21,8 @@ class Comm(object):
class Intracomm(Comm):
def __init__(self):
def __init__(self, name):
self.name = name
self.rank = 0
self.size = 1
......@@ -115,4 +116,4 @@ COMPLEX = _datatype("MPI_COMPLEX")
DOUBLE_COMPLEX = _datatype("MPI_DOUBLE_COMPLEX")
COMM_WORLD = Intracomm()
COMM_WORLD = Intracomm('MPI_dummy_COMM_WORLD')
keepers/nifty_configuration.py
View file @ 8bc3a2d0
......@@ -71,7 +71,7 @@ class configuration(object):
self.set_path(path=path, path_section=path_section)
try:
self.load()
except ValueError:
except:
pass
def __getitem__(self, key):
......
keepers/nifty_default_config.py
View file @ 8bc3a2d0
......@@ -23,13 +23,28 @@ variable_fft_module = variable('fft_module',
variable_lm2gl = variable('lm2gl',
[True, False],
lambda z: z is True or z is False,
genus = 'boolean')
lambda z: isinstance(z, bool),
genus='boolean')
variable_use_healpy = variable(
'use_healpy',
[True, False],
lambda z: (('healpy' in global_dependency_injector)
if z else True) and isinstance(z, bool),
genus='boolean')
variable_use_libsharp = variable('use_libsharp',
[True, False],
lambda z: (('libsharp_wrapper_gl' in
global_dependency_injector)
if z else True) and
isinstance(z, bool),
genus='boolean')
variable_verbosity = variable('verbosity',
[1],
lambda z: z == abs(int(z)),
genus = 'int')
genus='int')
variable_mpi_module = variable('mpi_module',
['MPI', 'MPI_dummy'],
......@@ -45,6 +60,8 @@ variable_default_distribution_strategy = variable(
global_configuration = configuration(
[variable_fft_module,
variable_lm2gl,
variable_use_healpy,
variable_use_libsharp,
variable_verbosity,
variable_mpi_module,
variable_default_distribution_strategy
......@@ -60,4 +77,7 @@ variable_default_comm = variable(
global_configuration.register(variable_default_comm)
try:
global_configuration.load()
except:
pass
keepers/nifty_dependency_injector.py
View file @ 8bc3a2d0
# -*- coding: utf-8 -*-
import imp
import sys
class dependency_injector(object):
......@@ -40,7 +39,7 @@ class dependency_injector(object):
# print pathname
# loaded_module = \
# imp.load_module(module_name, fp, pathname, description)
#
#
# print loaded_module
self.registry[key_name] = loaded_module
except ImportError:
......@@ -52,10 +51,15 @@ class dependency_injector(object):
# except (UnboundLocalError, AttributeError):
# pass
def unregister(self, module_name):
try:
del self.registry['module_name']
except KeyError:
pass
def recursive_import(name):
m = __import__(name)
for n in name.split(".")[1:]:
m = getattr(m, n)
return m
lm/nifty_lm.py
View file @ 8bc3a2d0
......@@ -55,10 +55,6 @@ from nifty.nifty_random import random
gl = gdi['libsharp_wrapper_gl']
hp = gdi['healpy']
if gl is None and gc['lm2gl']:
about.warnings.cprint("WARNING: global setting 'about.lm2gl' corrected.")
gc['lm2gl'] = False
LM_DISTRIBUTION_STRATEGIES = []
GL_DISTRIBUTION_STRATEGIES = []
HP_DISTRIBUTION_STRATEGIES = []
......@@ -154,9 +150,9 @@ class lm_space(point_space):
"""
# check imports
if 'libsharp_wrapper_gl' not in gdi and 'healpy' not in gdi:
if not gc['use_libsharp'] and not gc['use_healpy']:
raise ImportError(about._errors.cstring(
"ERROR: neither libsharp_wrapper_gl nor healpy available."))
"ERROR: neither libsharp_wrapper_gl nor healpy activated."))
self.paradict = lm_space_paradict(lmax=lmax, mmax=mmax)
......@@ -214,9 +210,9 @@ class lm_space(point_space):
dtype=self.dtype)
def get_shape(self):
mmax = self.paradict['mmax']
lmax = self.paradict['lmax']
return (np.int((mmax + 1) * (lmax + 1) - ((lmax + 1) * lmax) // 2),)
mmax = self.paradict['mmax']
return (np.int((mmax + 1) * (lmax + 1) - ((mmax + 1) * mmax) // 2),)
def get_dof(self, split=False):
"""
......@@ -234,8 +230,8 @@ class lm_space(point_space):
symmetry, which is assumed for the spherical harmonics components.
"""
# dof = 2*dim-(lmax+1) = (lmax+1)*(2*mmax+1)*(mmax+1)*mmax
mmax = self.paradict['lmax']
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
dof = np.int((lmax + 1) * (2 * mmax + 1) - (mmax + 1) * mmax)
if split:
return (dof, )
......@@ -284,7 +280,7 @@ class lm_space(point_space):
**kwargs)
complexity_mask = np.iscomplex(casted_x[:self.paradict['lmax']+1])
if np.any(complexity_mask):
about.warnings.cprint("WARNING: Taking only the real parts for " +
about.warnings.cprint("WARNING: Taking the absolute values for " +
"all complex entries where lmax==0")
casted_x[complexity_mask] = np.abs(casted_x[complexity_mask])
return casted_x
......@@ -449,44 +445,43 @@ class lm_space(point_space):
if arg is None:
return np.zeros(self.get_shape(), dtype=self.dtype)
elif arg[0] == "pm1":
elif arg['random'] == "pm1":
x = random.pm1(dtype=self.dtype, shape=self.get_shape())
return self.cast(x)
elif arg[0] == "gau":
elif arg['random'] == "gau":
x = random.gau(dtype=self.dtype,
shape=self.get_shape(),
mean=arg[1],
dev=arg[2],
var=arg[3])
mean=arg['mean'],
std=arg['std'])
return self.cast(x)
elif arg[0] == "syn":
elif arg['random'] == "syn":
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
if self.dtype == np.dtype('complex64'):
if 'libsharp_wrapper_gl' in gdi:
x = gl.synalm_f(arg[1], lmax=lmax, mmax=mmax)
if gc['use_libsharp']:
x = gl.synalm_f(arg['spec'], lmax=lmax, mmax=mmax)
else:
x = hp.synalm(arg[1].astype(np.complex128),
x = hp.synalm(arg['spec'].astype(np.complex128),
lmax=lmax, mmax=mmax).astype(np.complex64)
else:
if 'healpy' in gdi:
x = hp.synalm(arg[1], lmax=lmax, mmax=mmax)
if gc['use_healpy']:
x = hp.synalm(arg['spec'], lmax=lmax, mmax=mmax)
else:
x = gl.synalm(arg[1], lmax=lmax, mmax=mmax)
x = gl.synalm(arg['spec'], lmax=lmax, mmax=mmax)
return x
elif arg[0] == "uni":
elif arg['random'] == "uni":
x = random.uni(dtype=self.dtype,
shape=self.get_shape(),
vmin=arg[1],
vmax=arg[2])
vmin=arg['vmin'],
vmax=arg['vmax'])
return self.cast(x)
else:
raise KeyError(about._errors.cstring(
"ERROR: unsupported random key '" + str(arg[0]) + "'."))
"ERROR: unsupported random key '" + str(arg['random']) + "'."))
def calc_dot(self, x, y):
"""
......@@ -508,7 +503,7 @@ class lm_space(point_space):
x = self.cast(x)
y = self.cast(y)
if 'libsharp_wrapper_gl' in gdi:
if gc['use_libsharp']:
lmax = self.paradict['lmax']
mmax = self.paradict['mmax']
if self.dtype == np.dtype('complex64'):
......@@ -616,7 +611,7 @@ class lm_space(point_space):
sigma = np.sqrt(2) * np.pi / (self.paradict['lmax'] + 1)
elif sigma < 0:
raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
if 'healpy' in gdi:
if gc['use_healpy']:
return hp.smoothalm(x, fwhm=0.0, sigma=sigma, invert=False,
pol=True, mmax=self.paradict['mmax'],
verbose=False, inplace=False)
......@@ -646,14 +641,14 @@ class lm_space(point_space):
# power spectrum
if self.dtype == np.dtype('complex64'):
if 'libsharp_wrapper_gl' in gdi:
if gc['use_libsharp']:
return gl.anaalm_f(x, lmax=lmax, mmax=mmax)
else:
return hp.alm2cl(x.astype(np.complex128), alms2=None,
lmax=lmax, mmax=mmax, lmax_out=lmax,
nspec=None).astype(np.float32)
else:
if 'healpy' in gdi:
if gc['use_healpy']:
return hp.alm2cl(x, alms2=None, lmax=lmax, mmax=mmax,
lmax_out=lmax, nspec=None)
else:
......@@ -758,7 +753,7 @@ class lm_space(point_space):
ax0.set_title(title)
else:
x = self._enforce_shape(np.array(x))
x = self.cast(x)
if(np.iscomplexobj(x)):
if(title):
title += " "
......@@ -841,15 +836,6 @@ class lm_space(point_space):
else:
fig.canvas.draw()
def __repr__(self):
return "<nifty_lm.lm_space>"
def __str__(self):
return ("nifty_lm.lm_space instance\n- lmax = " +
str(self.paradict['lmax']) +
"\n- mmax = " + str(self.paradict['mmax']) +
"\n- dtype = " + str(self.dtype))
def getlm(self): # > compute all (l,m)
index = np.arange(self.get_dim())
n = 2 * self.paradict['lmax'] + 1
......@@ -940,9 +926,9 @@ class gl_space(point_space):
"""
# check imports
if 'libsharp_wrapper_gl' not in gdi:
if not gc['use_libsharp']:
raise ImportError(about._errors.cstring(
"ERROR: libsharp_wrapper_gl not available."))
"ERROR: libsharp_wrapper_gl not loaded."))
self.paradict = gl_space_paradict(nlat=nlat, nlon=nlon)
......@@ -1065,12 +1051,12 @@ class gl_space(point_space):
Compatible codomains are instances of :py:class:`gl_space` and
:py:class:`lm_space`.
"""
if not isinstance(codomain, space):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if codomain is None:
return False
if not isinstance(codomain, space):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if self.datamodel is not codomain.datamodel:
return False
......@@ -1158,38 +1144,40 @@ class gl_space(point_space):
if(arg is None):
x = np.zeros(self.get_shape(), dtype=self.dtype)
elif(arg[0] == "pm1"):
elif(arg['random'] == "pm1"):
x = random.pm1(dtype=self.dtype, shape=self.get_shape())
elif(arg[0] == "gau"):
elif(arg['random'] == "gau"):
x = random.gau(dtype=self.dtype,
shape=self.get_shape(),
mean=arg[1], dev=arg[2], var=arg[3])
mean=arg['mean'],
std=arg['std'])
elif(arg[0] == "syn"):
elif(arg['random'] == "syn"):
nlat = self.paradict['nlat']
nlon = self.paradict['nlon']
lmax = nlat - 1
if self.dtype == np.dtype('float32'):
x = gl.synfast_f(arg[1],
x = gl.synfast_f(arg['syn'],
nlat=nlat, nlon=nlon,
lmax=lmax, mmax=lmax, alm=False)
else:
x = gl.synfast(arg[1],
x = gl.synfast(arg['syn'],
nlat=nlat, nlon=nlon,
lmax=lmax, mmax=lmax, alm=False)
# weight if discrete
if self.discrete:
x = self.calc_weight(x, power=0.5)
elif(arg[0] == "uni"):
elif(arg['random'] == "uni"):
x = random.uni(dtype=self.dtype,
shape=self.get_shape(),
vmin=arg[1], vmax=arg[2])
vmin=arg['vmin'],
vmax=arg['vmax'])
else:
raise KeyError(about._errors.cstring(
"ERROR: unsupported random key '" + str(arg[0]) + "'."))
"ERROR: unsupported random key '" + str(arg['random']) + "'."))
return x
......@@ -1455,7 +1443,7 @@ class gl_space(point_space):
ax0.set_title(title)
else:
x = self._enforce_shape(np.array(x, dtype=self.dtype))
x = self.cast(x)
if(vmin is None):
vmin = np.min(x, axis=None, out=None)
if(vmax is None):
......@@ -1503,12 +1491,6 @@ class gl_space(point_space):
else:
fig.canvas.draw()
def __repr__(self):
return "<nifty_lm.gl_space>"
def __str__(self):
return "nifty_lm.gl_space instance\n- nlat = " + str(self.para[0]) + "\n- nlon = " + str(self.para[1]) + "\n- dtype = numpy." + str(np.result_type(self.dtype))
class hp_space(point_space):
"""
......@@ -1583,7 +1565,7 @@ class hp_space(point_space):
"""
# check imports
if 'healpy' not in gdi:
if not gc['use_healpy']:
raise ImportError(about._errors.cstring(
"ERROR: healpy not available."))
......@@ -1697,12 +1679,12 @@ class hp_space(point_space):
Compatible codomains are instances of :py:class:`hp_space` and
:py:class:`lm_space`.
"""
if not isinstance(codomain, space):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if codomain is None:
return False
if not isinstance(codomain, space):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
if self.datamodel is not codomain.datamodel:
return False
......@@ -1781,29 +1763,31 @@ class hp_space(point_space):
if arg is None:
x = np.zeros(self.get_shape(), dtype=self.dtype)
elif arg[0] == "pm1":
elif arg['random'] == "pm1":
x = random.pm1(dtype=self.dtype, shape=self.get_shape())
elif arg[0] == "gau":
elif arg['random'] == "gau":
x = random.gau(dtype=self.dtype, shape=self.get_shape(),
mean=arg[1], dev=arg[2], var=arg[3])
mean=arg['mean'],
std=arg['std'])
elif arg[0] == "syn":
elif arg['random'] == "syn":
nside = self.paradict['nside']
lmax = 3*nside-1
x = hp.synfast(arg[1], nside, lmax=lmax, mmax=lmax, alm=False,
x = hp.synfast(arg['spec'], nside, lmax=lmax, mmax=lmax, alm=False,
pol=True, pixwin=False, fwhm=0.0, sigma=None)
# weight if discrete
if self.discrete:
x = self.calc_weight(x, power=0.5)
elif arg[0] == "uni":
elif arg['random'] == "uni":
x = random.uni(dtype=self.dtype, shape=self.get_shape(),
vmin=arg[1], vmax=arg[2])
vmin=arg['vmin'],
vmax=arg['vmax'])
else:
raise KeyError(about._errors.cstring(
"ERROR: unsupported random key '" + str(arg[0]) + "'."))
"ERROR: unsupported random key '" + str(arg['random']) + "'."))
return x
......@@ -2038,7 +2022,7 @@ class hp_space(point_space):
ax0.set_title(title)
else:
x = self._enforce_shape(np.array(x, dtype=self.dtype))
x = self.cast(x)
if(norm == "log"):
if(vmin is not None):
if(vmin <= 0):
......@@ -2060,11 +2044,3 @@ class hp_space(point_space):
pl.close(fig)
else:
fig.canvas.draw()
def __repr__(self):
return "<nifty_lm.hp_space>"
def __str__(self):
return "nifty_lm.hp_space instance\n- nside = " + str(self.para[0])
nifty_core.py
View file @ 8bc3a2d0
......@@ -407,8 +407,11 @@ class space(object):
check : bool
Whether or not the given codomain is compatible to the space.
"""
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'check_codomain'."))
if codomain is None:
return False
else:
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'check_codomain'."))
def get_codomain(self, **kwargs):
"""
......@@ -735,11 +738,13 @@ class space(object):
"ERROR: no generic instance method 'get_plot'."))
def __repr__(self):
return "<nifty_core.space>"
string = ""
string += str(type(self)) + "\n"
string += "paradict: " + str(self.paradict) + "\n"
return string
def __str__(self):
return "nifty_core.space instance\n- para = " + str(self.para) + \
"\n- dtype = " + str(self.dtype.type)
return self.__repr__()
class point_space(space):
......@@ -1808,13 +1813,16 @@ class point_space(space):
fig.canvas.draw()
def __repr__(self):
return "<nifty_core.point_space>"
def __str__(self):
result = "nifty_core.point_space instance\n- num = " + \
str(self.para[0]) + "\n- dtype = numpy." + \
str(np.dtype(self.dtype))
return result
string = ""
string += str(type(self)) + "\n"
string += "paradict: " + str(self.paradict)
string += 'dtype: ' + str(self.dtype) + "\n"
string += 'datamodel: ' + str(self.datamodel) + "\n"
string += 'comm: ' + self.comm.name + "\n"
string += 'discrete: ' + str(self.discrete) + "\n"
string += 'harmonic: ' + str(self.harmonic) + "\n"
string += 'distances: ' + str(self.distances) + "\n"
return string
class field(object):
......
nifty_mpi_data.py
View file @ 8bc3a2d0
......@@ -281,6 +281,9 @@ class distributed_data_object(object):
def __gt__(self, other):
return self._compare_helper(other, '__gt__')
def __iter__(self):
return self.distributor.get_iter(self)
def equal(self, other):
if other is None:
return False
......@@ -2427,6 +2430,8 @@ class _slicing_distributor(distributor):
raise ImportError(about._errors.cstring(
"ERROR: h5py is not available"))
def get_iter(self, d2o):
return d2o_slicing_iter(d2o)
def _equal_slicer(comm, global_shape):
rank = comm.rank
......@@ -2675,6 +2680,9 @@ class _not_distributor(distributor):
raise ImportError(about._errors.cstring(