rm refine module
Oli
3 years ago
parent
f152649da5
commit
3210b268aa
@ -1,200 +0,0 @@
|
||||
/* 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);
|
||||
}
|
||||
@ -1,27 +0,0 @@
|
||||
PyArray_NewFromDescr(<PyTypeObject *> np.ndarray, np.dtype('<H'), 2,dims, strides,pBuf, np.NPY_C_CONTIGUOUS, None)
|
||||
|
||||
PyArray_NewFromDescr( &PyArray_Type, PyArray_DescrFromType(m_table.datatypeid()), m_table.ndim(), shape, strides, data, NPY_ARRAY_WRITEABLE, nullptr));
|
||||
|
||||
|
||||
|
||||
|
||||
/******INFOS******/
|
||||
PyObject* PyArray_NewFromDescr(PyTypeObject* subtype, PyArray_Descr* descr, int nd, npy_intp* dims, npy_intp* strides, void* data, int flags, PyObject* obj)
|
||||
|
||||
|
||||
nd dimensions
|
||||
described by dims
|
||||
data-type descriptor of the new array is descr.
|
||||
|
||||
If subtype is of an array subclass instead of the base &PyArray_Type, then obj is the object to pass to the __array_finalize__ method of the subclass.
|
||||
|
||||
If data is NULL, then new memory will be allocated and flags can be non-zero to indicate a Fortran-style contiguous array. If data is not NULL, then it is assumed to point to the memory to be used for the array and the flags argument is used as the new flags for the array (except the state of NPY_OWNDATA, NPY_ARRAY_WRITEBACKIFCOPY and NPY_ARRAY_UPDATEIFCOPY flags of the new array will be reset).
|
||||
|
||||
In addition, if data is non-NULL, then strides can also be provided. If strides is NULL, then the array strides are computed as C-style contiguous (default) or Fortran-style contiguous (flags is nonzero for data = NULL or flags & NPY_ARRAY_F_CONTIGUOUS is nonzero non-NULL data). Any provided dims and strides are copied into newly allocated dimension and strides arrays for the new array object.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
https://searchcode.com/file/101698327/numpy/core/src/multiarray/compiled_base.c#l-1251
|
||||
Binary file not shown.
@ -1,96 +0,0 @@
|
||||
#include <Python.h>
|
||||
#include <numpy/arrayobject.h>
|
||||
|
||||
#define NBDIM 3
|
||||
|
||||
/* wrapped refine function */
|
||||
static PyObject* refine(PyObject* self, PyObject* args)
|
||||
{
|
||||
PyArrayObject* in_array;
|
||||
int rx, ry, rz;
|
||||
PyObject* out_array;
|
||||
NpyIter* in_iter;
|
||||
NpyIter* out_iter;
|
||||
NpyIter_IterNextFunc* in_iternext;
|
||||
NpyIter_IterNextFunc* out_iternext;
|
||||
|
||||
/* parse single numpy array argument */
|
||||
if (!PyArg_ParseTuple(args, "O!iii", &PyArray_Type, &in_array, &rz, &ry, &rx))
|
||||
return NULL;
|
||||
if (PyArray_NDIM(in_array) != NBDIM)
|
||||
return NULL;
|
||||
|
||||
PyArray_Descr out_descr;
|
||||
out_descr=*PyArray_DESCR(in_array);
|
||||
|
||||
|
||||
|
||||
/* construct the output array, like the input array */
|
||||
out_array = PyArray_NewFromDescr(&PyArray_Type, out_descr, NBDIM, dims, NULL, data, NPY_ARRAY_OWNDATA, NULL);
|
||||
if (out_array == NULL)
|
||||
return NULL;
|
||||
|
||||
/* create the iterators */
|
||||
in_iter = NpyIter_New(in_array, NPY_ITER_READONLY, NPY_KEEPORDER, NPY_NO_CASTING, NULL);
|
||||
if (in_iter == NULL)
|
||||
goto fail;
|
||||
|
||||
out_iter = NpyIter_New((PyArrayObject*)out_array, NPY_ITER_READWRITE, NPY_KEEPORDER, NPY_NO_CASTING, NULL);
|
||||
if (out_iter == NULL) {
|
||||
NpyIter_Deallocate(in_iter);
|
||||
goto fail;
|
||||
}
|
||||
|
||||
in_iternext = NpyIter_GetIterNext(in_iter, NULL);
|
||||
out_iternext = NpyIter_GetIterNext(out_iter, NULL);
|
||||
if (in_iternext == NULL || out_iternext == NULL) {
|
||||
NpyIter_Deallocate(in_iter);
|
||||
NpyIter_Deallocate(out_iter);
|
||||
goto fail;
|
||||
}
|
||||
|
||||
|
||||
|
||||
double** in_dataptr = (double**) NpyIter_GetDataPtrArray(in_iter);
|
||||
double** out_dataptr = (double**) NpyIter_GetDataPtrArray(out_iter);
|
||||
|
||||
/* iterate over the arrays */
|
||||
do {
|
||||
**out_dataptr = cos(**in_dataptr);
|
||||
} while(in_iternext(in_iter) && out_iternext(out_iter));
|
||||
|
||||
/* clean up and return the result */
|
||||
NpyIter_Deallocate(in_iter);
|
||||
NpyIter_Deallocate(out_iter);
|
||||
Py_INCREF(out_array);
|
||||
return out_array;
|
||||
|
||||
/* in case bad things happen */
|
||||
fail:
|
||||
Py_XDECREF(out_array);
|
||||
return NULL;
|
||||
}
|
||||
|
||||
/* define functions in module */
|
||||
static PyMethodDef RefineMethods[] =
|
||||
{
|
||||
{"refine", refine, METH_VARARGS, "Refines a numpy array"},
|
||||
{NULL, NULL, 0, NULL}
|
||||
};
|
||||
|
||||
static struct PyModuleDef crefineDef =
|
||||
{
|
||||
PyModuleDef_HEAD_INIT,
|
||||
"refine", "",
|
||||
-1,
|
||||
RefineMethods
|
||||
};
|
||||
|
||||
|
||||
import_array();
|
||||
PyMODINIT_FUNC
|
||||
PyInit_refine(void)
|
||||
{
|
||||
return PyModule_Create(&crefineDef);
|
||||
}
|
||||
|
||||
@ -1,96 +0,0 @@
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
#include <time.h>
|
||||
|
||||
|
||||
int refine(double* array, long nx, long ny, long nz, int rx, int ry, int rz, double* r)//////, FILE* pf)
|
||||
{
|
||||
long x, y, z, curr=0, cura=0, curraux=0, currnext=0, curraux2=0, currnext2=0;
|
||||
int i, j, k;
|
||||
///////double* r=(double*)malloc(sizeof(double)*nz*nx*ny*rz*rx*ry);
|
||||
///////if(r==(double*)NULL) puts("not enough memory output");
|
||||
/////////////////npywheader(rx*nx, ry*ny, rz*nz, pf);
|
||||
|
||||
|
||||
|
||||
clock_t t1, t2;
|
||||
t1=clock();
|
||||
|
||||
|
||||
|
||||
for(z=0;z<nz;z++)
|
||||
{
|
||||
curraux2=curr;
|
||||
currnext2=curr+rz*ny*ry*nx*rx;
|
||||
////////////////curr=0;
|
||||
for(y=0;y<ny;y++)
|
||||
{
|
||||
curraux=curr;
|
||||
currnext=curr+ry*nx*rx;
|
||||
for(x=0;x<nx;x++)
|
||||
{
|
||||
for(i=0;i<rx;i++)
|
||||
r[curr++]=array[cura];
|
||||
cura++;
|
||||
}
|
||||
for(curr;curr<currnext;curr++)
|
||||
r[curr]=r[curraux++];
|
||||
}
|
||||
for(curr;curr<currnext2;curr++)
|
||||
r[curr]=r[curraux2++];
|
||||
//////////////////fwrite(r, sizeof(double), nx*ny*rx*ry, pf);
|
||||
|
||||
}
|
||||
t2=clock();
|
||||
printf("%lf seconds.\n", ((double)(t2-t1))/CLOCKS_PER_SEC);
|
||||
|
||||
|
||||
|
||||
|
||||
free(r);
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
long i, j, k, nx, ny, nz, rx, ry, rz;
|
||||
double* a;
|
||||
double* r;
|
||||
if(argc!=7) return -1;
|
||||
nx=atol(argv[1]);
|
||||
ny=atol(argv[2]);
|
||||
nz=atol(argv[3]);
|
||||
rx=atol(argv[4]);
|
||||
ry=atol(argv[5]);
|
||||
rz=atol(argv[6]);
|
||||
a=(double*)malloc(sizeof(double)*nx*ny*nz);
|
||||
if(a==(double*)NULL) return -1;
|
||||
for(i=0;i<nx*ny*nz;i++)
|
||||
a[i]=(double)i;
|
||||
r=(double*)malloc(sizeof(double)*nx*rx*ny*ry*nz*rz);
|
||||
if(r==(double*)NULL) return -1;
|
||||
refine(a, nx, ny, nz, rx, ry, rz, r);
|
||||
/*
|
||||
for(i=0;i<nx*rx*ny*ry*nz*rz;i++)
|
||||
printf("%lf ", r[i]);
|
||||
puts("\n");
|
||||
for(k=0;k<rz*nz;k++){
|
||||
for(j=0;j<ry*ny;j++){
|
||||
for(i=0;i<rx*nx;i++)
|
||||
printf("%lf ", r[k*ny*ry*nx*rx+j*nx*rx+i]);
|
||||
puts("\n");}
|
||||
puts("\n\n\n");}
|
||||
*/
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
Binary file not shown.
@ -1,84 +0,0 @@
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
#include <time.h>
|
||||
|
||||
|
||||
int refine(double* array, long nx, long ny, long nz, long rx, long ry, long rz, double* r)
|
||||
{
|
||||
long x, y, z, cur, cura, curnextlayerx, curlayery, curnextlayery, curlayerz, curnextlayerz;
|
||||
// double* r=(double*)malloc(sizeof(double)*nx*rx*ny*ry*nz*rz);
|
||||
// if(r==(double*)NULL) return NULL;
|
||||
cur=0;
|
||||
cura=0;
|
||||
int i, rxint=(int)rx;
|
||||
|
||||
clock_t t1, t2;
|
||||
t1=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++)
|
||||
{
|
||||
/*for(i=0;i<rxint;i++)
|
||||
r[cur++]=array[cura];
|
||||
cura++;*/
|
||||
curnextlayerx=cur+rx;
|
||||
r[cur++]=array[cura++];
|
||||
while(cur<curnextlayerx)
|
||||
r[cur++]=r[cur-2];
|
||||
}
|
||||
while(cur<curnextlayery)
|
||||
r[cur++]=r[curlayery++];
|
||||
}
|
||||
while(cur<curnextlayerz)
|
||||
r[cur++]=r[curlayerz++];
|
||||
}
|
||||
|
||||
|
||||
t2=clock();
|
||||
printf("%lf seconds.\n", ((double)(t2-t1))/CLOCKS_PER_SEC);
|
||||
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
long i, j, k, nx, ny, nz, rx, ry, rz;
|
||||
double* a;
|
||||
double* r;
|
||||
if(argc!=7) return -1;
|
||||
nx=atol(argv[1]);
|
||||
ny=atol(argv[2]);
|
||||
nz=atol(argv[3]);
|
||||
rx=atol(argv[4]);
|
||||
ry=atol(argv[5]);
|
||||
rz=atol(argv[6]);
|
||||
a=(double*)malloc(sizeof(double)*nx*ny*nz);
|
||||
if(a==(double*)NULL) return -1;
|
||||
for(i=0;i<nx*ny*nz;i++)
|
||||
a[i]=(double)i;
|
||||
r=(double*)malloc(sizeof(double)*nx*rx*ny*ry*nz*rz);
|
||||
if(r==(double*)NULL) return -1;
|
||||
refine(a, nx, ny, nz, rx, ry, rz, r);
|
||||
/*
|
||||
for(i=0;i<nx*rx*ny*ry*nz*rz;i++)
|
||||
printf("%lf ", r[i]);
|
||||
puts("\n");
|
||||
for(k=0;k<rz*nz;k++){
|
||||
for(j=0;j<ry*ny;j++){
|
||||
for(i=0;i<rx*nx;i++)
|
||||
printf("%lf ", r[k*ny*ry*nx*rx+j*nx*rx+i]);
|
||||
puts("\n");}
|
||||
puts("\n\n\n");}
|
||||
*/
|
||||
return 0;
|
||||
}
|
||||
@ -1,7 +0,0 @@
|
||||
from distutils.core import setup, Extension
|
||||
|
||||
|
||||
module = Extension("refine", sources=["FINALrefine.c"])
|
||||
|
||||
|
||||
setup(ext_modules=[module])
|
||||
@ -1,12 +0,0 @@
|
||||
from time import time
|
||||
import numpy as np
|
||||
import refine
|
||||
|
||||
size = 420
|
||||
a = np.arange(size ** 3).astype("f8").reshape((size, size, size))
|
||||
ti = time()
|
||||
b = refine.refine(a, 2, 2, 2)
|
||||
tf = time()
|
||||
dt = tf - ti
|
||||
|
||||
raw_input("")
|
||||
@ -1,56 +0,0 @@
|
||||
import numpy as np
|
||||
import sys
|
||||
from refine import refine as ref
|
||||
|
||||
|
||||
def get_p(pn, pdir, pprefix):
|
||||
|
||||
p = np.load(pdir + pprefix + "0" + ".npy")
|
||||
for i in range(1, pn):
|
||||
p = np.concatenate((p, np.load(pdir + pprefix + str(i) + ".npy")), axis=0)
|
||||
|
||||
return p
|
||||
|
||||
|
||||
def get_k(pn, kdir, kprefix):
|
||||
|
||||
k = (np.load(kdir + kprefix + "0" + ".npy"))[1:-1, :, :]
|
||||
for i in range(1, pn):
|
||||
k = np.concatenate(
|
||||
(k, (np.load(kdir + kprefix + str(i) + ".npy"))[1:-1, :, :]), axis=0
|
||||
)
|
||||
return ref(k, 2, 2, 2)
|
||||
|
||||
|
||||
def kef(P, K, i, j, k, pbc):
|
||||
# tx=2*K[:,:,i]*K[:,:,i+1]/(K[:,:,i]+K[:,:,i+1])
|
||||
# ty=2*K[:,j,:]*K[:,j+1,:]/(K[:,j,:]+K[:,j+1,:])
|
||||
tz = 2 * K[k, :, :] * K[k + 1, :, :] / (K[k, :, :] + K[k + 1, :, :])
|
||||
|
||||
# qx=tx*(P[:,:,i+1]-P[:,:,i])
|
||||
# qy=ty*(P[:,j+1,:]-P[:,j,:])
|
||||
qz = -tz * (P[k + 1, :, :] - P[k, :, :])
|
||||
|
||||
kz = qz.sum() * (K.shape[0] + 1) / (pbc * K.shape[1] * K.shape[2])
|
||||
return kz
|
||||
|
||||
|
||||
def test(pn, kdir, pdir, kprefix, pprefix):
|
||||
|
||||
K = get_k(pn, kdir, kprefix)
|
||||
print(K.shape)
|
||||
P = get_p(pn, pdir, pprefix)
|
||||
print(P)
|
||||
print(P.shape)
|
||||
print(kef(P, K, 1, 1, 1, 1000))
|
||||
return
|
||||
|
||||
|
||||
pn = 4
|
||||
kdir = "./test/"
|
||||
pdir = "./test/"
|
||||
kprefix = "k"
|
||||
pprefix = "P"
|
||||
|
||||
|
||||
test(pn, kdir, pdir, kprefix, pprefix)
|
||||
Loading…
Reference in New Issue