cupy integrate

genax.py

-import numpy, matplotlib.pyplot as plt
+import cupy as cp, matplotlib.pyplot as plt
 from numba import cuda
 from utils import expSpace
 
@@ -6,41 +6,57 @@ def gen_angels_to_pick(n, plot = False):
     if n % 2 == 0:
         raise "n%2==0"
     N = int((n + 3) / 2)
-    X1 = expSpace(0.0, numpy.pi/2.0, N)
-    X3 = expSpace(numpy.pi/2.0, numpy.pi, N)
-    X3 = -1*X3 + 3.0*numpy.pi/2.0
+    X1 = expSpace(0.0, cp.pi/2.0, N)
+    X3 = expSpace(cp.pi/2.0, cp.pi, N)
+    X3 = -1*X3 + 3.0*cp.pi/2.0
 
-    anglestopick = numpy.concatenate((X1, X3), axis=None)
-    anglestopick = numpy.unique(anglestopick) # Vigyázni vele!
+    anglestopick = cp.concatenate((X1, X3), axis=None)
+    anglestopick = cp.unique(anglestopick) # Vigyázni vele!
     anglestopick = anglestopick[1:-1]
 
     if plot:
-        Y = numpy.zeros(n)
+        Y = cp.zeros(n)
         plt.plot(anglestopick, Y, '|')
         plt.show()
     
     return anglestopick
 
-def angles_alap(anglestopick, plot = False):
-    alpha_arr = numpy.array([])
-    beta_arr = numpy.array([])
+parosit = cp.RawKernel(r'''
+extern "C" 
+__global__ 
+void my_add(const double* x1, const double* x2, double* a, double* b, const size_t m, const float PI) {
+    int tid = blockDim.x * blockIdx.x + threadIdx.x;
+    float alpha = x1[tid];
+    for (size_t i = 0; i < m; i++) {
+        float betha = x2[i];
+        if ((alpha + betha) < PI && betha > alpha && alpha > 0.0) {
+            a[tid*m+i] = alpha;
+            b[tid*m+i] = betha;
+        } else {
+            a[tid*m+i] = -1.0;
+            b[tid*m+i] = -1.0;
+        }
+    }
+}
+''', 'parosit')
 
-    for alpha in anglestopick:
-        for beta in anglestopick:
-            if alpha + beta < numpy.pi and beta >= alpha and alpha != 0.0: # vigyázni
-                alpha_arr = numpy.insert(alpha_arr, 0, alpha)
-                beta_arr = numpy.insert(beta_arr, 0, beta)
+def angles_alap(anglestopick, plot = False):
+    m = anglestopick.size
+    alpha_arr = cp.zeros((m, m), dtype=cp.float64)
+    beta_arr = cp.zeros((m, m), dtype=cp.float64)
+    
+    parosit((1,), (m,), (anglestopick, anglestopick, alpha_arr, beta_arr, m, cp.pi))
 
-    tompa_beta_arr = beta_arr[beta_arr > numpy.pi]
-    tompa_beta_mpi_arr = tompa_beta_arr - numpy.pi/2
-    hegyes_beta_arr = beta_arr[beta_arr <= numpy.pi]
-    tompa_alpha_arr = alpha_arr[beta_arr > numpy.pi]
-    hegyes_alpha_arr = alpha_arr[alpha_arr <= numpy.pi]
+    tompa_beta_arr = beta_arr[beta_arr > cp.pi]
+    tompa_beta_mpi_arr = tompa_beta_arr - cp.pi/2
+    hegyes_beta_arr = beta_arr[beta_arr <= cp.pi]
+    tompa_alpha_arr = alpha_arr[beta_arr > cp.pi]
+    hegyes_alpha_arr = alpha_arr[alpha_arr <= cp.pi]
 
-    tghB = numpy.tan(hegyes_beta_arr)
-    tghA = numpy.tan(hegyes_alpha_arr)
-    tgtB_mpi = numpy.tan(tompa_beta_mpi_arr)
-    tgtA = numpy.tan(tompa_alpha_arr)
+    tghB = cp.tan(hegyes_beta_arr)
+    tghA = cp.tan(hegyes_alpha_arr)
+    tgtB_mpi = cp.tan(tompa_beta_mpi_arr)
+    tgtA = cp.tan(tompa_alpha_arr)
 
     # hegyes
     sztgh = tghB*tghA
