/* Wrap an array refinement function using the Numpy-C-API. */

#include <Python.h>
#include <numpy/arrayobject.h>
#define NDIM 3


/////
/*
int refine(double** in_dataptr, NpyIter* in_iter, NpyIter_IterNextFunc* in_iternext, long nz, long ny, long nx, long rz, long ry, long rx, double* r)
{
	long x, y, z, cur, curnextlayerx, curlayery, curnextlayery, curlayerz, curnextlayerz;
	cur=0;

	clock_t ti, tf;
	ti=clock();

	for(z=0;z<nz;z++)
	{
		curlayerz=cur;
		curnextlayerz=cur+rz*ny*ry*nx*rx;
		for(y=0;y<ny;y++)
		{
			curlayery=cur;
			curnextlayery=cur+ry*nx*rx;
			for(x=0;x<nx;x++)
			{
				curnextlayerx=cur+rx;
				r[cur++]=**in_dataptr;
				in_iternext(in_iter);
				while(cur<curnextlayerx)
					r[cur++]=r[cur-2];
			}
			while(cur<curnextlayery)
				r[cur++]=r[curlayery++];
		}
		while(cur<curnextlayerz)
			r[cur++]=r[curlayerz++];
	}

	tf=clock();
	printf("%lf\n",((double)(tf-ti))/CLOCKS_PER_SEC);

	return 0;
}

*/
/////




/////

int refine(double** in_dataptr, NpyIter* in_iter, NpyIter_IterNextFunc* in_iternext, long n0, long n1, long n2, long r0, long r1, long r2, double* r)
{
	long cur, i0, i1, i2, curlayer0, curnextlayer0, curlayer1, curnextlayer1, curnextlayer2;
	cur=0;

	//clock_t ti, tf;
	//ti=clock();

	for(i0=0;i0<n0;i0++)
	{
		curlayer0=cur;
		curnextlayer0=cur+r0*n1*r1*n2*r2;
		for(i1=0;i1<n1;i1++)
		{
			curlayer1=cur;
			curnextlayer1=cur+r1*n2*r2;
			for(i2=0;i2<n2;i2++)
			{
				curnextlayer2=cur+r2;
				r[cur++]=**in_dataptr;
				in_iternext(in_iter);
				while(cur<curnextlayer2)
					r[cur++]=r[cur-2];
			}
			while(cur<curnextlayer1)
				r[cur++]=r[curlayer1++];
		}
		while(cur<curnextlayer0)
			r[cur++]=r[curlayer0++];
	}

	//tf=clock();
	//printf("%lf\n",((double)(tf-ti))/CLOCKS_PER_SEC);

	return 0;
}


/////











/* Wrapped refinement function */
static PyObject* refineFunc(PyObject* self, PyObject* args)
{
	npy_intp r0, r1, r2, out_dims[NPY_MAXDIMS];
	PyArrayObject* in_array;
	PyObject* out_array;
	NpyIter* in_iter;
	NpyIter_IterNextFunc* in_iternext;
	double** in_dataptr;
	double* ref_array;

	/* parse single numpy array argument */
	if(!PyArg_ParseTuple(args, "O!lll", &PyArray_Type, &in_array, &r0, &r1, &r2)) return NULL;

	out_dims[0]=PyArray_DIMS(in_array)[0]*r0;
	out_dims[1]=PyArray_DIMS(in_array)[1]*r1;
	out_dims[2]=PyArray_DIMS(in_array)[2]*r2;

	ref_array=(double*)malloc(sizeof(double)*(long)out_dims[0]*(long)out_dims[1]*(long)out_dims[2]);
	if(ref_array==(double*)NULL) return NULL;





	/* create the iterators */
	in_iter=NpyIter_New(in_array, NPY_ITER_READONLY/*|NPY_ITER_EXTERNAL_LOOP*/, NPY_KEEPORDER, NPY_NO_CASTING, NULL);
	if(in_iter==NULL) goto fail;


	in_iternext=NpyIter_GetIterNext(in_iter, NULL);
	if(in_iternext==NULL)
	{
		NpyIter_Deallocate(in_iter);
		goto fail;
	}
	in_dataptr=(double**)NpyIter_GetDataPtrArray(in_iter);



	/* refine the array */
	refine	(
		in_dataptr,
		in_iter,
		in_iternext,
		(long)PyArray_DIMS(in_array)[0],
		(long)PyArray_DIMS(in_array)[1],
		(long)PyArray_DIMS(in_array)[2],
		(long)r0,
		(long)r1,
		(long)r2,
		ref_array
		);


	out_array = PyArray_SimpleNewFromData(NDIM,out_dims,NPY_DOUBLE,ref_array);

	if (out_array == NULL)
		return NULL;

	/* clean up and return the result */
	NpyIter_Deallocate(in_iter);
	//Py_INCREF(out_array);
	PyArray_ENABLEFLAGS(out_array, NPY_ARRAY_OWNDATA);
	return out_array;

	/* in case bad things happen */
	fail:
		Py_XDECREF(out_array);
		return NULL;
}




/* define functions in module */
static PyMethodDef RefineMethods[]=
{
	{"refine", refineFunc, METH_VARARGS, "Refines a 3D numpy-array."},
	{NULL, NULL, 0, NULL}
};

static struct PyModuleDef crefineDef =
{
    PyModuleDef_HEAD_INIT,
    "refine", "",
    -1,
    RefineMethods
};

PyMODINIT_FUNC
PyInit_refine(void)
{
	import_array();
	return PyModule_Create(&crefineDef);
}