How to use a wrapper function that i generated using swig in python?
Question: Okay, I am new to swig. I have finally successfully wrapped the most expensive
part of my python program using swig and also numpy.i . The program is a
finite difference scheme for the 2D wave PDE . My question is how do I use it
now ? I can see it after I import it within IPython.

    
    
    In [1]: import wave2
    
    In [2]: wave2.wave_prop
    Out[2]: <function _wave2.wave_prop>
    

But, when I go to use it I get an error saying:

    
    
    TypeError: in method 'wave_prop', argument 1 of type 'float **'
    

How can I transform my 2D numpy arrays to some form that will enable me to use
this. There is another stackoverflow that is very similar that did not help
me, although I have found a lot of help with this matter along the way.

Here is the header:

    
    
    void wave_prop(float** u_prev ,int Lx, int Ly,float** u ,int Lx2, int Ly2,float** u_next,int Lx3,int Ly3  );
    

Here is the c code:

    
    
    #include<string.h>
    #include<stdlib.h>
    #include<stdio.h>
    #include<math.h>
    
    #define n 100
    
    
    
    void wave_prop(float** u_prev ,int Lx,int Ly,float** u ,int Lx2,int Ly2,float** u_next,int Lx3,int Ly3 ){
    
    int dx=1;
    int dy=1;
    float c=1;
    float dt =1;
    int t_old=0;int t=0;int t_end=150;
    int x[Lx];
    int y[Ly];
    
    for(int i=0;i<=99;i++){
            x[i]=i;
            y[i]=i;
        }
    
    
    while(t<t_end){
        t_old=t; t +=dt;
        //the wave steps through time
        for (int i=1;i<99;i++){
            for (int j=1;j<99;j++){
                    u_next[i][j] = - u_prev[i][j] + 2*u[i][j] + \
                            (c*dt/dx)*(c*dt/dx)*u[i-1][j] - 2*u[i][j] + u[i+1][j] + \
                    (c*dt/dx)*(c*dt/dx)*u[i][j-1] - 2*u[i][j] + u[i][j+1];
                    }
                 }
    
        //set boundary conditions to 0
    
        for (int j=0;j<=99;j++){ u_next[0][j] = 0;}
        for (int i=0;i<=99;i++){ u_next[i][0] = 0;}
        for (int j=0;j<=99;j++){ u_next[Lx-1][j] = 0;}
        for (int i=0;i<=99;i++){ u_next[i][Ly-1] = 0;}
    
        //memcpy(dest, src, sizeof (mytype) * rows * coloumns);
        memcpy(u_prev, u, sizeof (float) * Lx * Ly);
        memcpy(u, u_next, sizeof (float) * Lx * Ly);
    
        }
    }
    

And here is my interface:

    
    
       %module wave2
    %{
        #define SWIG_FILE_WITH_INIT
        #include "wave2.h"
    
    %}
    
    %include "numpy.i"
    
    %init %{
        import_array();
    %}
    
    %include "wave2.h"
    %apply (float** INPLACE_ARRAY2, int DIM1, int DIM2) { (float** u_prev,int Lx,int Ly ),(float** u,int Lx2,int Ly2),(float* u_next,int Lx3,int Ly3)}
    

These are the commands I used to compile and link:

    
    
    $ swig -python wave2.i 
    $ gcc -c -fpic wave2.c wave2_wrap.c -I/usr/include/python2.7 -std=c99
    $ gcc -shared wave2.o wave2_wrap.o -o _wave2.so
    

Without any errors or warnings. There is a lack of intermediate examples like
this on the internet, trust me I have scoured!, so if WE can get this working
it could serve as a good tutorial for someone. Please do not mark my question
down then go away into the night . If you think some of my coding needs
improvement please let me know I am trying to basically teach myself
everything right now... Thank you very much for your help

