BUG: may fixes

gyral_structure/basis/evaluator.py

@@ -278,6 +278,8 @@ class FuncRequestEvaluator(RequestEvaluator):
             arr = res[0][0] * res[0][1]
             for weight, part_res in res[1:]:
                 arr += weight * part_res
+            if not inverse and arr.shape[0] > self.request.npos:
+                arr = arr[:self.request.npos, :]
             results[key] = arr
 
         if self.potential:
@@ -553,6 +555,7 @@ class MultEvaluator(object):
                 from pycuda.autoinit import context
                 context.synchronize()
 
+            from .radial import RadialGPUArrays
             for req, evaluators in zip(self.request_list, self.evaluators.T):
                 for (bf, slc), eval in zip(self.basis_list, evaluators):
                     eval.evaluate_results(

gyral_structure/basis/radial.py

@@ -148,11 +148,8 @@ class RadialBasis(BasisFunc):
         :param return_compressed: return the request indices in a compressed format
         :return: tuple with the request and centroid indices in compressed format
         """
-        if self._precomputed_grids is not None:
-            print('check grid', req.radius(), self._precomputed_grids[0])
-            print(req)
-        if self._precomputed_grids is not None and req.radius() <= self._precomputed_grids[0]:
-            if not hasattr(self, '_ref_list_of_lists'):
+        if self._precomputed_grids is not None and np.amax(req.radius()) <= self._precomputed_grids[0]:
+            if not hasattr(self, '_ref_list_of_lists') or req.npos > self._ref_list_of_lists.size:
                 self._ref_list_of_lists = np.zeros(req.npos, dtype='object')
                 empty_arr = np.zeros(0, dtype='i4')
                 for idx in range(req.npos):
@@ -617,13 +614,15 @@ __global__ void matrix_mult_nd_invert({dtype} *derparam, {dtype} *derfield, int
         RadialGPUArrays.register(self.basis)
         RadialGPUArrays.register(self.request)
 
-    def update_indices(self, only_forward=False):
+    def update_indices(self, request=None, only_forward=False):
         del self.forward_idx
+        if request is None:
+            request = self.request
         if only_forward:
-            forward_compressed, forward_idx = self.basis.within_range(self.request, return_compressed=True)
+            forward_compressed, forward_idx = self.basis.within_range(request, return_compressed=True)
         else:
-            idx_req, idx_centroids = self.basis.within_range(self.request)
-            forward_compressed = sp.append(0, sp.cumsum(sp.bincount(idx_req, minlength=self.request.npos)))
+            idx_req, idx_centroids = self.basis.within_range(request)
+            forward_compressed = sp.append(0, sp.cumsum(sp.bincount(idx_req, minlength=request.npos)))
             forward_idx = idx_centroids  # [sp.argsort(idx_req)]; idx_req is always sorted
         self.forward_idx = (cuda.to_gpu_correct(forward_compressed),
                             cuda.to_gpu_correct(forward_idx))
@@ -631,7 +630,7 @@ __global__ void matrix_mult_nd_invert({dtype} *derparam, {dtype} *derfield, int
             return
 
         del self.backward_idx
-        backward_compressed = sp.append(0, sp.cumsum(sp.bincount(idx_centroids, minlength=self.request.npos)))
+        backward_compressed = sp.append(0, sp.cumsum(sp.bincount(idx_centroids, minlength=request.npos)))
         backward_idx = idx_req[sp.argsort(idx_centroids)]
         self.backward_idx = (cuda.to_gpu_correct(backward_compressed),
                              cuda.to_gpu_correct(backward_idx))
@@ -661,6 +660,7 @@ __global__ void matrix_mult_nd_invert({dtype} *derparam, {dtype} *derfield, int
             self.forward(*args, block=(cuda.nthreads, 1, 1),
                          grid=(int(sp.ceil(32 * self.request.npos / cuda.nthreads)), 1, 1))
 
+
     def combine_results(self, inverse=False):
         """Retrieves the vector field at the requested positions as computed by evaluate_results
 
@@ -675,8 +675,13 @@ __global__ void matrix_mult_nd_invert({dtype} *derparam, {dtype} *derfield, int
         RadialGPUArrays.clean()
 
     def update_pos(self, new_request):
-        RadialGPUArrays.set(self.request, RadialGPUArrays.get(self.request)[:new_request.npos])
+        new_arr = RadialGPUArrays.get(self.request)
+        if new_request.npos < self.request.npos:
+            new_arr = new_arr[:new_request.npos]
+        RadialGPUArrays.arrays[new_request] = new_arr
+        self.hidden_request = new_request
+           
         idx = self.req_params_cuda_names.index('all_pos')
         self.global_cuda_params[self.request][idx][:new_request.positions.size] = new_request.positions.astype(cuda.dtype).flatten()
-        self.update_indices(only_forward=True)
+        self.update_indices(new_request, only_forward=True)
 

gyral_structure/test/test_computation.py

@@ -58,6 +58,7 @@ def check_computation(basis_func, npos=13, noder=False, nodiv=True):
             #print(map1.method, method)
             assert map1.method == method
         res1 = map1(parameters)
+        print('running param_evaluator')
         resalt = basis_func.param_evaluator(parameters, method=method, nsim=npos//2)(req)
         if hasattr(basis_func, '_precomputed_grid'):
             basis_func._precomputed_grid = None
@@ -97,17 +98,17 @@ def check_computation(basis_func, npos=13, noder=False, nodiv=True):
 @pytest.mark.slowtest
 def test_computation():
     sp.random.seed(12345)
-    check_computation(basis.Fourier(150, 3, 10))
+    check_computation(basis.RadialBasis(basis.buhmann(3),
+                                        sp.randn(15, 3), size=sp.rand(15) * 3 + 0.1))
+    check_computation(basis.RadialBasis(basis.wendland(2), sp.randn(14, 2), size=1.82))
     check_computation(basis.ChargeDistribution(sp.randn(11, 3)), nodiv=False)
+    check_computation(basis.Fourier(150, 3, 10))
     mb = basis.SumBase([basis.Fourier(150, 3, 10), basis.ChargeDistribution(sp.randn(91, 3))])
     check_computation(mb, nodiv=False)
     for ndim in (2, 3):
         for order in (0, 1):
             print('gradient of ndim %i and order %i' % (ndim, order))
             check_computation(basis.Gradient(order, ndim=ndim))
-    check_computation(basis.RadialBasis(basis.wendland(2), sp.randn(14, 2), size=1.82))
-    check_computation(basis.RadialBasis(basis.buhmann(3),
-                                        sp.randn(15, 3), size=sp.rand(15) * 3 + 0.1))
     mb = basis.SumBase([basis.Fourier(150, 3, 10), basis.RadialBasis(basis.wendland(3), sp.randn(17, 3),
                                                                      size=sp.rand(17))])
     check_computation(mb)