SWIG python wrapper (should install)

Author: johnh2o2

.gitignore

@@ -1,2 +1,3 @@
 *o
 culsp
+*swp

Makefile

@@ -26,10 +26,11 @@ CC=gcc
 cuda_lib=/usr/local/cuda/lib64
 cuda_inc=/usr/local/cuda/include
 
+python_inc=/usr/include/python2.7 
 BLOCK_SIZE=256
 
-NVCCFLAGS := --ptxas-options=-v -DBLOCK_SIZE=$(BLOCK_SIZE) -arch $(ARCH)
-CXXFLAGS := -DBLOCK_SIZE=$(BLOCK_SIZE)
+NVCCFLAGS := -Xcompiler -fpic --ptxas-options=-v -DBLOCK_SIZE=$(BLOCK_SIZE) -arch $(ARCH)
+CXXFLAGS := -fPIC -DBLOCK_SIZE=$(BLOCK_SIZE)
 LINK := -largtable2 -lm -lcudart -L$(cuda_lib)
 
 all : $(EXECUTABLE)
@@ -41,5 +42,23 @@ $(EXECUTABLE): culsp.o periodogram.o
 periodogram.o : periodogram.cpp
 	$(CXX) $(CXXFLAGS) -c -o $@ $^
 
+periodogram_nomain.o : periodogram.cpp
+	$(CXX) -Dmain=oldmain $(CXXFLAGS) -c -o $@ $^
+
 culsp.o : culsp.cu
-	$(NVCC) $(NVCCFLAGS) -c -o $@ $^ -I$(cuda_inc)
+	$(NVCC) -Xcompiler -fpic $(NVCCFLAGS) -c -o $@ $^ -I$(cuda_inc)
+
+culsp_wrap.o : culsp_wrap.cpp
+	$(CXX) -fPIC $(CXXFLAGS) -c -o $@ $^ -I$(python_inc)
+
+python : culsp_wrap.o culsp.o periodogram_nomain.o
+	$(CXX) -fPIC $(CXXFLAGS) -shared -o _culspy.so $^ $(LINK) -lpython2.7 
+	#mkdir culspy
+	#mv culspy.py culspy/
+	#mv _culspy.so culspy/
+	#touch culspy/__init__.py
+
+clean :
+	rm -f *o *so *pyc $(EXECUTABLE)
+	rm -f -r culspy/
+	rm -r -f build/

culsp.cu

@@ -20,7 +20,7 @@
 #include <stdlib.h>
 #include <string.h>
 #include <argtable2.h>
-
+#include "culsp.h"
 #include "periodogram.h"
 #include "culsp_kernel.cu"
 
@@ -36,9 +36,9 @@
 
 // Forward declarations
 
-void initialize (int, char **, char **, char **, float *, float *, int *);
-void initialize_cuda (int, int);
-void eval_LS_periodogram (int, int, float, float *, float *, float *);
+//void initialize (int, char **, char **, char **, float *, float *, int *);
+//void initialize_cuda (int, int);
+//void eval_LS_periodogram (int, int, float, float *, float *, float *);
 
 // Main program
 
@@ -60,6 +60,7 @@ main( int argc, char** argv)
 
   float df;
   float *P;
+  float minf = 0.0;
 
   // Initialize
 
@@ -79,15 +80,16 @@ main( int argc, char** argv)
 
   // Initialize CUDA
 
-  initialize_cuda(device, N_f);
+  initialize_cuda(device);
 
   // Start the timer
 
   double time_a = get_time();
 
   // Evaluate the Lomb-Scargle periodogram
- 
-  eval_LS_periodogram(N_t, N_f, df, t, X, P);
+  // set minf to 0 here (simplicity; I'm interacting with this
+  // through python anyway.
+  eval_LS_periodogram(N_t, N_f, df, minf,  t, X, P);
 
   // Stop the timer
 
