C-extension in Python -- returning results
Keith Farmer
deoradh at yahoo.com
Wed Nov 7 06:55:31 EST 2001
Any help on the following would be appreciated. After all -- it's for
science.
I've got an extension that I'm writing -- a specialized math/database
package for a project of mine. I'm not the most well-practiced when
it comes to C, so I'm having some difficulty picking out the problem.
While I could use Numpy for this, I'd rather for efficiency's sake do
this directly, to avoid having to repackage the data to and from
Python several times.
The gist: the database represents a several-thousand term equation
giving planetary positions in time, returning the positions in
spherical polar coordinates. It's possible (I've done so in Numpy) to
sample this equation and fit smaller-period slices (~weeks, year at
most) to a set of 5th degree polynomials. This involves solving a set
of linear equations using standard techniques.
My error lies somewhere in the linear equation solver.
PythonWin 2.1 (#15, Apr 16 2001, 18:25:49) [MSC 32 bit (Intel)] on
win32.
Portions Copyright 1994-2001 Mark Hammond (MarkH at ActiveState.com) -
see 'Help/About PythonWin' for further copyright information.
>>> import AstroMath
>>> AstroMath.mkqvsop(2448976.5,2449076.5,10,3)
[-1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND,
-1.#IND,
-1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND, -1.#IND,
-1.#IND,
-1.#IND, -1.#IND]
mkqvsop samples the positions and performs the matrix math to yield
the polynomial fit. I've already verified that the function providing
the coodinates already works(it's used at
http://www.thuban.org/astro). I don't understand what the '-1.#IND'
returns indicate.
Following is the VC++ source code for the dll:
/* astromath.cpp -- astronomical calculations
*/
#import "c:\program files\common files\system\ado\msado15.dll" \
rename ( "EOF", "adoEOF" ) no_namespace
#include "Python.h"
#include <stdio.h>
#include <math.h>
#define SUCCESS 0
#define FAIL 1
#define PI 3.141592653589
#define RAD(x) (x*180/PI)
int linsys6x3(double [6][6], double [6][3]);
int calcVsop(double, int, double []);
int makeQVsop(double, double, int, int, double [3][6]);
static PyObject *AstroMathError;
static PyObject *AstroMath_version(PyObject *, PyObject *);
static PyObject *AstroMath_vsop(PyObject *, PyObject *);
static PyObject *AstroMath_mkqvsop(PyObject *, PyObject *);
// Initialize OLE
#define VERSION 2001.0904
#define TINY 1.0e-20
struct InitOle {
InitOle() { ::CoInitialize(NULL); }
~InitOle() { ::CoUninitialize(); }
} _init_InitOle_;
// Python init
static PyMethodDef AstroMathMethods[] = {
{ "vsop", AstroMath_vsop, METH_VARARGS },
{ "mkqvsop", AstroMath_mkqvsop, METH_VARARGS },
{ "version", AstroMath_version, METH_VARARGS },
{ NULL, NULL } /* Sentinel */
};
#ifdef MS_WIN32
__declspec(dllexport)
#endif
void initAstroMath()
{
PyObject *m, *d;
m = Py_InitModule("AstroMath", AstroMathMethods);
d = PyModule_GetDict(m);
AstroMathError = PyErr_NewException("AstroMath.error", NULL, NULL);
PyDict_SetItemString(d, "error", AstroMathError);
}
// Python stubs
static PyObject *AstroMath_version(PyObject *self, PyObject *args)
{
return Py_BuildValue("d", VERSION );
}
static PyObject *AstroMath_vsop(PyObject *self, PyObject *args)
{
double results[3] = {0, 0, 0};
int body;
double t0;
if (!PyArg_ParseTuple(args, "di", &t0, &body))
return NULL;
if (calcVsop(t0, body, results) == FAIL)
/* Error */
return AstroMathError;
else
/* Return value */
return Py_BuildValue("[ddd]", results[0], results[1], results[2]);
}
static PyObject *AstroMath_mkqvsop(PyObject *self, PyObject *args)
{
double results[3][6] = {{0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0}};
int N, body;
double startJD, endJD;
if (!PyArg_ParseTuple(args, "ddii", &startJD, &endJD, &N, &body))
return NULL;
if (makeQVsop(startJD, endJD, N, body, results) == FAIL)
/* Error */
return AstroMathError;
else
/* Return value */
return Py_BuildValue("[dddddddddddddddddd]",
results[0][0],
results[0][1],
results[0][2],
results[0][3],
results[0][4],
results[0][5],
results[1][0],
results[1][1],
results[1][2],
results[1][3],
results[1][4],
results[1][5],
results[2][0],
results[2][1],
results[2][2],
results[2][3],
results[2][4],
results[2][5]
);
}
// Main functions
int calcVsop(double t0, int body, double results[])
{
char sql[200];
double storage[3][6] = {{0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0}};
long cnt = 0;
long tot = 0; /* Define return variable */
double a, b, c = 0;
int coord;
int powerOfT;
t0 = (t0 - 2451545.0) / 365250;
sprintf(
sql,
"SELECT [a], [b], [c], [coord], [powerOfT] FROM vsop WHERE [body] =
%d",
body
);
HRESULT hr = S_OK;
_RecordsetPtr Rs = NULL;
_bstr_t Connect("DSN=Astronomy");
hr = Rs.CreateInstance(__uuidof(Recordset));
Rs->Open(
sql,
Connect,
adOpenForwardOnly,
adLockReadOnly,
adCmdUnknown
);
Rs->MoveFirst();
while ( VARIANT_FALSE == Rs->GetadoEOF() ) {
a = (double) Rs->Fields->GetItem( _variant_t("a"))->Value;
b = (double) Rs->Fields->GetItem( _variant_t("b"))->Value;
c = (double) Rs->Fields->GetItem( _variant_t("c"))->Value;
coord = (long) Rs->Fields->GetItem( _variant_t("coord"))->Value - 1;
powerOfT = (long) Rs->Fields->GetItem(
_variant_t("PowerOfT"))->Value;
storage[coord][powerOfT] += a * cos(b + c*t0);
Rs->MoveNext();
}
for (coord = 0; coord < 3; coord++)
{
results[coord] = 0;
for (powerOfT = 0; powerOfT < 6; powerOfT++)
results[coord] += storage[coord][powerOfT]*pow(t0,powerOfT);
}
Rs->Close();
Rs = NULL;
/* Return value */
return SUCCESS;
}
int makeQVsop(double startJD, double endJD, int N, int bodyID, double
results[3][6])
{
int X = 0;
int YX = 1;
int i = 0;
int j = 0;
int coord = 0;
int parity = 0;
double sumsX[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
double sumsYX[6][3] = {
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0}
};
double xarray[6][6] = {
{0,0,0,0,0,0},
{0,0,0,0,0,0},
{0,0,0,0,0,0},
{0,0,0,0,0,0},
{0,0,0,0,0,0},
{0,0,0,0,0,0}
};
double stepJD = (endJD - startJD)/(N-1);
double t = startJD;
int c = 1;
double t0 = 0;
double tpow[] = {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
double vsop[3] = { 0, 0, 0 };
while (t <= endJD)
{
for (coord = 0; coord < 3; coord++) vsop[coord]=0;
if (calcVsop(t, bodyID, vsop) == FAIL) return FAIL;
vsop[0] *= 180.0/PI;
vsop[1] *= 180.0/PI;
t0 = (t - startJD)/365.0;
sumsX[0]++;
for (j = 1; j < 11; j++)
tpow[j] = tpow[j-1] * t0; /* tpow[j] = t0^j = tpow[j-1] * t0 */
sumsX[j] += tpow[j]; /* sumsX[j] = Sum(t0^j) */
for (coord = 0; coord < 3; coord++)
{
for (j = 0; j < 6; j++)
sumsYX[j][coord] += tpow[j] * vsop[coord];
/* sumsYX[j][c] = vsop[c] * t0^j */
}
t += stepJD;
c += 1;
}
for (i = 0; i < 6; i++)
for (j = 0; j < 6; j++)
xarray[i][j] = sumsX[i+j];
/*
xarray = Sum(t0^...): index of sumsX
0 1 2 3 4 5
1 2 3 4 5 6
2 3 4 5 6 7
3 4 5 6 7 8
4 5 6 7 8 9
5 6 7 8 9 10
*/
if (linsys6x3(xarray, sumsYX) == FAIL)
return FAIL;
else
{
for (coord = 0; coord < 3; coord++)
for (i = 0; i < 6; i++)
results[coord][i] = sumsYX[i][coord];
return SUCCESS;
}
}
int linsys6x3(double a[6][6], double b[6][3])
{
int m = 3; /* cols in b (# coordinates) */
int n = 6; /* rows in a and b (# terms) */
int i, j, k;
double d;
/* reduce to upper triangular form */
for(i = 0; i < n; i++) /* for each row */
{
/* find maximal element */
j = i; /* store row index i in j as default */
for(k = i+1; k < n; k++) /* k = row index */
if(fabs(a[k][i]) > fabs(a[j][i])) /* use i as column index */
/* for shortcut */
j = k; /* update j with new row index */
/* j contains row index containing max element */
/* swap row i and j in both a and b */
for(k = i; k < n; k++)
{
d = a[i][k];
a[i][k] = a[j][k];
a[j][k] = d;
}
for(k = 0; k < m; k++)
{
d = b[i][k];
b[i][k] = b[j][k];
b[j][k] = d;
}
/* check if the system is singular */
/* if(fabs(a[i][i]) < 0.001)
{
// (void) fprintf(stderr,"Singular matrix, column %d\n", i+1);
return FAIL;
}
*/
/* elimination of column below diagonal */
for(j = i+1; j < n; j++) /* for each row j > i */
{
d = a[j][i]/a[i][i]; /* row subtraction factor */
a[j][i] = 0.; /* zero the element (result of exact subtraction) */
for(k = i+1; k < n; k++) /* for each col k > i */
a[j][k] -= d*a[i][k]; /* subtract factored row from a */
for(k = 0; k < m; k++) /* for each col k */
b[j][k] -= d*b[i][k]; /* subtract factored row from b */
}
}
/* normalization */
for(i = 0; i < n; i++) /* for each row */
{
for(j = i+1; j < n; j++) /* for each column right of the diagonal */
a[i][j] /= a[i][i]; /* divide by the diagonal */
for(j = 0; j < m; j++) /* for each col in b */
b[i][j] /= a[i][i]; /* divide by the diagonal */
a[i][i] = 1.; /* diagonal divided by itself is 1 */
}
/* reduction from upper trigonal form to diagonal form */
for(i = n-1; i >= 0; i--) /* for each row in reverse order */
{
for(j = 0; j < i; j++) /* for each row above i */
{
for(k = 0; k < m; k++) /* for each column in b */
b[j][k] -= b[i][k]*a[j][i]; /* subtract the row, factored by a */
a[j][i] = 0; /* exact subtraction */
}
}
return SUCCESS; /* a is identity matrix, b is solution */
}
More information about the Python-list
mailing list