Commit e94b3f41 authored by Ultima's avatar Ultima
Browse files

Migration of point_space to d2o completed.

Reworked init of d2o.
Implemented a distributor cache.
parent a43a3b04
......@@ -29,7 +29,8 @@ from nifty_cmaps import ncmap
from nifty_core import space,\
point_space,\
nested_space,\
field
field
from nifty_mpi_data import distributed_data_object
from nifty_power import *
from nifty_random import random
......
......@@ -36,7 +36,7 @@ from nifty import * # version 0.8.0
about.warnings.off()
# some signal space; e.g., a two-dimensional regular grid
x_space = rg_space([128, 128]) # define signal space
x_space = rg_space([1280, 1280], datamodel = 'd2o') # define signal space
#x_space = rg_space(512)
#x_space = hp_space(32)
#x_space = gl_space(96)
......@@ -56,6 +56,7 @@ d_space = R.target # get data spac
N = diagonal_operator(d_space, diag=s.var(), bare=True) # define noise covariance
n = N.get_random_field(domain=d_space) # generate noise
d = R(s) + n # compute data
j = R.adjoint_times(N.inverse_times(d)) # define information source
......
......@@ -128,7 +128,7 @@ class lm_space(point_space):
vol : numpy.ndarray
Pixel volume of the :py:class:`lm_space`, which is always 1.
"""
def __init__(self,lmax,mmax=None,datatype=None):
def __init__(self, lmax, mmax=None, datatype=None):
"""
Sets the attributes for an lm_space class instance.
......@@ -186,7 +186,10 @@ class lm_space(point_space):
self.paradict['lmax'] = x[0]
self.paradict['mmax'] = x[1]
def copy(self):
return lm_space(lmax = self.paradict['lmax'],
mmax = self.paradict['mmax'],
datatype = self.datatype)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -213,12 +216,12 @@ class lm_space(point_space):
"""
return self.paradict['mmax']
def shape(self):
def get_shape(self):
mmax = self.paradict['mmax']
lmax = self.paradict['lmax']
return np.array([(mmax+1)*(lmax+1)-((lmax+1)*lmax)//2], dtype=int)
def dim(self,split=False):
def get_dim(self, split=False):
"""
Computes the dimension of the space, i.e.\ the number of spherical
harmonics components that are stored.
......@@ -242,15 +245,15 @@ class lm_space(point_space):
"""
## dim = (mmax+1)*(lmax-mmax/2+1)
if(split):
return self.shape()
return self.get_shape()
#return np.array([(self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2],dtype=np.int)
else:
return np.prod(self.shape())
return np.prod(self.get_shape())
#return (self.para[0]+1)*(self.para[1]+1)-(self.para[1]+1)*self.para[1]//2
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def dof(self):
def get_dof(self):
"""
Computes the number of degrees of freedom of the space, taking into
account symmetry constraints and complex-valuedness.
......@@ -325,7 +328,7 @@ class lm_space(point_space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _getlm(self): ## > compute all (l,m)
index = np.arange(self.dim(split=False))
index = np.arange(self.get_dim(split=False))
n = 2*self.para[0]+1
m = np.ceil((n-np.sqrt(n**2-8*(index-self.para[0])))/2).astype(np.int)
l = index-self.para[0]*m+m*(m-1)//2
......@@ -397,7 +400,7 @@ class lm_space(point_space):
else:
if(np.size(x)==1):
if(extend):
x = self.datatype(x)*np.ones(self.dim(split=True),dtype=self.datatype,order='C')
x = self.datatype(x)*np.ones(self.get_dim(split=True),dtype=self.datatype,order='C')
else:
if(np.isscalar(x)):
x = np.array([x],dtype=self.datatype)
......@@ -459,13 +462,13 @@ class lm_space(point_space):
arg = random.parse_arguments(self,**kwargs)
if(arg is None):
return np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
return np.zeros(self.get_dim(split=True),dtype=self.datatype,order='C')
elif(arg[0]=="pm1"):
x = random.pm1(datatype=self.datatype,shape=self.dim(split=True))
x = random.pm1(datatype=self.datatype,shape=self.get_dim(split=True))
elif(arg[0]=="gau"):
x = random.gau(datatype=self.datatype,shape=self.dim(split=True),mean=None,dev=arg[2],var=arg[3])
x = random.gau(datatype=self.datatype,shape=self.get_dim(split=True),mean=None,dev=arg[2],var=arg[3])
elif(arg[0]=="syn"):
lmax = self.para[0] ## lmax
......@@ -483,7 +486,7 @@ class lm_space(point_space):
return x
elif(arg[0]=="uni"):
x = random.uni(datatype=self.datatype,shape=self.dim(split=True),vmin=arg[1],vmax=arg[2])
x = random.uni(datatype=self.datatype,shape=self.get_dim(split=True),vmin=arg[1],vmax=arg[2])
else:
raise KeyError(about._errors.cstring("ERROR: unsupported random key '"+str(arg[0])+"'."))
......@@ -612,9 +615,9 @@ class lm_space(point_space):
each determine another component with negative :math:`m`.
"""
if(total):
return self.dof()
return self.get_dof()
else:
mol = np.ones(self.dim(split=True),dtype=self.vol.dtype,order='C')
mol = np.ones(self.get_dim(split=True),dtype=self.vol.dtype,order='C')
mol[self.para[0]+1:] = 2 ## redundant in (l,m) and (l,-m)
return mol
......@@ -1060,6 +1063,10 @@ class gl_space(point_space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def copy(self):
return gl_space(nlat = self.paradict['nlat'],
nlon = self.paradict['nlon'],
datatype = self.datatype)
def nlat(self):
"""
Returns the number of latitudinal bins.
......@@ -1082,10 +1089,10 @@ class gl_space(point_space):
"""
return self.paradict['nlon']
def shape(self):
def get_shape(self):
return np.array([(self.paradict['nlat']*self.paradict['nlon'])], dtype=np.int)
def dim(self,split=False):
def get_dim(self,split=False):
"""
Computes the dimension of the space, i.e.\ the number of pixels.
......@@ -1102,15 +1109,15 @@ class gl_space(point_space):
"""
## dim = nlat*nlon
if(split):
return self.shape()
return self.get_shape()
#return np.array([self.para[0]*self.para[1]],dtype=np.int)
else:
return np.prod(self.shape())
return np.prod(self.get_shape())
#return self.para[0]*self.para[1]
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def dof(self):
def get_dof(self):
"""
Computes the number of degrees of freedom of the space.
......@@ -1236,13 +1243,13 @@ class gl_space(point_space):
arg = random.parse_arguments(self,**kwargs)
if(arg is None):
x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
x = np.zeros(self.get_dim(split=True),dtype=self.datatype,order='C')
elif(arg[0]=="pm1"):
x = random.pm1(datatype=self.datatype,shape=self.dim(split=True))
x = random.pm1(datatype=self.datatype,shape=self.get_dim(split=True))
elif(arg[0]=="gau"):
x = random.gau(datatype=self.datatype,shape=self.dim(split=True),mean=None,dev=arg[2],var=arg[3])
x = random.gau(datatype=self.datatype,shape=self.get_dim(split=True),mean=None,dev=arg[2],var=arg[3])
elif(arg[0]=="syn"):
lmax = self.para[0]-1 ## nlat-1
......@@ -1255,7 +1262,7 @@ class gl_space(point_space):
x = self.calc_weight(x,power=0.5)
elif(arg[0]=="uni"):
x = random.uni(datatype=self.datatype,shape=self.dim(split=True),vmin=arg[1],vmax=arg[2])
x = random.uni(datatype=self.datatype,shape=self.get_dim(split=True),vmin=arg[1],vmax=arg[2])
else:
raise KeyError(about._errors.cstring("ERROR: unsupported random key '"+str(arg[0])+"'."))
......@@ -1345,7 +1352,7 @@ class gl_space(point_space):
if(total):
return self.datatype(4*pi)
else:
mol = np.ones(self.dim(split=True),dtype=self.datatype,order='C')
mol = np.ones(self.get_dim(split=True),dtype=self.datatype,order='C')
return self.calc_weight(mol,power=1)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -1689,7 +1696,7 @@ class hp_space(point_space):
"""
niter = 0 ## default number of iterations used for transformations
def __init__(self,nside):
def __init__(self, nside):
"""
Sets the attributes for a hp_space class instance.
......@@ -1734,6 +1741,9 @@ class hp_space(point_space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def copy(self):
return hp_space(nside = self.paradict['nside'])
def nside(self):
"""
Returns the resolution parameter.
......@@ -1745,10 +1755,10 @@ class hp_space(point_space):
"""
return self.paradict['nside']
def shape(self):
def get_shape(self):
return np.array([12*self.paradict['nside']**2], dtype=np.int)
def dim(self,split=False):
def get_dim(self,split=False):
"""
Computes the dimension of the space, i.e.\ the number of pixels.
......@@ -1765,15 +1775,15 @@ class hp_space(point_space):
"""
## dim = 12*nside**2
if(split):
return self.shape()
return self.get_shape()
#return np.array([12*self.para[0]**2],dtype=np.int)
else:
return np.prod(self.shape())
return np.prod(self.get_shape())
#return 12*self.para[0]**2
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def dof(self):
def get_dof(self):
"""
Computes the number of degrees of freedom of the space.
......@@ -1899,13 +1909,13 @@ class hp_space(point_space):
arg = random.parse_arguments(self,**kwargs)
if(arg is None):
x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
x = np.zeros(self.get_dim(split=True),dtype=self.datatype,order='C')
elif(arg[0]=="pm1"):
x = random.pm1(datatype=self.datatype,shape=self.dim(split=True))
x = random.pm1(datatype=self.datatype,shape=self.get_dim(split=True))
elif(arg[0]=="gau"):
x = random.gau(datatype=self.datatype,shape=self.dim(split=True),mean=None,dev=arg[2],var=arg[3])
x = random.gau(datatype=self.datatype,shape=self.get_dim(split=True),mean=None,dev=arg[2],var=arg[3])
elif(arg[0]=="syn"):
lmax = 3*self.para[0]-1 ## 3*nside-1
......@@ -1915,7 +1925,7 @@ class hp_space(point_space):
x = self.calc_weight(x,power=0.5)
elif(arg[0]=="uni"):
x = random.uni(datatype=self.datatype,shape=self.dim(split=True),vmin=arg[1],vmax=arg[2])
x = random.uni(datatype=self.datatype,shape=self.get_dim(split=True),vmin=arg[1],vmax=arg[2])
else:
raise KeyError(about._errors.cstring("ERROR: unsupported random key '"+str(arg[0])+"'."))
......@@ -2001,7 +2011,7 @@ class hp_space(point_space):
if(total):
return self.datatype(4*pi)
else:
mol = np.ones(self.dim(split=True),dtype=self.datatype,order='C')
mol = np.ones(self.get_dim(split=True),dtype=self.datatype,order='C')
return self.calc_weight(mol,power=1)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
......@@ -149,6 +149,8 @@ from nifty_paradict import space_paradict,\
from nifty_about import about
from nifty_random import random
from nifty.nifty_mpi_data import distributed_data_object
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
......@@ -214,7 +216,7 @@ class space(object):
An array of pixel volumes, only one component if the pixels all
have the same volume.
"""
def __init__(self,para=0,datatype=None):
def __init__(self, para=0, datatype=None):
"""
Sets the attributes for a space class instance.
......@@ -260,7 +262,9 @@ class space(object):
"""
return frozenset(dictionary.items())
def copy(self):
return space(para = self.para,
datatype = self.datatype)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def getitem(self, data, key):
......@@ -289,16 +293,16 @@ class space(object):
"ERROR: no generic instance method 'binary_operation'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def norm(self, x, q):
def get_norm(self, x, q):
raise NotImplementedError(about._errors.cstring(\
"ERROR: no generic instance method 'norm'."))
def shape(self):
def get_shape(self):
raise NotImplementedError(about._errors.cstring(\
"ERROR: no generic instance method 'shape'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def dim(self,split=False):
def get_dim(self,split=False):
"""
Computes the dimension of the space, i.e.\ the number of pixels.
......@@ -313,11 +317,12 @@ class space(object):
dim : {int, numpy.ndarray}
Dimension(s) of the space.
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'dim'."))
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'dim'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def dof(self):
def get_dof(self):
"""
Computes the number of degrees of freedom of the space.
......@@ -326,7 +331,8 @@ class space(object):
dof : int
Number of degrees of freedom of the space.
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'dof'."))
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'dof'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -371,7 +377,8 @@ class space(object):
(default: 0).
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'enforce_power'."))
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'enforce_power'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -404,7 +411,8 @@ class space(object):
get_power_indices
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'set_power_indices'."))
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'set_power_indices'."))
def get_power_indices(self,**kwargs):
"""
......@@ -456,7 +464,8 @@ class space(object):
set_power_indices
"""
raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'get_power_indices'."))
raise NotImplementedError(about._errors.cstring(
"ERROR: no generic instance method 'get_power_indices'."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -1059,7 +1068,7 @@ class point_space(space):
vol : numpy.ndarray
Pixel volume of the :py:class:`point_space`, which is always 1.
"""
def __init__(self,num,datatype=None):
def __init__(self, num, datatype=None, datamodel='d2o'):
"""
Sets the attributes for a point_space class instance.
......@@ -1077,15 +1086,30 @@ class point_space(space):
self.paradict = point_space_paradict(num=num)
## check datatype
if(datatype is None):
if (datatype is None):
datatype = np.float64
elif(datatype not in [np.int8,np.int16,np.int32,np.int64,np.float16,np.float32,np.float64,np.complex64,np.complex128]):
elif (datatype not in [np.int8,
np.int16,
np.int32,
np.int64,
np.float16,
np.float32,
np.float64,
np.complex64,
np.complex128]):
about.warnings.cprint("WARNING: data type set to default.")
datatype = np.float64
self.datatype = datatype
if datamodel not in ['np', 'd2o']:
about.warnings.cprint("WARNING: datamodel set to default.")
self.datamodel = 'd2o'
else:
self.datamodel = datamodel
self.discrete = True
self.vol = np.real(np.array([1],dtype=self.datatype))
self.vol = np.real(np.array([1], dtype=self.datatype))
@property
......@@ -1096,7 +1120,13 @@ class point_space(space):
@para.setter
def para(self, x):
self.paradict['num'] = x
def copy(self):
return point_space(num = self.paradict['num'],
datatype = self.datatype,
datamodel = self.datamodel)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def getitem(self, data, key):
return data[key]
......@@ -1104,21 +1134,29 @@ class point_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def setitem(self, data, update, key):
data[key]=update
data[key] = update
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def apply_scalar_function(self, x, function, inplace=False):
if inplace == False:
try:
return function(x)
except:
return np.vectorize(function)(x)
if self.datamodel == 'np':
if inplace == False:
try:
return function(x)
except:
return np.vectorize(function)(x)
else:
try:
x[:] = function(x)
except:
x[:] = np.vectorize(function)(x)
return x
elif self.datamodel == 'd2o':
return x.apply_scalar_function(function, inplace=inplace)
else:
try:
x[:] = function(x)
except:
x[:] = np.vectorize(function)(x)
return x
raise NotImplementedError(about._errors.cstring(
"ERROR: function is not implemented for given datamodel."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -1128,44 +1166,67 @@ class point_space(space):
Valid operations are
"""
def _argmin(z, **kwargs):
ind = np.argmin(z, **kwargs)
if np.isscalar(ind):
ind = np.unravel_index(ind, z.shape, order='C')
if(len(ind)==1):
return ind[0]
return ind
def _argmax(z, **kwargs):
ind = np.argmax(z, **kwargs)
if np.isscalar(ind):
ind = np.unravel_index(ind, z.shape, order='C')
if(len(ind)==1):
return ind[0]
return ind
translation = {"pos" : lambda y: getattr(y, '__pos__')(),
if self.datamodel == 'np':
def _argmin(z, **kwargs):
ind = np.argmin(z, **kwargs)
if np.isscalar(ind):
ind = np.unravel_index(ind, z.shape, order='C')
if(len(ind)==1):
return ind[0]
return ind
def _argmax(z, **kwargs):
ind = np.argmax(z, **kwargs)
if np.isscalar(ind):
ind = np.unravel_index(ind, z.shape, order='C')
if(len(ind)==1):
return ind[0]
return ind
translation = {"pos" : lambda y: getattr(y, '__pos__')(),
"neg" : lambda y: getattr(y, '__neg__')(),
"abs" : lambda y: getattr(y, '__abs__')(),
"nanmin" : np.nanmin,
"min" : np.amin,
"nanmax" : np.nanmax,
"max" : np.amax,
"med" : np.median,
"mean" : np.mean,
"std" : np.std,
"var" : np.var,
"argmin" : _argmin,
"argmin_flat" : np.argmin,
"argmax" : _argmax,
"argmax_flat" : np.argmax,
"conjugate" : np.conjugate,
"sum" : np.sum,
"prod" : np.prod,
"None" : lambda y: y}
elif self.datamodel == 'd2o':
translation = {"pos" : lambda y: getattr(y, '__pos__')(),
"neg" : lambda y: getattr(y, '__neg__')(),
"abs" : lambda y: getattr(y, '__abs__')(),
"nanmin" : np.nanmin,
"min" : np.amin,
"nanmax" : np.nanmax,
"max" : np.amax,
"med" : np.median,
"mean" : np.mean,
"std" : np.std,
"var" : np.var,
"argmin" : _argmin,
"argmin_flat" : np.argmin,
"argmax" : _argmax,
"argmax_flat" : np.argmax,
"conjugate" : np.conjugate,
"sum" : np.sum,
"prod" : np.prod,
"nanmin" : lambda y: getattr(y, 'nanmin')(),
"min" : lambda y: getattr(y, 'amin')(),
"nanmax" : lambda y: getattr(y, 'nanmax')(),
"max" : lambda y: getattr(y, 'amax')(),
"median" : lambda y: getattr(y, 'median')(),
"mean" : lambda y: getattr(y, 'mean')(),
"std" : lambda y: getattr(y, 'std')(),
"var" : lambda y: getattr(y, 'var')(),
"argmin" : lambda y: getattr(y, 'argmin')(),
"argmin_flat" : lambda y: getattr(y, 'argmin_flat')(),
"argmax" : lambda y: getattr(y, 'argmax')(),
"argmax_flat" : lambda y: getattr(y, 'argmax_flat')(),
"conjugate" : lambda y: getattr(y, 'conjugate')(),
"sum" : lambda y: getattr(y, 'sum')(),
"prod" : lambda y: getattr(y, 'prod')(),
"None" : lambda y: y}
else:
raise NotImplementedError(about._errors.cstring(
"ERROR: function is not implemented for given datamodel."))
return translation[op](x, **kwargs)
......@@ -1216,7 +1277,7 @@ class point_space(space):
"""
if(q == 2):
if q == 2:
result = self.calc_dot(x,x)
else:
y = x**(q-1)
......@@ -1228,7 +1289,7 @@ class point_space(space):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def num(self):
def get_num(self):