@@ -166,9 +168,9 @@ initialize (int argc, char **argv, char **filename_in, char **filename_out,
 ////
 
 void
-initialize_cuda (int device, int N_f)
-{
-
+initialize_cuda (int device)
+{ 
+  
   // Select the device
 
   CUDA_CALL(cudaSetDevice(device));
@@ -187,7 +189,7 @@ initialize_cuda (int device, int N_f)
 ////
 
 void
-eval_LS_periodogram (int N_t, int N_f, float df, 
+eval_LS_periodogram (int N_t, int N_f, float df, float minf, 
 		     float *t, float *X, float *P)
 {
 
@@ -209,13 +211,13 @@ eval_LS_periodogram (int N_t, int N_f, float df,
   dim3 grid_dim(N_f/BLOCK_SIZE, 1, 1);
   dim3 block_dim(BLOCK_SIZE, 1, 1);
 
-  printf("Grid of %d frequency blocks of size %d threads\n", N_f/BLOCK_SIZE, BLOCK_SIZE);
+  //printf("Grid of %d frequency blocks of size %d threads\n", N_f/BLOCK_SIZE, BLOCK_SIZE);
 
   // Launch the kernel
 
-  printf("Launching kernel...\n");
+  //printf("Launching kernel...\n");
 
-  culsp_kernel<<<grid_dim, block_dim>>>(d_t, d_X, d_P, df, N_t);
+  culsp_kernel<<<grid_dim, block_dim>>>(d_t, d_X, d_P, df, N_t, minf);
 
   cudaError_t err = cudaGetLastError();
   if(err != cudaSuccess) {
@@ -226,7 +228,7 @@ eval_LS_periodogram (int N_t, int N_f, float df,
 
   CUDA_CALL(cudaThreadSynchronize());
 
-  printf("Completed!\n");
+  //printf("Completed!\n");
 
   // Copy data from the device
 

culsp.h

@@ -0,0 +1,6 @@
+#ifndef CULSP_H
+#define CULSP_H
+void initialize(int, char **, char **, char **, float *, float *, int *);
+void initialize_cuda(int);
+void eval_LS_periodogram(int, int, float, float, float *, float *, float *);
+#endif

culsp.i

@@ -0,0 +1,83 @@
+%module culspy
+%include typemaps.i
+
+%{
+extern void initialize_cuda(int);
+extern void set_frequency_params(int, float *, float, float, int *, float *);
+extern void eval_LS_periodogram(int, int, float,float, float *, float *, float *);
+%}
+
+extern void initialize_cuda(int);
+
+%apply int *OUTPUT { int *OUTPUT1 };
+extern void set_frequency_params(int, float *, float, float, int *OUTPUT1, float *OUTPUT);
+extern void eval_LS_periodogram(int, int, float,float,  float *, float *, float *);
+
+%inline %{
+
+float *get_float_array(int n){
+        float *x = (float *)malloc(n * sizeof(float));
+        return x;
+}
+float get_val(float *x,int i){
+        return x[i];
+}
+void set_val(float *x, int i, float val){
+        x[i] = val;
+}
+
+int get_block_size(){
+    return BLOCK_SIZE;
+}
+%}
+
+%pythoncode %{
+BLOCK_SIZE = _culspy.get_block_size()
+
+def _convert_to_c(arr):
+    N = len(arr);
+    carr = _culspy.get_float_array(N);
+    for i, val in enumerate(arr):
+        _culspy.set_val(carr, i, val)
+    return carr
+
+def _convert_to_py(carr, N):
+    return [ _culspy.get_val(carr, i) for i in range(N) ]
+
+def _insettings(arr, settings):
+    return [ v in settings for v in arr ]
+
+def correct_nf(Nf):
+    if Nf % BLOCK_SIZE != 0:
+        return (BLOCK_SIZE - Nf % BLOCK_SIZE) % BLOCK_SIZE
+    return Nf
+
+
+def LSP(t, x, f_over=1.0, f_high=1.0, minf=0.0, maxf=None, Nf=None ):
+    Nt = len(t)
+    ct = _convert_to_c(t)
+    cx = _convert_to_c(x)
+
+    if maxf is not None and Nf is not None:
+        Nf = correct_nf(Nf)
+        df = (maxf - minf)/Nf
+    else:
+        ct = _convert_to_c(t)
+        Nf0, df0 = _culspy.set_frequency_params(Nt, ct, f_over, f_high)
+        if not Nf is None:
+            Nf = correct_nf(Nf)
+            df = (Nf0 * df0)/Nf
+        else:
+            Nf = Nf0
+            df = df0
+    
+    cpower = _culspy.get_float_array(Nf)
+    _culspy.eval_LS_periodogram(Nt, Nf, df, minf, ct, cx, cpower)
+
+    freqs = [ minf + df * i for i in range(Nf) ]
+    power = _convert_to_py(cpower, Nf)
+    
+    
+    return freqs, power
+
+%}

culsp_kernel.cu

@@ -24,15 +24,15 @@
 
 __global__ void
 __launch_bounds__(BLOCK_SIZE)
-culsp_kernel(float *d_t, float *d_X, float *d_P, float df, int N_t)
+culsp_kernel(float *d_t, float *d_X, float *d_P, float df, int N_t, float minf)
 {
 
   __shared__ float s_t[BLOCK_SIZE];
   __shared__ float s_X[BLOCK_SIZE];
 
   // Calculate the frequency
 
-  float f = (blockIdx.x*BLOCK_SIZE+threadIdx.x+1)*df;
+  float f = (blockIdx.x*BLOCK_SIZE+threadIdx.x+1)*df + minf;
 
   // Calculate the various sums