@@ -54,32 +70,32 @@ def angles_alap(anglestopick, plot = False):
     tCy = tgtA / mtgt
     tCx = tCy / tgtA
 
-    Cy = numpy.concatenate((hCy, tCy), axis=None)
-    Cx = numpy.concatenate((hCx, tCx), axis=None)
+    Cy = cp.concatenate((hCy, tCy), axis=None)
+    Cx = cp.concatenate((hCx, tCx), axis=None)
 
     return Cx, Cy
 
 
 def angles_ratet(anglestopick, plot = False):
-    alpha_arr = numpy.array([])
-    beta_arr = numpy.array([])
+    alpha_arr = cp.array([])
+    beta_arr = cp.array([])
 
     for alpha in anglestopick:
         for beta in anglestopick:
-            if alpha + beta < numpy.pi and alpha != 0.0: # vigyázni
-                alpha_arr = numpy.insert(alpha_arr, 0, alpha)
-                beta_arr = numpy.insert(beta_arr, 0, beta)
+            if alpha + beta < cp.pi and alpha != 0.0: # vigyázni
+                alpha_arr = cp.put(alpha_arr, 0, alpha)
+                beta_arr = cp.put(beta_arr, 0, beta)
 
-    tompa_beta_arr = beta_arr[beta_arr > numpy.pi]
-    tompa_beta_mpi_arr = tompa_beta_arr - numpy.pi/2
-    hegyes_beta_arr = beta_arr[beta_arr <= numpy.pi]
-    tompa_alpha_arr = alpha_arr[beta_arr > numpy.pi]
-    hegyes_alpha_arr = alpha_arr[alpha_arr <= numpy.pi]
+    tompa_beta_arr = beta_arr[beta_arr > cp.pi]
+    tompa_beta_mpi_arr = tompa_beta_arr - cp.pi/2
+    hegyes_beta_arr = beta_arr[beta_arr <= cp.pi]
+    tompa_alpha_arr = alpha_arr[beta_arr > cp.pi]
+    hegyes_alpha_arr = alpha_arr[alpha_arr <= cp.pi]
 
-    tghB = numpy.tan(hegyes_beta_arr)
-    tghA = numpy.tan(hegyes_alpha_arr)
-    tgtB_mpi = numpy.tan(tompa_beta_mpi_arr)
-    tgtA = numpy.tan(tompa_alpha_arr)
+    tghB = cp.tan(hegyes_beta_arr)
+    tghA = cp.tan(hegyes_alpha_arr)
+    tgtB_mpi = cp.tan(tompa_beta_mpi_arr)
+    tgtA = cp.tan(tompa_alpha_arr)
 
     # hegyes
     sztgh = tghB*tghA
@@ -93,14 +109,14 @@ def angles_ratet(anglestopick, plot = False):
     tCy = tgtA / mtgt
     tCx = tCy / tgtA
 
-    Ey = numpy.concatenate((hCy, tCy), axis=None)
-    Ex = numpy.concatenate((hCx, tCx), axis=None)
+    Ey = cp.concatenate((hCy, tCy), axis=None)
+    Ex = cp.concatenate((hCx, tCx), axis=None)
 
-    cos = numpy.cos(anglestopick)
-    sin = numpy.sin(anglestopick)
+    cos = cp.cos(anglestopick)
+    sin = cp.sin(anglestopick)
 
-    Dx = numpy.outer(numpy.full(anglestopick.size, 1.0), Ex).flatten()
-    Dy = numpy.outer(cos, Ey).flatten()
-    Dz = numpy.outer(sin, Ey).flatten()
+    Dx = cp.outer(cp.full(anglestopick.size, 1.0), Ex).flatten()
+    Dy = cp.outer(cos, Ey).flatten()
+    Dz = cp.outer(sin, Ey).flatten()
 
     return Dx, Dy, Dz
\ No newline at end of file

utils.py

-import numpy
+import cupy as cp
 
 def expSpace(min, max, N, exponentialliness = 20.0):
-    LinVec = numpy.linspace(0, numpy.log10(exponentialliness+1),N)
+    LinVec = cp.linspace(0, cp.log10(exponentialliness+1),N)
     return (max-min)/exponentialliness * (10.0**LinVec - 1) + min
\ No newline at end of file