Commit 6909b700 authored by Christoph Groth's avatar Christoph Groth
Browse files

initial version

parents
*~
*.pyc
*.so
/build
/dist
import tinyarray
import numpy
from time import time
################################################################
# Mock a numpy like "module" using Python tuples.
def tuple_array(seq, dtype):
return tuple(seq)
def tuple_zeros(shape, dtype):
return (0,) * shape
def tuple_dot(a, b):
return a[0] * b[0] + a[1] * b[1]
class Empty(object):
pass
tuples = Empty()
tuples.__name__ = 'tuples'
tuples.array = tuple_array
tuples.zeros = tuple_zeros
tuples.dot = tuple_dot
################################################################
def zeros(module, n=100000):
zeros = module.zeros
return list(zeros(2, int) for i in xrange(n))
def make_from_list(module, n=100000):
array = module.array
l = range(3)
return list(array(l, int) for i in xrange(n))
def dot(module, n=1000000):
dot = module.dot
a = module.zeros(2, int)
b = module.zeros(2, int)
for i in xrange(n):
c = dot(a, b)
def dot_tuple(module, n=100000):
dot = module.dot
a = module.zeros(2, int)
b = (0, 0)
for i in xrange(n):
c = dot(a, b)
def compare(function, modules):
print '****', function.__name__, '****'
for module in modules:
t = time()
try:
function(module)
except:
print "{0:15} failed".format(module.__name__)
else:
print '{0:15} {1:.4f} s'.format(module.__name__, time() - t)
def main():
modules = [tuples, tinyarray, numpy]
compare(zeros, modules)
compare(make_from_list, modules)
compare(dot, modules)
compare(dot_tuple, modules)
if __name__ == '__main__':
main()
#!/usr/bin/env python
from distutils.core import setup, Extension
module = Extension('tinyarray',
language = 'c++',
extra_compile_args = ['--std=c++0x'],
sources = ['src/main.cc', 'src/array.cc',
'src/functions.cc'])
setup (name = 'tinyarray',
version = '0.0',
ext_modules = [module])
#include <Python.h>
#include <limits>
#include <cassert>
#include "array.hh"
int dtype_converter(const PyObject *ob, Dtype *dtype)
{
if (ob == Py_None) {
*dtype = default_dtype;
} else if (ob == (PyObject *)(&PyFloat_Type)) {
*dtype = Dtype::DOUBLE;
} else if (ob == (PyObject *)(&PyComplex_Type)) {
*dtype = Dtype::COMPLEX;
} else if (ob == (PyObject *)(&PyInt_Type) ||
ob == (PyObject *)(&PyLong_Type)) {
*dtype = Dtype::LONG;
} else {
PyErr_SetString(PyExc_TypeError, "Invalid dtype.");
return 0;
}
return 1;
}
static Py_ssize_t load_seq_as_long(PyObject *obj, long *out, int maxlen)
{
assert(maxlen >= 1);
Py_ssize_t len;
if (PySequence_Check(obj)) {
obj = PySequence_Fast(obj, "Bug in tinyarray, load_seq_as_long");
if (!obj) return -1;
len = PySequence_Fast_GET_SIZE(obj);
if (len > maxlen) {
PyErr_Format(PyExc_ValueError, "Sequence too long."
" Maximum length is %d.", maxlen);
Py_DECREF(obj);
return -1;
}
for (PyObject **p = PySequence_Fast_ITEMS(obj), **e = p + len; p < e;
++p, ++out) {
*out = PyInt_AsLong(*p);
if (*out == -1 and PyErr_Occurred()) {
Py_DECREF(obj);
return -1;
}
}
} else {
len = 1;
*out = PyInt_AsLong(obj);
if (*out == -1 and PyErr_Occurred()) return -1;
}
return len;
}
Array *array_make(PyObject *shape_obj, Dtype dtype)
{
// Read shape.
long shape[max_ndim];
int ndim = load_seq_as_long(shape_obj, shape, max_ndim);
if (ndim == -1) return 0;
// Check shape and calculate ob_size, the total number of elements.
Py_ssize_t ob_size = 1;
#ifndef NO_OVERFLOW_DANGER
// `reserve' allows to detect overflow.
Py_ssize_t reserve = PY_SSIZE_T_MAX;
#endif
for (int d = 0; d < ndim; ++d) {
long elem = shape[d];
if (elem < 0) {
PyErr_SetString(PyExc_ValueError,
"Negative dimensions are not allowed.");
return 0;
}
if (elem > std::numeric_limits<Shape_t>::max()) {
PyErr_Format(PyExc_ValueError, "shape[%d] is too big.", d);
return 0;
}
#ifndef NO_OVERFLOW_DANGER
if (elem > reserve) {
PyErr_SetString(PyExc_ValueError, "Array would be too big.");
return 0;
}
#endif
ob_size *= elem;
#ifndef NO_OVERFLOW_DANGER
if (elem) reserve /= elem;
#endif
}
if (dtype != Dtype::LONG) {
PyErr_SetString(PyExc_ValueError, "dtype must be int for now");
return 0;
}
Array *result = PyObject_NewVar(Array, &Array_type, ob_size);
if (!result) return 0;
result->buffer_internal = 0;
result->ndim = ndim;
for (int d = 0; d < ndim; ++d) result->shape[d] = shape[d];
return result;
}
static void dealloc(Array *self)
{
if (self->buffer_internal) PyMem_Free(self->buffer_internal);
PyObject_Del(self);
}
static PyObject *repr(Array *self)
{
Py_ssize_t n = Py_SIZE(self);
char buf[30 * n + 200];
char *s = buf;
s += sprintf(s, "< ");
for (Py_ssize_t i = 0; i < n; ++i) {
if (i)
s += sprintf(s, ", %ld", self->ob_item[i]);
else
s += sprintf(s, "%ld", self->ob_item[i]);
}
s += sprintf(s, " shape=(");
for (int d = 0; d < self->ndim; ++d)
if (d)
s += sprintf(s, ", %d", self->shape[d]);
else
s += sprintf(s, "%d", self->shape[d]);
s += sprintf(s, ") >");
return PyString_FromString(buf);
}
PyDoc_STRVAR(doc, "array docstring: to be written\n");
static Py_ssize_t len(Array *self)
{
if (self->ndim == 0) {
PyErr_SetString(PyExc_TypeError, "len() of unsized object.");
return -1;
}
return self->shape[0];
}
static Py_ssize_t index_from_key(Array *self, PyObject *key)
{
long indices[max_ndim];
int ndim = load_seq_as_long(key, indices, max_ndim);
if (ndim == -1) {
PyErr_SetString(PyExc_IndexError,
"Invalid index.");
return -1;
}
if (ndim != self->ndim) {
PyErr_SetString(PyExc_IndexError, "Number of indices "
"must be equal to number of dimensions.");
return -1;
}
Py_ssize_t s = self->shape[0];
Py_ssize_t index = indices[0];
if (index < 0) index += s;
if (index < 0 || index >= s) {
PyErr_Format(PyExc_IndexError, "Index %ld out of range "
"(-%ld <= index < %ld) in dimension 0.",
indices[0], s, s);
return -1;
}
for (int d = 1; d < ndim; ++d) {
s = self->shape[d];
Py_ssize_t i = indices[d];
if (i < 0) i += s;
if (i < 0 || i >= s) {
PyErr_Format(PyExc_IndexError, "Index %ld out of range "
"(-%ld <= index < %ld) in dimension %d.",
indices[d], s, s, d);
return -1;
}
index *= s;
index += i;
}
return index;
}
static PyObject *getitem(Array *self, PyObject *key)
{
if (PySlice_Check(key)) {
PyErr_SetString(PyExc_NotImplementedError,
"slices are not implemented");
return 0;
} else {
Py_ssize_t index = index_from_key(self, key);
if (index == -1) return 0;
return PyInt_FromLong(self->ob_item[index]);
}
}
// The memoryview implementation is broken in Python <3.3. See
// http://bugs.python.org/issue7433 and http://bugs.python.org/issue10181.
// This means that the function bf_releasebuffer is useless for releasing
// ressources as it always gets the same copy of the view, even when called
// several times. To work around this, we allocate the needed memory in
// getbuffer and store a pointer to it in the array data structure. The memory
// gets freed only when the array is deleted.
//
// Once we require Python 3.3 we can return to a sane way of handling things,
// restore bf_releasebuffer and get rid of the buffer_internal field in struct
// Array.
int getbuffer(Array *self, Py_buffer *view, int flags)
{
assert(view);
if ((flags & PyBUF_F_CONTIGUOUS) == PyBUF_F_CONTIGUOUS) {
PyErr_SetString(PyExc_BufferError,
"Tinyarrays are not Fortran contiguous.");
goto fail;
}
if ((flags & PyBUF_WRITEABLE) == PyBUF_WRITEABLE) {
PyErr_SetString(PyExc_BufferError, "Tinyarrays are not writeable");
goto fail;
}
view->buf = self->ob_item;
view->itemsize = sizeof(long);
view->len = Py_SIZE(self) * view->itemsize;
view->readonly = 1;
if ((flags & PyBUF_FORMAT) == PyBUF_FORMAT)
view->format = (char*)"l";
else
view->format = 0;
if ((flags & PyBUF_ND) == PyBUF_ND) {
// From the documentation it's not clear whether it is allowed not to
// set strides to NULL (for C continuous arrays), but it works!
// However, there is a bug in current numpy
// (http://projects.scipy.org/numpy/ticket/2197) which requires strides
// if view->len == 0. Once this bug is fixed, we can remove the code
// which sets strides.
Py_ssize_t *shape, *strides;
if (!self->buffer_internal) {
if (view->len == 0 && (flags & PyBUF_STRIDES) == PyBUF_STRIDES) {
self->buffer_internal = shape = (Py_ssize_t*)
PyMem_Malloc(2 * self->ndim * sizeof(Py_ssize_t));
strides = shape + self->ndim;
} else {
self->buffer_internal = shape =
(Py_ssize_t*)PyMem_Malloc(self->ndim * sizeof(Py_ssize_t));
strides = 0;
}
if (!shape) {
PyErr_SetNone(PyExc_MemoryError);
goto fail;
}
for (int d = 0; d < self->ndim; ++d) shape[d] = self->shape[d];
if (strides) {
Py_ssize_t l = view->len;
for (int d = 0; d < self->ndim; ++d) {
if (shape[d]) l /= shape[d];
strides[d] = l;
}
}
} else {
shape = self->buffer_internal;
if (view->len == 0 && (flags & PyBUF_STRIDES) == PyBUF_STRIDES)
strides = shape + self->ndim;
else
strides = 0;
}
view->ndim = self->ndim;
view->shape = shape;
view->strides = strides;
} else {
view->ndim = 0;
view->shape = 0;
view->strides = 0;
}
view->internal = 0;
view->suboffsets = 0;
// Success.
Py_INCREF(self);
view->obj = (PyObject*)self;
return 0;
fail:
view->obj = 0;
return -1;
}
// This is modelled on Python's tuple hash function.
static long hash(Array *self)
{
long mult = 1000003, r = 0x345678;
Py_ssize_t len = Py_SIZE(self);
long *p = self->ob_item;
while (--len >= 0) {
r = (r ^ *p++) * mult;
mult += long(82520 + len + len);
}
r += 97531;
if (r == -1) r = -2;
return r;
}
static PyMappingMethods as_mapping = {
(lenfunc)len,
(binaryfunc)getitem,
0
};
static PyBufferProcs as_buffer = {
// We only support the new buffer protocol.
0, // bf_getreadbuffer
0, // bf_getwritebuffer
0, // bf_getsegcount
0, // bf_getcharbuffer
(getbufferproc)getbuffer, // bf_getbuffer
0 // bf_releasebuffer
};
PyTypeObject Array_type = {
PyVarObject_HEAD_INIT(&PyType_Type, 0)
"tinyarray.ndarray",
sizeof(Array) - sizeof(long), // tp_basicsize
sizeof(long), // tp_itemsize
(destructor)dealloc, // tp_dealloc
0, // tp_print
0, // tp_getattr
0, // tp_setattr
0, // tp_compare
(reprfunc)repr, // tp_repr
0, // tp_as_number
0, // tp_as_sequence
&as_mapping, // tp_as_mapping
(hashfunc)hash, // tp_hash
0, // tp_call
0, // tp_str
PyObject_GenericGetAttr, // tp_getattro
0, // tp_setattro
&as_buffer, // tp_as_buffer
Py_TPFLAGS_DEFAULT
| Py_TPFLAGS_HAVE_NEWBUFFER, // tp_flags
doc, // tp_doc
0, // tp_traverse
0, // tp_clear
// richcompare, // tp_richcompare
0, // tp_richcompare
0, // tp_weaklistoffset
// iter, // tp_iter
0, // tp_iter
0, // tp_iternext
0, // tp_methods
0, // tp_members
0, // tp_getset
0, // tp_base
0, // tp_dict
0, // tp_descr_get
0, // tp_descr_set
0, // tp_dictoffset
0, // tp_init
0, // tp_alloc
0, // tp_new
PyObject_Del, // tp_free
};
#ifndef ARRAY_HH
#define ARRAY_HH
#include <climits>
enum class Dtype : char {DOUBLE = 0, COMPLEX, LONG, NONE};
const Dtype default_dtype = Dtype::DOUBLE;
int dtype_converter(const PyObject *ob, Dtype *dtype);
const int max_ndim = 3;
typedef unsigned short Shape_t;
// Set NO_OVERFLOW_DANGER if the linear index of an array entry calculated for
// a valid shape can never overflow.
// KEEP THE FOLLOWING CONSISTENT WITH THE PREVIOUS DEFINITIONS.
#if ((((size_t)-1)>>1) >= 3 * USHRT_MAX)
#ifdef NDEBUG
#define NO_OVERFLOW_DANGER
#endif
#endif
struct Array {
PyObject_VAR_HEAD
Py_ssize_t *buffer_internal; // Can be removed for Python 3.3, see array.cc
Dtype dtype;
unsigned char ndim;
Shape_t shape[max_ndim];
long ob_item[1];
};
PyAPI_DATA(PyTypeObject) Array_type;
#define array_check_exact(op) (Py_TYPE(op) == &Array_type)
Array *array_make(PyObject *shape_obj, Dtype dtype);
#endif // !ARRAY_HH
#include <Python.h>
#include <limits>
#include "array.hh"
#include "functions.hh"
using namespace std;
static PyObject *empty(PyObject *module, PyObject *args)
{
PyObject *shape;
Dtype dtype = default_dtype;
if (!PyArg_ParseTuple(args, "O|O&", &shape, dtype_converter, &dtype))
return 0;
return (PyObject *)array_make(shape, dtype);
}
static PyObject *zeros(PyObject *module, PyObject *args)
{
Array *a = (Array*)empty(module, args);
if (a)
for (long *p = a->ob_item, *e = p + a->ob_size; p < e; ++p)
*p = 0;
return (PyObject *)a;
}
static PyObject *ones(PyObject *module, PyObject *args)
{
Array *a = (Array*)empty(module, args);
if (a)
for (long *p = a->ob_item, *e = p + a->ob_size; p < e; ++p)
*p = 1;
return (PyObject *)a;
}
static PyObject *identity(PyObject *module, PyObject *args)
{
long n;
Dtype dtype = default_dtype;
if (!PyArg_ParseTuple(args, "l|O&", &n, dtype_converter, &dtype))
return 0;
if (n < 0) {
PyErr_SetString(PyExc_ValueError,
"Negative dimensions are not allowed.");
return 0;
}
if (n > std::numeric_limits<Shape_t>::max()) {
PyErr_SetString(PyExc_ValueError,
"Requested size is too big.");
return 0;
}
Array *a = PyObject_NewVar(Array, &Array_type, n * n);
if (!a) return 0;
a->buffer_internal = 0;
a->ndim = 2;
a->shape[0] = a->shape[1] = n;
long *p = a->ob_item;
for (int i = 1; i < n; ++i) {
*p++ = 1;
for (long *e = p + n; p < e; ++p)
*p = 0;
}
if (n) *p = 1;
return (PyObject *)a;
}
static Array *array_from_sequence(PyObject *top_level_seq,
Dtype dtype=Dtype::NONE)
{
assert(PySequence_Check(top_level_seq));
int d = 0, ndim = 0;
long shape[max_ndim];
// seqs is the stack of sequences being processed, all returned by
// PySequence_Fast. ps[d] and es[d] are the begin and end of the elements
// of seqs[d - 1].
PyObject *seqs[max_ndim], **ps[max_ndim], **es[max_ndim];
es[0] = 1 + (ps[0] = &top_level_seq);
Py_ssize_t ob_size = 1;
long *dest = 0; // Init. only needed to suppress a warning.
Array *result = 0;
while (true) {
// See http://projects.scipy.org/numpy/ticket/2199.
const char *msg = "A sequence does not support sequence protocol - "
"this is probably due to a bug in numpy for 0-d arrays.";
seqs[d] = PySequence_Fast(*ps[d]++, msg);
if (!seqs[d]) {--d; goto fail_gracefully;}
Py_ssize_t len = PySequence_Fast_GET_SIZE(seqs[d]);
if (result) {
// Check that the length of current sequence agrees with the shape.
if (len != shape[d]) {
PyErr_SetString(PyExc_ValueError,
"Input has irregular shape.");
goto fail_gracefully;
}
} else {
// It's the first time we visit a sequence at this depth.
if (len > std::numeric_limits<Shape_t>::max()) {
PyErr_SetString(PyExc_ValueError,
"Source sequence too long for tinyarray.");
goto fail_gracefully;
}
shape[d] = len;
ob_size *= len;
++ndim;
if (ndim > max_ndim) {
PyErr_SetString(PyExc_ValueError, "Too many dimensions.");
goto fail_gracefully;
}
}
PyObject **p = PySequence_Fast_ITEMS(seqs[d]);
PyObject **e = p + len;
if (len && PySequence_Check(p[0])) {
if (d + 1 == ndim && result) {
PyErr_SetString(PyExc_ValueError,
"Input has irregular nesting depth.");
goto fail_gracefully;
}
++d;