memcof and evlmem implementation Help Needed


McBuell
11-17-2008, 02:14 PM
I'm interested in developping a Python wrapper for the memcof and evlmem c++ functions.

After some minor troubles in getting started (too bad there seems to be no examples .cpp files anymore in the NR3 release), I managed (with hints of the old NR2 xmemcof.cpp and evlmem.cpp) to write some code for the NR3 functions.

Here is what I've got so far: pymesa.cpp


#include "nr3.h"
#include "roots_poly.h"
#include "linpredict.h"
using namespace std;


int mesa(VecDoub_I &data,int N, int M, int NFDT, VecDoub_O &return_x, VecDoub_O &return_y)
{
int i;
Doub fdt, pm, temp_evlmem;
VecDoub cof(M);

memcof(data,pm,cof);

//cout << "Coefficients for spectral estimation of spctrl.dat" << endl;
//cout << fixed << setprecision(6);
//for (i=0;i<M;i++) {
// cout << "a[" << setw(2) << i << "] = ";
// cout << setw(12) << cof[i] << endl;
//}
//cout << endl << "a0 = " << setw(12) << pm << endl;

//cout << "Power spectum estimate of data in spctrl.dat" << endl;
//cout << " f*delta power" << endl;
//cout << fixed << setprecision(6);

i=0;
for (fdt=0.0;fdt<=0.5;fdt+=0.5/NFDT) {
//cout << setw(12) << fdt;
return_x[i] = fdt;
temp_evlmem = evlmem(fdt,cof,pm);
//cout << setw(13) << temp_evlmem << endl;
return_y[i] = temp_evlmem;
i += 1;
}
return 0;

}

int main(void)
{
cout << "starting from prompt... take some initial parameters...." << endl;
const int N=1000,M=10,NFDT=16;
int result;
VecDoub data(N);
VecDoub return_x(NFDT+1);
VecDoub return_y(NFDT+1);
int i;
ifstream fp("spctrl.dat");

if (fp.fail())
throw("Data file spctrl.dat not found");
for (i=0;i<N;i++) {
fp >> data[i];
//cout << "data[" << i <<"]: "<<data[i]<<endl;
}
fp.close();

mesa(data,N,M,NFDT, return_x, return_y);

for (i=0;i<=NFDT;i++) {
cout << return_x[i] << ": " << return_y[i] << endl;
}

}



I can compile pymesa.cpp without any errors.
The resulting output is executable and seems to work fine.

Now I've compiled it with the -shared parameter and wrapped it in Python this way:


import ctypes, numpy
ctypes.cdll.LoadLibrary("_pymesa.dll")
pymesa = ctypes.CDLL("_pymesa.dll")

Ns=2048

x=numpy.arange(0.0,1.0*Ns,1.0)
A=numpy.sin(2.0*numpy.pi*x/150.0)


return_x = ctypes.POINTER(ctypes.c_double) * NFDT
return_x1 = return_x()
return_y = ctypes.POINTER(ctypes.c_double) * NFDT
return_y1 = return_x()

#Create input array
input_array = ctypes.POINTER(ctypes.c_double) * N
input_array1 = input_array()

#fill input array with the values of the sin array A
i = 0
while i < N:
B = ctypes.c_double(A[i])
input_array1[i] = ctypes.pointer(B)
i += 1

# execute the dll mesa function
pymesa.mesa(input_array1, ctypes.c_int(1000),ctypes.c_int(10),ctypes.c_int(1 6), ctypes.pointer(return_x1), ctypes.pointer(return_y1))




I'm just not able to get the right input_array data structure in ctypes for the VecDoub, the mesa call ends with memory read violation.

The pymesa.dll is working within python, because I can call pymesa.main() with proper results.

Could anybody help me understanding the VecDoub structure ?
Or maybe a python/ctypes guru is able to declare the propper ctypes and pointer property (and how to construct them).

Any Help is realy appriciated :D

Greetings from Switzerland

McBuell
11-17-2008, 04:30 PM
hiho again...

problem is solved !

Now I pass a normal c_double * N array to the c DLL function, and declare the VecDoub Instance within the c function, rather then trying to mimic VecDoub in Python ctypes...

work great !

Cheers