Oh and also here is a script that I am trying to use it in. I have also tried
to use the function in other ways within IPython...

    
    
    '''George Lees Jr.
    2D Wave pde '''
    
    from numpy import *
    import numpy as np
    import matplotlib.pyplot as plt
    from wave2 import * 
    import wave2
    
    #declare variables
    #need 3 arrays u_prev is for previous time step due to d/dt
    
    Lx=Ly = (100)
    n=100
    dx=dy = 1
    x=y = np.array(xrange(Lx))
    u_prev = np.array(zeros((Lx,Ly),float))
    u = np.array(zeros((Lx,Ly),float))
    u_next = np.array(zeros((Lx,Ly),float))
    c = 1 #constant velocity
    dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2))
    t_old=0;t=0;t_end=150
    
    #set Initial Conditions and Boundary Points
    #I(x) is initial shape of the wave
    #f(x,t) is outside force that creates waves set =0
    
    def I(x,y): return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0)
    def f(x,t,y): return 0
    
    #set up initial wave shape
    
    for i in xrange(100):
        for j in xrange(100):
            u[i,j] = I(x[i],y[j])
    
    #copy initial wave shape for printing later
    
    u1=u.copy()
    
    #set up previous time step array
    
    for i in xrange(1,99):
        for j in xrange(1,99):
                u_prev[i,j] = u[i,j] + 0.5*((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
                0.5*((c*dt/dy)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \
                dt*dt*f(x[i], y[j], t)
    
    #set boundary conditions to 0
    
    for j in xrange(100): u_prev[0,j] = 0
    for i in xrange(100): u_prev[i,0] = 0
    for j in xrange(100): u_prev[Lx-1,j] = 0
    for i in xrange(100): u_prev[i,Ly-1] = 0
    
    wave2.wave_prop( u_prev ,Lx ,Ly , u , Lx, Ly, u_next,Lx,Ly )
    
    #while t<t_end:
    #   t_old=t; t +=dt
        #the wave steps through time
    #   for i in xrange(1,99):
    #       for j in xrange(1,99):
    #               u_next[i,j] = - u_prev[i,j] + 2*u[i,j] + \
    #                       ((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
    #               ((c*dt/dx)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) + \
    #                       dt*dt*f(x[i], y[j], t_old)
    #
    #   #set boundary conditions to 0
    #
    #   for j in xrange(100): u_next[0,j] = 0
    #   for i in xrange(100): u_next[i,0] = 0
    #   for j in xrange(100): u_next[Lx-1,j] = 0
    #   for i in xrange(100): u_next[i,Ly-1] = 0
    
        #set prev time step equal to current one
    #   u_prev = u.copy(); u = u_next.copy(); 
    
    fig = plt.figure()
    plt.imshow(u,cmap=plt.cm.ocean)
    plt.colorbar()
    plt.show()
    print u_next
    

Also yes I checked to make sure that the arrays were all numpy nd array types


Answer: Okay so i was able to accomplish what i wanted to in Cython, thank god. In the
process I found Cython much more powerful of a tool !

So the result is a 2D wave PDE solved using a finite difference. The main
computation has been exported to a Cythonized function. The function takes in
three 2D np.ndarrays and returns one back. Much easier to work with numpy and
other data types as well.

Here is the Cythonized function:

    
    
    from numpy import *
    cimport numpy as np
    
    def cwave_prop( np.ndarray[double,ndim=2] u_prev, np.ndarray[double,ndim=2] u, np.ndarray[double,ndim=2] u_next):
    
        cdef double t = 0
        cdef double t_old = 0
        cdef double t_end = 100
        cdef int i,j
        cdef double c = 1
        cdef double Lx = 100
        cdef double Ly = 100
        cdef double dx = 1
        cdef double dy = 1
        cdef double dt = (1/(c))*(1/(sqrt(1/dx**2 + 1/dy**2)))
    
        while t<t_end:
            t_old=t; t +=dt
    
            #wave steps through time and space
    
            for i in xrange(1,99):
                for j in xrange(1,99):
                    u_next[i,j] = - u_prev[i,j] + 2*u[i,j] + \
                            ((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
                    ((c*dt/dx)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1])
    
            #set boundary conditions of grid to 0
    
            for j in xrange(100): u_next[0,j] = 0
            for i in xrange(100): u_next[i,0] = 0
            for j in xrange(100): u_next[Lx-1,j] = 0
            for i in xrange(100): u_next[i,Ly-1] = 0
    
            #set prev time step equal to current one
            for i in xrange(100):
                for j in xrange(100):       
                    u_prev[i,j] = u[i,j]; 
                    u[i,j] = u_next[i,j]; 
    
    
    
        print u_next
    

And here is my python script that calls it and also returns it back then plots
the result. Any suggestions on writing better code are welcome...

    
    
    '''George Lees Jr.
    2D Wave pde '''
    
    from numpy import *
    import numpy as np
    import matplotlib.pyplot as plt 
    import cwave2
    np.set_printoptions(threshold=np.nan)
    
    #declare variables
    #need 3 arrays u_prev is for previous time step due to time derivative
    
    Lx=Ly = (100)                       #Length of x and y dims of grid
    dx=dy = 1                       #derivative of x and y respectively
    x=y = np.array(xrange(Lx))              #linspace to set the initial condition of wave
    u_prev=np.ndarray(shape=(Lx,Ly), dtype=np.double)   #u_prev 2D grid for previous time step needed bc of time derivative
    u=np.ndarray(shape=(Lx,Ly), dtype=np.double)        #u 2D grid
    u_next=np.ndarray(shape=(Lx,Ly), dtype=np.double)   #u_next for advancing the time step #also these are all numpy ndarrays
    c = 1                           #setting constant velocity of the wave
    dt = (1/float(c))*(1/sqrt(1/dx**2 + 1/dy**2))       #we have to set dt specifically to this or numerical approx will fail!
    print dt
    
    #set Initial Conditions and Boundary Points
    #I(x) is initial shape of the wave
    
    def I(x,y): return exp(-(x-Lx/2.0)**2/2.0 -(y-Ly/2.0)**2/2.0)
    
    #set up initial wave shape
    
    for i in xrange(100):
        for j in xrange(100):
            u[i,j] = I(x[i],y[j])
    
    
    #set up previous time step array
    
    for i in xrange(1,99):
        for j in xrange(1,99):
                u_prev[i,j] = u[i,j] + 0.5*((c*dt/dx)**2)*(u[i-1,j] - 2*u[i,j] + u[i+1,j]) + \
                0.5*((c*dt/dy)**2)*(u[i,j-1] - 2*u[i,j] + u[i,j+1]) 
    
    #set boundary conditions to 0
    
    for j in xrange(100): u_prev[0,j] = 0
    for i in xrange(100): u_prev[i,0] = 0
    for j in xrange(100): u_prev[Lx-1,j] = 0
    for i in xrange(100): u_prev[i,Ly-1] = 0
    
    #call C function from Python
    cwave2.cwave_prop( u_prev , u , u_next )
    #returned u (2D np.ndarray)
    
    from tempfile import TemporaryFile
    outfile = TemporaryFile()
    np.save(outfile,u)
    fig = plt.figure()
    plt.imshow(u,cmap=plt.cm.ocean)
    plt.colorbar()
    plt.show()
    

