python - 如何使用我在 python 中使用 swig 生成的包装函数？

Question

好吧，我是新来的。我终于使用 swig 和 numpy.i 成功包装了 python 程序中最昂贵的部分。该程序是二维波偏微分方程的有限差分格式。我的问题是我现在如何使用它？在 IPython 中导入后我可以看到它。

In [1]: import wave2

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

但是，当我使用它时，我收到一条错误消息:

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

如何将我的 2D numpy 数组转换为某种形式，使我能够使用它。还有另一个非常相似的 stackoverflow 对我没有帮助，尽管我一路上在这件事上找到了很多帮助。

这是标题:

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

这是c代码:

#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);

    }
}

这是我的界面:

   %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)}

这些是我用来编译和链接的命令:

$ 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

没有任何错误或警告。互联网上缺乏这样的中间示例，相信我，我已经搜索过了!，所以如果我们能让这个工作正常运行，它可以作为某人的一个很好的教程。请不要记下我的问题然后消失在夜色中。如果您认为我的某些编码需要改进，请告诉我，我现在基本上正在尝试自学一切......非常感谢您的帮助

哦，这里还有一个我正在尝试使用它的脚本。我还尝试在 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

另外，是的，我检查过以确保数组都是 numpy nd 数组类型

Answer 1

好吧，我能够在 Cython 中完成我想要的事情，感谢上帝。在这个过程中我发现 Cython 是一个更强大的工具!

因此结果是使用有限差分求解的二维波偏微分方程。主要计算已导出到 Cythonized 函数。该函数接收三个 2D np.ndarray 并返回一个。使用 numpy 和其他数据类型也更容易。

这是 Cythonized 函数:

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

这是我的 python 脚本，它调用它并返回它然后绘制结果。欢迎任何有关编写更好代码的建议...

'''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()