Skip to content
Snippets Groups Projects
Commit 020eadac authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into mmm4

# Conflicts:
#	nifty/spaces/lm_space/lm_space.py
parents 09bba27a bb9d467c
Branches
Tags
2 merge requests!69Master,!67Martin's monster merge part 4/N: revamp dependencies, switch to pyHealpix
Pipeline #
......@@ -31,10 +31,8 @@ def _math_helper(x, function):
result_val = x.val.apply_scalar_function(function)
result = x.copy_empty(dtype=result_val.dtype)
result.val = result_val
elif isinstance(x, distributed_data_object):
result = x.apply_scalar_function(function, inplace=False)
else:
result = function(np.asarray(x))
......
......@@ -77,8 +77,6 @@ class InformationStore(object):
self._sy_store = {}
self._yy_store = {}
# self.dot_matrix = {}
@property
def history_length(self):
return min(self.k, self.max_history_length)
......@@ -193,48 +191,6 @@ class InformationStore(object):
self.last_x = x.copy()
self.last_gradient = gradient.copy()
#
# k = self.k
# m = self.actual_history_length
# big_m = self.history_length
#
# # compute dot products
# for i in xrange(k-1, k-m-1, -1):
# # new_s with s
# key = (big_m+m, big_m+1+i)
# self.dot_matrix[key] = new_s.dot(self.s[i])
#
# # new_s with y
# key = (big_m+m, i+1)
# self.dot_matrix[key] = new_s.dot(self.y[i])
#
# # new_y with s
# if i != k-1:
# key = (big_m+1+i, k)
# self.dot_matrix[key] = new_y.dot(self.s[i])
#
# # new_y with y
# # actually key = (i+1, k) but the convention is that the first
# # index is larger than the second one
# key = (k, i+1)
# self.dot_matrix[key] = new_y.dot(self.y[i])
#
# # gradient with s
# key = (big_m+1+i, 0)
# self.dot_matrix[key] = gradient.dot(self.s[i])
#
# # gradient with y
# key = (i+1, 0)
# self.dot_matrix[key] = gradient.dot(self.y[i])
#
# # gradient with gradient
# key = (0, 0)
# self.dot_matrix[key] = gradient.dot(gradient)
#
# self.last_x = x
# self.last_gradient = gradient
#
class LimitedList(object):
def __init__(self, history_length):
......
import numpy as np
def buildLm(inp, **kwargs):
if inp.dtype == np.dtype('float32'):
return _buildLm_f(inp, **kwargs)
else:
return _buildLm(inp, **kwargs)
def buildIdx(inp, **kwargs):
if inp.dtype == np.dtype('complex64'):
return _buildIdx_f(inp, **kwargs)
else:
return _buildIdx(inp, **kwargs)
def buildLm(nr, lmax):
new_dtype = np.result_type(nr.dtype, np.complex64)
def _buildIdx_f(nr, lmax):
size = (lmax+1)*(lmax+1)
size = (len(nr)-lmax-1)/2+lmax+1
res = np.zeros([size], dtype=new_dtype)
res[0:lmax+1] = nr[0:lmax+1]
res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
return res
final=np.zeros([size], dtype=np.float32)
final[0:lmax+1] = nr[0:lmax+1].real
final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
return final
def _buildIdx(nr, lmax):
size = (lmax+1)*(lmax+1)
def buildIdx(nr, lmax):
if nr.dtype == np.dtype('complex64'):
new_dtype = np.float32
elif nr.dtype == np.dtype('complex128'):
new_dtype = np.float64
elif nr.dtype == np.dtype('complex256'):
new_dtype = np.float128
else:
raise TypeError("dtype of nr not supported.")
final=np.zeros([size], dtype=np.float64)
size = (lmax+1)*(lmax+1)
final = np.zeros([size], dtype=new_dtype)
final[0:lmax+1] = nr[0:lmax+1].real
final[lmax+1::2] = np.sqrt(2)*nr[lmax+1:].real
final[lmax+2::2] = np.sqrt(2)*nr[lmax+1:].imag
return final
def _buildLm_f(nr, lmax):
size = (len(nr)-lmax-1)/2+lmax+1
res=np.zeros([size], dtype=np.complex64)
res[0:lmax+1] = nr[0:lmax+1]
res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
return res
def _buildLm(nr, lmax):
size = (len(nr)-lmax-1)/2+lmax+1
res=np.zeros([size], dtype=np.complex128)
res[0:lmax+1] = nr[0:lmax+1]
res[lmax+1:] = np.sqrt(0.5)*(nr[lmax+1::2] + 1j*nr[lmax+2::2])
return res
......@@ -24,13 +24,6 @@ from nifty.spaces.space import Space
from d2o import arange
def _distance_array_helper(index_array, lmax):
u=2*lmax+1
index_half=(index_array+np.minimum(lmax,index_array)+1)//2
m = (np.ceil((u - np.sqrt(u*u - 8 * (index_half - lmax))) / 2)).astype(int)
res = (index_half - m * (u - m) // 2).astype(np.float64)
return res
class LMSpace(Space):
"""
......@@ -149,11 +142,19 @@ class LMSpace(Space):
distribution_strategy=distribution_strategy)
dists = dists.apply_scalar_function(
lambda x: _distance_array_helper(x, self.lmax),
lambda x: self._distance_array_helper(x, self.lmax),
dtype=np.float)
return dists
@staticmethod
def _distance_array_helper(index_array, lmax):
u = 2*lmax + 1
index_half = (index_array+np.minimum(lmax, index_array)+1)//2
m = (np.ceil((u - np.sqrt(u*u - 8*(index_half - lmax)))/2)).astype(int)
res = (index_half - m*(u - m)//2).astype(np.float64)
return res
def get_fft_smoothing_kernel_function(self, sigma):
return lambda x: np.exp(-0.5 * x * (x + 1) * sigma**2)
......
......@@ -17,7 +17,6 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from setuptools import setup, find_packages
import os
from Cython.Build import cythonize
import numpy
......
......@@ -55,16 +55,14 @@ def _distance_array_helper(index_arr, lmax):
index_half = index_arr
else:
if (index_arr - lmax) % 2 == 0:
index_half = (index_arr + lmax) / 2
index_half = (index_arr + lmax)/2
else:
index_half = (index_arr + lmax + 1) / 2
index_half = (index_arr + lmax + 1)/2
m = (
np.ceil(((2 * lmax + 1) - np.sqrt((2 * lmax + 1)**2 -
8 * (index_half - lmax))) / 2)
).astype(int)
m = np.ceil(((2*lmax + 1) - np.sqrt((2*lmax + 1)**2 -
8*(index_half - lmax)))/2).astype(int)
return index_half - m * (2 * lmax + 1 - m) // 2
return index_half - m*(2*lmax + 1 - m)//2
def get_distance_array_configs():
da_0 = [_distance_array_helper(idx, 5) for idx in np.arange(36)]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment