Merge branch 'master' into mmm3

nifty/basic_arithmetics.py

@@ -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))
 

nifty/minimization/vl_bfgs.py

@@ -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):

nifty/operators/fft_operator/transformations/lm_transformation_factory.py

 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

nifty/spaces/lm_space/lm_space.py

@@ -30,13 +30,6 @@ from d2o import arange
 gl = gdi.get('libsharp_wrapper_gl')
 hp = gdi.get('healpy')
 
-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):
     """
@@ -182,11 +175,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):
         if sigma is None:
             sigma = np.sqrt(2) * np.pi / (self.lmax + 1)

setup.py

@@ -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