Commit 6e1782e5 authored by theos's avatar theos
Browse files

Corrected bugs in the distance_array functions of lm/gl/hp space and in...

Corrected bugs in the distance_array functions of lm/gl/hp space and in HPLMTransformation and LMHPTransformation
parent 1c4d7918
...@@ -44,9 +44,8 @@ class HPLMTransformation(SlicingTransformation): ...@@ -44,9 +44,8 @@ class HPLMTransformation(SlicingTransformation):
"ERROR: domain needs to be a HPSpace")) "ERROR: domain needs to be a HPSpace"))
lmax = 3 * domain.nside - 1 lmax = 3 * domain.nside - 1
mmax = lmax
result = LMSpace(lmax=lmax, mmax=mmax, dtype=np.dtype('float64')) result = LMSpace(lmax=lmax, dtype=np.dtype('float64'))
cls.check_codomain(domain, result) cls.check_codomain(domain, result)
return result return result
...@@ -62,21 +61,16 @@ class HPLMTransformation(SlicingTransformation): ...@@ -62,21 +61,16 @@ class HPLMTransformation(SlicingTransformation):
nside = domain.nside nside = domain.nside
lmax = codomain.lmax lmax = codomain.lmax
mmax = codomain.mmax
if 3 * nside - 1 != lmax: if 3 * nside - 1 != lmax:
raise ValueError(about._errors.cstring( raise ValueError(about._errors.cstring(
'ERROR: codomain has 3*nside-1 != lmax.')) 'ERROR: codomain has 3*nside-1 != lmax.'))
if lmax != mmax:
raise ValueError(about._errors.cstring(
'ERROR: codomain has lmax != mmax.'))
return None return None
def _transformation_of_slice(self, inp, **kwargs): def _transformation_of_slice(self, inp, **kwargs):
lmax = self.codomain.lmax lmax = self.codomain.lmax
mmax = self.codomain.mmax mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating): if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [hp.map2alm(x.astype(np.float64, [resultReal, resultImag] = [hp.map2alm(x.astype(np.float64,
......
...@@ -64,11 +64,6 @@ class LMHPTransformation(SlicingTransformation): ...@@ -64,11 +64,6 @@ class LMHPTransformation(SlicingTransformation):
nside = codomain.nside nside = codomain.nside
lmax = domain.lmax lmax = domain.lmax
mmax = domain.mmax
if lmax != mmax:
raise ValueError(about._errors.cstring(
'ERROR: domain has lmax != mmax.'))
if 3*nside - 1 != lmax: if 3*nside - 1 != lmax:
raise ValueError(about._errors.cstring( raise ValueError(about._errors.cstring(
...@@ -79,7 +74,7 @@ class LMHPTransformation(SlicingTransformation): ...@@ -79,7 +74,7 @@ class LMHPTransformation(SlicingTransformation):
def _transformation_of_slice(self, inp, **kwargs): def _transformation_of_slice(self, inp, **kwargs):
nside = self.codomain.nside nside = self.codomain.nside
lmax = self.domain.lmax lmax = self.domain.lmax
mmax = self.domain.mmax mmax = lmax
if issubclass(inp.dtype.type, np.complexfloating): if issubclass(inp.dtype.type, np.complexfloating):
[resultReal, resultImag] = [ltf.buildLm(x, lmax=lmax) [resultReal, resultImag] = [ltf.buildLm(x, lmax=lmax)
......
...@@ -160,22 +160,24 @@ class GLSpace(Space): ...@@ -160,22 +160,24 @@ class GLSpace(Space):
) )
dists = dists.apply_scalar_function( dists = dists.apply_scalar_function(
lambda x: self._distance_array_helper(divmod(int(x), self.nlon)), lambda x: self._distance_array_helper(divmod(x, self.nlon)),
dtype=np.float dtype=np.float
) )
return dists return dists
def _distance_array_helper(self, qr_tuple): def _distance_array_helper(self, qr_tuple):
numerator = np.sqrt(np.sin(qr_tuple[1])**2 + lat = qr_tuple[0]*(np.pi/self.nlat)
(np.sin(qr_tuple[0]) * np.cos(qr_tuple[1]))**2) lon = qr_tuple[1]*(2*np.pi/self.nlon)
denominator = np.cos(qr_tuple[0]) * np.cos(qr_tuple[1]) numerator = np.sqrt(np.sin(lat)**2 +
(np.sin(lon) * np.cos(lat))**2)
denominator = np.cos(lon) * np.cos(lat)
return np.arctan(numerator / denominator) return np.arctan(numerator / denominator)
def get_smoothing_kernel_function(self, sigma): def get_smoothing_kernel_function(self, sigma):
if sigma is None: if sigma is None:
sigma = np.sqrt(2) * np.pi / self.nlat sigma = np.sqrt(2) * np.pi
return lambda x: np.exp((-0.5 * x**2) / sigma**2) return lambda x: np.exp((-0.5 * x**2) / sigma**2)
......
...@@ -176,20 +176,18 @@ class HPSpace(Space): ...@@ -176,20 +176,18 @@ class HPSpace(Space):
distribution_strategy=distribution_strategy distribution_strategy=distribution_strategy
) )
# setting the center to fixed value # translate distances to 3D unit vectors on a sphere,
center_vec = (1, 0, 0) # extract the first entry (simulates the scalar product with (1,0,0))
# and apply arccos
dists = dists.apply_scalar_function( dists = dists.apply_scalar_function(
lambda x: np.arccos(np.dot(hp.pix2vec(self.nside, int(x)), lambda z: np.arccos(hp.pix2vec(self.nside, z)[0]),
center_vec)), dtype=np.float)
dtype=np.float
)
return dists return dists
def get_smoothing_kernel_function(self, sigma): def get_smoothing_kernel_function(self, sigma):
if sigma is None: if sigma is None:
sigma = np.sqrt(2) * np.pi / self.nlat sigma = np.sqrt(2) * np.pi
return lambda x: np.exp((-0.5 * x**2) / sigma**2) return lambda x: np.exp((-0.5 * x**2) / sigma**2)
......
...@@ -113,15 +113,12 @@ class LMSpace(Space): ...@@ -113,15 +113,12 @@ class LMSpace(Space):
self._lmax = self._parse_lmax(lmax) self._lmax = self._parse_lmax(lmax)
def distance_array(self, distribution_strategy): def distance_array(self, distribution_strategy):
dists = arange( dists = arange(start=0, stop=self.shape[0],
start=0, stop=self.shape[0], distribution_strategy=distribution_strategy)
distribution_strategy=distribution_strategy
)
dists = dists.apply_scalar_function( dists = dists.apply_scalar_function(
lambda x: _distance_array_helper(x, self.lmax), lambda x: _distance_array_helper(x, self.lmax),
dtype=np.float dtype=np.float)
)
return dists return dists
......
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