Commit afca20a4 authored by Ultima's avatar Ultima
Browse files

Finished update power_operator and vecvec_operator.

Fixed a small error in random number generation in nifty_lm.py
parent 4a3fa523
...@@ -456,7 +456,7 @@ class lm_space(point_space): ...@@ -456,7 +456,7 @@ class lm_space(point_space):
vmax : float, *optional* vmax : float, *optional*
Upper limit for a uniform distribution (default: 1). Upper limit for a uniform distribution (default: 1).
""" """
arg = random.arguments(self,**kwargs) arg = random.parse_arguments(self,**kwargs)
if(arg is None): if(arg is None):
return np.zeros(self.dim(split=True),dtype=self.datatype,order='C') return np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
...@@ -1233,7 +1233,7 @@ class gl_space(point_space): ...@@ -1233,7 +1233,7 @@ class gl_space(point_space):
vmax : float, *optional* vmax : float, *optional*
Upper limit for a uniform distribution (default: 1). Upper limit for a uniform distribution (default: 1).
""" """
arg = random.arguments(self,**kwargs) arg = random.parse_arguments(self,**kwargs)
if(arg is None): if(arg is None):
x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C') x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
...@@ -1896,7 +1896,7 @@ class hp_space(point_space): ...@@ -1896,7 +1896,7 @@ class hp_space(point_space):
vmax : float, *optional* vmax : float, *optional*
Upper limit for a uniform distribution (default: 1). Upper limit for a uniform distribution (default: 1).
""" """
arg = random.arguments(self,**kwargs) arg = random.parse_arguments(self,**kwargs)
if(arg is None): if(arg is None):
x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C') x = np.zeros(self.dim(split=True),dtype=self.datatype,order='C')
...@@ -2124,7 +2124,7 @@ class hp_space(point_space): ...@@ -2124,7 +2124,7 @@ class hp_space(point_space):
if(self.discrete): if(self.discrete):
x = self.calc_weight(x,power=-0.5) x = self.calc_weight(x,power=-0.5)
## power spectrum ## power spectrum
return hp.anafast(x,map2=None,nspec=None,lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,iter=kwargs.get("iter",self.niter),alm=False,pol=True,use_weights=False,regression=False,datapath=None) return hp.anafast(x,map2=None,nspec=None,lmax=3*self.para[0]-1,mmax=3*self.para[0]-1,iter=kwargs.get("iter",self.niter),alm=False,pol=True,use_weights=False,datapath=None)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
...@@ -174,11 +174,11 @@ class operator(object): ...@@ -174,11 +174,11 @@ class operator(object):
@property @property
def val(self): def val(self):
return self.__val return self._val
@val.setter @val.setter
def val(self, x): def val(self, x):
self.__val = self.domain.cast(x) self._val = self.domain.cast(x)
def set_val(self, new_val): def set_val(self, new_val):
...@@ -1989,7 +1989,7 @@ class power_operator(diagonal_operator): ...@@ -1989,7 +1989,7 @@ class power_operator(diagonal_operator):
The space wherein the operator output lives The space wherein the operator output lives
""" """
def __init__(self,domain,spec=1,bare=True,pindex=None,**kwargs): def __init__(self, domain, spec=1, bare=True, **kwargs):
""" """
Sets the diagonal operator's standard properties Sets the diagonal operator's standard properties
...@@ -2034,9 +2034,29 @@ class power_operator(diagonal_operator): ...@@ -2034,9 +2034,29 @@ class power_operator(diagonal_operator):
(default: 0). (default: 0).
""" """
if(not isinstance(domain,space)): ## Set the domain
if isinstance(domain,space) == False:
raise TypeError(about._errors.cstring("ERROR: invalid input.")) raise TypeError(about._errors.cstring("ERROR: invalid input."))
self.domain = domain self.domain = domain
## Set the target
self.target = self.domain
## Set imp
self.imp = True
## Save the kwargs
self.kwargs = kwargs
## Set the diag
self.set_power(new_spec = spec, bare = bare, **kwargs)
self.sym = True
## check whether identity
if(np.all(spec==1)):
self.uni = True
else:
self.uni = False
"""
## check implicit pindex ## check implicit pindex
if(pindex is None): if(pindex is None):
try: try:
...@@ -2065,20 +2085,21 @@ class power_operator(diagonal_operator): ...@@ -2065,20 +2085,21 @@ class power_operator(diagonal_operator):
else: else:
self.val = diag self.val = diag
self.sym = True """
@property
## check whether identity def val(self):
if(np.all(spec==1)): return self._val
self.uni = True
else: ## The domain is used for calculations of the power-spectrum, not for
self.uni = False ## actual field values. Therefore the casting of self.val must be switched
## off.
self.imp = True @val.setter
self.target = self.domain def val(self, x):
self._val = x
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_power(self,newspec,bare=True,pindex=None,**kwargs): def set_power(self, new_spec, bare=True, pindex=None, **kwargs):
""" """
Sets the power spectrum of the diagonal operator Sets the power spectrum of the diagonal operator
...@@ -2116,42 +2137,67 @@ class power_operator(diagonal_operator): ...@@ -2116,42 +2137,67 @@ class power_operator(diagonal_operator):
(default: None). (default: None).
""" """
# if(bare is None):
# about.warnings.cprint("WARNING: bare keyword set to default.")
# bare = True ## Cast the pontentially given pindex. If no pindex was given,
## check implicit pindex ## extract it from self.domain using the supplied kwargs.
if(pindex is None): pindex = self._cast_pindex(pindex, **kwargs)
try:
self.domain.set_power_indices(**kwargs)
except: ## Cast the new powerspectrum function
raise ValueError(about._errors.cstring("ERROR: invalid domain.")) temp_spec = self.domain.enforce_power(new_spec)
else:
pindex = self.domain.power_indices.get("pindex") ## Calculate the diagonal
## check explicit pindex
else:
pindex = np.array(pindex,dtype=np.int)
if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
## set diagonal
try: try:
diag = self.domain.enforce_power(newspec,size=np.max(pindex,axis=None,out=None)+1)[pindex] diag = pindex.apply_scalar_function(lambda x: temp_spec[x],
except(AttributeError): dtype = temp_spec.dtype.type)
raise ValueError(about._errors.cstring("ERROR: invalid input.")) diag.hermitian = True
## weight if ... except(AttributeError): ##TODO: update all pindices to d2o's
if(not self.domain.discrete)and(bare): diag = temp_spec[pindex]
self.val = np.real(self.domain.calc_weight(diag,power=1))
## Weight if necessary
if self.domain.discrete == False and bare == True:
self.val = self.domain.calc_weight(diag, power=1)
else: else:
self.val = diag self.val = diag
## check whether identity ## check whether identity
if(np.all(newspec==1)): if (self.val == 1).all() == True:
self.uni = True self.uni = True
else: else:
self.uni = False self.uni = False
return self.val
def _cast_pindex(self, pindex = None, **kwargs):
## Update the internal kwargs dict with the given one:
temp_kwargs = self.kwargs
temp_kwargs.update(kwargs)
## Case 1: no pindex given
if pindex is None:
try:
pindex = self.domain.power_indices.\
get_index_dict(temp_kwargs)['pindex']
except(AttributeError):
## TODO: update all spaces to use the power_indices class
try:
self.domain.set_power_indices(temp_kwargs)
except:
raise ValueError(about._errors.cstring(
"ERROR: Domain is not capable of returning a pindex"))
else:
pindex = self.domain.power_indices.get("pindex")
## Case 2: explicit pindex given
else:
## TODO: Pindex casting could be done here. No must-have.
assert(np.all(pindex.shape == self.domain.shape()))
return pindex
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_power(self,bare=True,pundex=None,pindex=None,**kwargs): def get_power(self, bare=True, **kwargs):
""" """
Computes the power spectrum Computes the power spectrum
...@@ -2192,35 +2238,29 @@ class power_operator(diagonal_operator): ...@@ -2192,35 +2238,29 @@ class power_operator(diagonal_operator):
(default: 0). (default: 0).
""" """
## weight if ...
if(not self.domain.discrete)and(bare): temp_kwargs = self.kwargs
diag = np.real(self.domain.calc_weight(self.val,power=-1)) temp_kwargs.update(kwargs)
## Weight the diagonal values if necessary
if self.domain.discrete == False and bare == True:
diag = self.domain.calc_weight(self.val, power = -1)
else: else:
diag = self.val diag = self.val
## check implicit pundex
if(pundex is None):
if(pindex is None):
try:
self.domain.set_power_indices(**kwargs)
except:
raise ValueError(about._errors.cstring("ERROR: invalid domain."))
else:
pundex = self.domain.power_indices.get("pundex")
else:
pindex = np.array(pindex,dtype=np.int)
if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
## quick pundex
pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
## check explicit pundex
else:
pundex = np.array(pundex,dtype=np.int)
return diag.flatten(order='C')[pundex] ## Use the calc_power routine of the domain in order to to stay
## independent of the implementation
diag = diag**(0.5)
power = self.domain.calc_power(diag, **temp_kwargs)
return power
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_projection_operator(self,pindex=None,**kwargs): def get_projection_operator(self, pindex=None, **kwargs):
""" """
Generates a spectral projection operator Generates a spectral projection operator
...@@ -2253,21 +2293,9 @@ class power_operator(diagonal_operator): ...@@ -2253,21 +2293,9 @@ class power_operator(diagonal_operator):
(default: 0). (default: 0).
""" """
## check implicit pindex
if(pindex is None):
try:
self.domain.set_power_indices(**kwargs)
except:
raise ValueError(about._errors.cstring("ERROR: invalid domain."))
else:
pindex = self.domain.power_indices.get("pindex")
## check explicit pindex
else:
pindex = np.array(pindex,dtype=np.int)
if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
return projection_operator(self.domain,assign=pindex) pindex = self._cast_pindex(pindex, **kwargs)
return projection_operator(self.domain, assign=pindex)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -2366,7 +2394,7 @@ class projection_operator(operator): ...@@ -2366,7 +2394,7 @@ class projection_operator(operator):
The space wherein the operator output lives The space wherein the operator output lives
""" """
def __init__(self,domain,assign=None,**kwargs): def __init__(self, domain, assign=None, **kwargs):
""" """
Sets the standard operator properties and `indexing`. Sets the standard operator properties and `indexing`.
...@@ -2707,7 +2735,7 @@ class vecvec_operator(operator): ...@@ -2707,7 +2735,7 @@ class vecvec_operator(operator):
target : space target : space
The space wherein the operator output lives. The space wherein the operator output lives.
""" """
def __init__(self,domain=None,val=1): def __init__(self, domain=None, val=1):
""" """
Sets the standard operator properties and `values`. Sets the standard operator properties and `values`.
...@@ -2725,20 +2753,21 @@ class vecvec_operator(operator): ...@@ -2725,20 +2753,21 @@ class vecvec_operator(operator):
------- -------
None None
""" """
if(domain is None)and(isinstance(val,field)): if domain is None and isinstance(val,field)==True:
domain = val.domain domain = val.domain
if(not isinstance(domain,space)): if isinstance(domain,space) == False:
raise TypeError(about._errors.cstring("ERROR: invalid input.")) raise TypeError(about._errors.cstring("ERROR: invalid input."))
self.domain = domain self.domain = domain
self.val = self.domain.enforce_values(val,extend=True) self.target = self.domain
self.val = self.domain.cast(val)
self.sym = True self.sym = True
self.uni = False self.uni = False
self.imp = False self.imp = False
self.target = self.domain
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_val(self,newval): def set_val(self, newval):
""" """
Sets the field values of the operator Sets the field values of the operator
...@@ -2753,17 +2782,18 @@ class vecvec_operator(operator): ...@@ -2753,17 +2782,18 @@ class vecvec_operator(operator):
------- -------
None None
""" """
self.val = self.domain.enforce_values(newval,extend=True) self.val = self.domain.cast(newval)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _multiply(self,x,**kwargs): ## > applies the operator to a given field def _multiply(self, x, **kwargs): ## > applies the operator to a given field
x_ = field(self.target,val=None,target=x.target) y = x.copy_empty(domain = self.target)
x_.val = self.val*self.domain.calc_dot(self.val,x.val) ## bypasses self.domain.enforce_values y.set_val(new_val = self.val * self.domain.calc_dot(self.val, x.val))
return x_ return y
def _inverse_multiply(self,x,**kwargs): def _inverse_multiply(self,x,**kwargs):
raise AttributeError(about._errors.cstring("ERROR: singular operator.")) raise AttributeError(about._errors.cstring(
"ERROR: singular operator."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -2806,19 +2836,22 @@ class vecvec_operator(operator): ...@@ -2806,19 +2836,22 @@ class vecvec_operator(operator):
of probing case. of probing case.
""" """
if(domain is None)or(domain==self.domain): if domain is None or domain == self.domain:
if(not self.domain.discrete): if self.domain.discrete == False:
return self.domain.calc_dot(self.val,self.domain.calc_weight(self.val,power=1)) return self.domain.calc_dot(self.val,
self.domain.calc_weight(self.val,power=1))
else: else:
return self.domain.calc_dot(self.val,self.val) return self.domain.calc_dot(self.val, self.val)
else: else:
return super(vecvec_operator,self).tr(domain=domain,**kwargs) ## probing ## probing
return super(vecvec_operator,self).tr(domain=domain,**kwargs)
def inverse_tr(self): def inverse_tr(self):
""" """
Inverse is ill-defined for this operator. Inverse is ill-defined for this operator.
""" """
raise AttributeError(about._errors.cstring("ERROR: singular operator.")) raise AttributeError(about._errors.cstring(
"ERROR: singular operator."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -2883,15 +2916,18 @@ class vecvec_operator(operator): ...@@ -2883,15 +2916,18 @@ class vecvec_operator(operator):
entries; e.g., as variance in case of an covariance operator. entries; e.g., as variance in case of an covariance operator.
""" """
if(domain is None)or(domain==self.domain): if (domain is None) or (domain==self.domain):
diag = np.real(self.val*np.conjugate(self.val)) ## bare diagonal diag = self.val * self.val.conjugate() ## bare diagonal
## weight if ... ## weight if ...
if(not self.domain.discrete)and(not bare): if self.domain.discrete == False and bare == False:
return self.domain.calc_weight(diag,power=1) return self.domain.calc_weight(diag, power=1)
else: else:
return diag return diag
else: else:
return super(vecvec_operator,self).diag(bare=bare,domain=domain,**kwargs) ## probing ## probing
return super(vecvec_operator,self).diag(bare=bare,
domain=domain,
**kwargs)
def inverse_diag(self): def inverse_diag(self):
""" """
......
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