Commit df3f01fd authored by eric pellegrini's avatar eric pellegrini

add extensions

parent 7fa7ff83
import cython
cimport numpy as np
from numpy cimport ndarray
cdef extern from "math.h":
double floor(double x)
double ceil(double x)
double sqrt(double x)
cdef inline double round(double r):
return floor(r + 0.5) if (r > 0.0) else ceil(r - 0.5)
def distance_histogram(ndarray[np.float64_t, ndim=2] config not None,
ndarray[np.float64_t, ndim=2] cell not None,
ndarray[np.float64_t, ndim=2] rcell not None,
ndarray[np.int32_t, ndim=1] indexes not None,
ndarray[np.int32_t, ndim=1] molindex not None,
ndarray[np.int32_t, ndim=1] symbolindex not None,
ndarray[np.float64_t, ndim=3] hintra not None,
ndarray[np.float64_t, ndim=3] hinter not None,
ndarray[np.float64_t, ndim=2] scaleconfig not None,
float rmin,
float dr):
# This computes the intra and intermolecular distances histogram.
# The algorithm is a Pyrex adaptation of the FORTRAN implementation
# made by Miguel Angel Gonzalez (Institut Laue Langevin).
cdef double x, y, z, sdx, sdy, sdz, rx, ry, rz, r
cdef int i, j, bin, nbins
nbins = hinter.shape[2]
for 0 <= i < indexes.shape[0]:
x = config[i,0]
y = config[i,1]
z = config[i,2]
scaleconfig[i,0] = x*rcell[0,0] + y*rcell[0,1] + z*rcell[0,2]
scaleconfig[i,1] = x*rcell[1,0] + y*rcell[1,1] + z*rcell[1,2]
scaleconfig[i,2] = x*rcell[2,0] + y*rcell[2,1] + z*rcell[2,2]
for 0 <= i < indexes.shape[0] - 1:
sx = scaleconfig[i,0]
sy = scaleconfig[i,1]
sz = scaleconfig[i,2]
for i + 1 <= j < indexes.shape[0]:
sdx = scaleconfig[j,0] - sx
sdy = scaleconfig[j,1] - sy
sdz = scaleconfig[j,2] - sz
sdx -= round(sdx)
sdy -= round(sdy)
sdz -= round(sdz)
rx = sdx*cell[0,0] + sdy*cell[0,1] + sdz*cell[0,2]
ry = sdx*cell[1,0] + sdy*cell[1,1] + sdz*cell[1,2]
rz = sdx*cell[2,0] + sdy*cell[2,1] + sdz*cell[2,2]
r = sqrt(rx*rx + ry*ry + rz*rz)
bin = <int>((r-rmin)/dr)
if ( (bin < 0) or (bin >= nbins)):
continue
if molindex[i] == molindex[j]:
hintra[symbolindex[i],symbolindex[j],bin] += 1.0
else:
hinter[symbolindex[i],symbolindex[j],bin] += 1.0
# -*- coding: utf-8 -*-
#################################################
# fast calculation module of trajectory-viewer #
# Gael Goret for the Institut Laue-Langevin #
# gael.goret@ill.eu #
#################################################
import cython
cimport numpy as np
from numpy cimport ndarray
cdef extern from "math.h":
double sqrt(double x)
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
def cpt_cluster_connectivity(ndarray[np.float64_t, ndim=2] coords,
ndarray[np.float64_t, ndim=1] covRadii,
np.float32_t tolerance,
bonds):
'''
Compute the connectivity of an atom cluster.
'''
cdef int i, j, nbat
cdef float radius, distance
nbat = len(covRadii)
for i in range(nbat-1):
radiusi = covRadii[i] + tolerance
for j in range(i+1,nbat):
radius = radiusi + covRadii[j]
distance = (coords[i,0] - coords[j,0])**2 + (coords[i,1] - coords[j,1])**2 + (coords[i,2] - coords[j,2])**2
if distance <= (radius*radius):
bonds.append([i,j])
return bonds
\ No newline at end of file
cimport numpy as np
import numpy as np
from numpy cimport ndarray
from libcpp.vector cimport vector
from libcpp.pair cimport pair
from libcpp cimport bool
ctypedef pair[ int, vector[int] ] equiv_t
ctypedef vector[double] point
cdef inline bool exclude_from_bounded_3Dbox(double x,
double y,
double z,
double Xmin,
double Xmax,
double Ymin,
double Ymax,
double Zmin,
double Zmax,
double border):
return not ((Xmin-border <=x) and (x <= Xmax+border) and (Ymin-border <= y) and (y <= Ymax+border) and (Zmin-border <= z) and (z <= Zmax+border))
cdef inline bool exclude_from_bounded_2Dbox(double x,
double y,
double Xmin,
double Xmax,
double Ymin,
double Ymax,
double border):
return not ((Xmin-border <=x) and (x <= Xmax+border) and (Ymin-border <= y) and (y <= Ymax+border))
def mic_generator_2D(double[:,:] pts not None, double border, double[:] box_param):
cdef int i, m, new_id, nb_init_points, nb_final_points
cdef double x, y, Xmin, Xmax, Ymin, Ymax, j, k, n0, n1
cdef vector[equiv_t] equivalence
cdef vector[point] new_points
cdef equiv_t e
cdef point n
cdef double lol[3]
cdef ndarray[np.float64_t, ndim = 2] res
nb_init_points = pts.shape[0]
Xmax = box_param[0]
Ymax = box_param[1]
Xmin = 0
Ymin = 0
lol[0] = -1
lol[1] = 0
lol[2] = 1
for i in range(nb_init_points):
x = pts[i,0]
y = pts[i,1]
for j in lol[0:3]:
for k in lol[0:3]:
if j == 0 and k == 0:
continue
n0 = x+j*(Xmax-Xmin)
n1 = y+k*(Ymax-Ymin)
if exclude_from_bounded_2Dbox(n0, n1, Xmin, Xmax, Ymin, Ymax, border):
continue
n = point()
n.push_back(n0)
n.push_back(n1)
new_points.push_back(n)
new_id = new_points.size() + nb_init_points
for m in range(equivalence.size()):
if equivalence[m].first == i:
equivalence[m].second.push_back(new_id)
continue
e = equiv_t()
e.first=i
e.second.push_back(new_id)
equivalence.push_back(e)
d = {}
for i in range(equivalence.size()):
d[equivalence[i].first] = <list> equivalence[i].second
nb_final_points = pts.shape[0] + new_points.size()
res = np.zeros((nb_init_points + new_points.size() ,2), dtype = np.float64)
for i in range(nb_init_points):
res[i,0] = pts[i,0]
res[i,1] = pts[i,1]
for i in range(new_points.size()):
j = i + nb_init_points
res[j,0] = new_points[i][0]
res[j,1] = new_points[i][1]
return res, d
def mic_generator_3D(double[:,:] pts not None, double border, double[:] box_param):
cdef int i, m, new_id, nb_init_points
cdef double x, y, z, Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, j, k, l, n0, n1, n2
cdef vector[equiv_t] equivalence
cdef vector[point] new_points
cdef equiv_t e
cdef point n
cdef double lol[3]
cdef ndarray[np.float64_t, ndim = 2] res
nb_init_points = pts.shape[0]
Xmax = box_param[0]
Ymax = box_param[1]
Zmax = box_param[2]
Xmin = 0
Ymin = 0
Zmin = 0
lol[0] = -1
lol[1] = 0
lol[2] = 1
for i in range(nb_init_points):
x = pts[i,0]
y = pts[i,1]
z = pts[i,2]
for j in lol[0:3]:
for k in lol[0:3]:
for l in lol[0:3]:
if j == 0 and k == 0 and l == 0:
continue
n0 = x+j*(Xmax-Xmin)
n1 = y+k*(Ymax-Ymin)
n2 = z+l*(Zmax-Zmin)
if exclude_from_bounded_3Dbox(n0, n1, n2, Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, border):
continue
n = point()
n.push_back(n0)
n.push_back(n1)
n.push_back(n2)
new_points.push_back(n)
new_id = new_points.size() + nb_init_points
for m in range(equivalence.size()):
if equivalence[m].first == i:
equivalence[m].second.push_back(new_id)
continue
e = equiv_t()
e.first=i
e.second.push_back(new_id)
equivalence.push_back(e)
d = {}
for i in range(equivalence.size()):
d[equivalence[i].first] = <list> equivalence[i].second
res = np.zeros((nb_init_points + new_points.size() ,3), dtype = np.float64)
for i in range(nb_init_points):
res[i,0] = pts[i,0]
res[i,1] = pts[i,1]
res[i,2] = pts[i,2]
for i in range(new_points.size()):
j = i + nb_init_points
res[j,0] = new_points[i][0]
res[j,1] = new_points[i][1]
res[j,2] = new_points[i][2]
return res, d
import cython
cimport numpy as np
import numpy as np
from numpy cimport ndarray
cdef extern from "math.h":
double floor(double x)
double ceil(double x)
@cython.boundscheck(False)
@cython.wraparound(False)
def mt(ndarray[np.float64_t, ndim = 2] config not None,
ndarray[np.int32_t, ndim = 3] grid not None,
double resolution,ndarray[np.float64_t, ndim = 1] mini not None):
cdef int at, nbatom , atom
cdef double Xpos, Ypos, Zpos, i, j, k , mx, my, mz
mx = mini[0]
my = mini[1]
mz = mini[2]
nbatom = config.shape[0]
for atom in range(nbatom):
# The of atoms |i| in the current configuration.
i = floor((config[atom,0]-mx)/resolution)
j = floor((config[atom,1]-my)/resolution)
k = floor((config[atom,2]-mz)/resolution)
grid[i,j,k] += 1
This diff is collapsed.
# -*-cython-*-
"""
Qhull shared definitions, for use by other Cython modules
"""
#
# Copyright (C) Pauli Virtanen, 2010.
#
# Distributed under the same BSD license as Scipy.
#
cdef extern from "numpy/ndarrayobject.h":
cdef enum:
NPY_MAXDIMS
ctypedef struct DelaunayInfo_t:
int ndim
int npoints
int nsimplex
double *points
int *simplices
int *neighbors
double *equations
double *transform
int *vertex_to_simplex
double paraboloid_scale
double paraboloid_shift
double *max_bound
double *min_bound
int *vertex_neighbors_indices
int *vertex_neighbors_indptr
cdef int _get_delaunay_info(DelaunayInfo_t *, obj,
int compute_transform,
int compute_vertex_to_simplex,
int compute_vertex_neighbors) except -1
#
# N-D geometry
#
cdef int _barycentric_inside(int ndim, double *transform,
double *x, double *c, double eps) nogil
cdef void _barycentric_coordinate_single(int ndim, double *transform,
double *x, double *c, int i) nogil
cdef void _barycentric_coordinates(int ndim, double *transform,
double *x, double *c) nogil
#
# N+1-D geometry
#
cdef void _lift_point(DelaunayInfo_t *d, double *x, double *z) nogil
cdef double _distplane(DelaunayInfo_t *d, int isimplex, double *point) nogil
#
# Finding simplices
#
cdef int _is_point_fully_outside(DelaunayInfo_t *d, double *x, double eps) nogil
cdef int _find_simplex_bruteforce(DelaunayInfo_t *d, double *c, double *x,
double eps, double eps_broad) nogil
cdef int _find_simplex_directed(DelaunayInfo_t *d, double *c, double *x,
int *start, double eps, double eps_broad) nogil
cdef int _find_simplex(DelaunayInfo_t *d, double *c, double *x, int *start,
double eps, double eps_broad) nogil
# -*- cython -*-
"""
List of sets of integers, low-level C implementation
Works similarly as
setlist = [set() for j in range(n)]
but with integer values.
"""
cimport libc.stdlib
cimport numpy as np
import numpy as np
cdef struct setlist_t:
size_t n
size_t *sizes
size_t *alloc_sizes
int **sets
cdef inline int init(setlist_t *setlist, size_t n, size_t size_guess) except -1:
"""
Initialise a list of `n` sets with a given guessed size
"""
cdef int j
setlist.n = n
setlist.sets = <int**>libc.stdlib.malloc(sizeof(int*) * n)
if setlist.sets == NULL:
raise MemoryError
setlist.sizes = <size_t*>libc.stdlib.malloc(sizeof(size_t) * n)
if setlist.sizes == NULL:
libc.stdlib.free(<void*>setlist.sets)
raise MemoryError
setlist.alloc_sizes = <size_t*>libc.stdlib.malloc(sizeof(size_t) * n)
if setlist.alloc_sizes == NULL:
libc.stdlib.free(<void*>setlist.sets)
libc.stdlib.free(<void*>setlist.sizes)
raise MemoryError
for j in xrange(n):
setlist.sizes[j] = 0
setlist.alloc_sizes[j] = size_guess
setlist.sets[j] = <int*>libc.stdlib.malloc(sizeof(int) * size_guess)
return 0
cdef inline void free(setlist_t *setlist):
"""
Free the set list
"""
cdef int j
for j in xrange(setlist.n):
libc.stdlib.free(<void*>setlist.sets[j])
libc.stdlib.free(<void*>setlist.sets)
libc.stdlib.free(<void*>setlist.sizes)
libc.stdlib.free(<void*>setlist.alloc_sizes)
setlist.sets = NULL
setlist.sizes = NULL
setlist.alloc_sizes = NULL
setlist.n = 0
cdef inline int add(setlist_t *setlist, int n, int value) nogil:
"""
Add a value to set `n`
"""
cdef size_t i, sz
cdef int *p
if n < 0 or n >= setlist.n:
return 1
for i in xrange(setlist.sizes[n]):
if setlist.sets[n][i] == value:
return 0
if setlist.sizes[n] >= setlist.alloc_sizes[n]:
sz = 2*setlist.alloc_sizes[n] + 1
p = <int*>libc.stdlib.realloc(<void*>setlist.sets[n], sz * sizeof(int))
if p == NULL:
return -1
setlist.sets[n] = p
setlist.alloc_sizes[n] = sz
setlist.sets[n][setlist.sizes[n]] = value
setlist.sizes[n] += 1
return 0
cdef inline object tocsr(setlist_t *setlist):
"""
Convert list of sets to CSR format
Integers for set `i` reside in data[indices[i]:indices[i+1]]
Returns
-------
indices
CSR indices
data
CSR data
"""
cdef size_t i, j, pos
cdef size_t total_size
cdef np.ndarray[np.npy_int, ndim=1] indices, data
total_size = 0
for j in xrange(setlist.n):
total_size += setlist.sizes[j]
indices = np.empty((setlist.n+1,), dtype=np.intc)
data = np.empty((total_size,), dtype=np.intc)
pos = 0
for i in xrange(setlist.n):
indices[i] = pos
for j in xrange(setlist.sizes[i]):
data[pos] = setlist.sets[i][j]
pos += 1
indices[setlist.n] = pos
return indices, data
This diff is collapsed.
/*<html><pre> -<a href="qh-geom.htm"
>-------------------------------</a><a name="TOP">-</a>
geom.h
header file for geometric routines
see qh-geom.htm and geom.c
Copyright (c) 1993-2012 The Geometry Center.
$Id: //main/2011/qhull/src/libqhull/geom.h#3 $$Change: 1464 $
$DateTime: 2012/01/25 22:58:41 $$Author: bbarber $
*/
#ifndef qhDEFgeom
#define qhDEFgeom 1
#include "libqhull.h"
/* ============ -macros- ======================== */
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="fabs_">-</a>
fabs_(a)
returns the absolute value of a
*/
#define fabs_( a ) ((( a ) < 0 ) ? -( a ):( a ))
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="fmax_">-</a>
fmax_(a,b)
returns the maximum value of a and b
*/
#define fmax_( a,b ) ( ( a ) < ( b ) ? ( b ) : ( a ) )
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="fmin_">-</a>
fmin_(a,b)
returns the minimum value of a and b
*/
#define fmin_( a,b ) ( ( a ) > ( b ) ? ( b ) : ( a ) )
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="maximize_">-</a>
maximize_(maxval, val)
set maxval to val if val is greater than maxval
*/
#define maximize_( maxval, val ) { if (( maxval ) < ( val )) ( maxval )= ( val ); }
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="minimize_">-</a>
minimize_(minval, val)
set minval to val if val is less than minval
*/
#define minimize_( minval, val ) { if (( minval ) > ( val )) ( minval )= ( val ); }
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="det2_">-</a>
det2_(a1, a2,
b1, b2)
compute a 2-d determinate
*/
#define det2_( a1,a2,b1,b2 ) (( a1 )*( b2 ) - ( a2 )*( b1 ))
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="det3_">-</a>
det3_(a1, a2, a3,
b1, b2, b3,
c1, c2, c3)
compute a 3-d determinate
*/
#define det3_( a1,a2,a3,b1,b2,b3,c1,c2,c3 ) ( ( a1 )*det2_( b2,b3,c2,c3 ) \
- ( b1 )*det2_( a2,a3,c2,c3 ) + ( c1 )*det2_( a2,a3,b2,b3 ) )
/*-<a href="qh-geom.htm#TOC"
>--------------------------------</a><a name="dX">-</a>
dX( p1, p2 )
dY( p1, p2 )
dZ( p1, p2 )
given two indices into rows[],
compute the difference between X, Y, or Z coordinates
*/
#define dX( p1,p2 ) ( *( rows[p1] ) - *( rows[p2] ))
#define dY( p1,p2 ) ( *( rows[p1]+1 ) - *( rows[p2]+1 ))
#define dZ( p1,p2 ) ( *( rows[p1]+2 ) - *( rows[p2]+2 ))
#define dW( p1,p2 ) ( *( rows[p1]+3 ) - *( rows[p2]+3 ))
/*============= prototypes in alphabetical order, infrequent at end ======= */
void qh_backnormal(realT **rows, int numrow, int numcol, boolT sign, coordT *normal, boolT *nearzero);
void qh_distplane(pointT *point, facetT *facet, realT *dist);
facetT *qh_findbest(pointT *point, facetT *startfacet,
boolT bestoutside, boolT isnewfacets, boolT noupper,
realT *dist, boolT *isoutside, int *numpart);
facetT *qh_findbesthorizon(boolT ischeckmax, pointT *point,
facetT *startfacet, boolT noupper, realT *bestdist, int *numpart);
facetT *qh_findbestnew(pointT *point, facetT *startfacet, realT *dist,
boolT bestoutside, boolT *isoutside, int *numpart);
void qh_gausselim(realT **rows, int numrow, int numcol, boolT *sign, boolT *nearzero);
realT qh_getangle(pointT *vect1, pointT *vect2);
pointT *qh_getcenter(setT *vertices);
pointT *qh_getcentrum(facetT *facet);
realT qh_getdistance(facetT *facet, facetT *neighbor, realT *mindist, realT *maxdist);
void qh_normalize(coordT *normal, int dim, boolT toporient);
void qh_normalize2 (coordT *normal, int dim, boolT toporient,
realT *minnorm, boolT *ismin);
pointT *qh_projectpoint(pointT *point, facetT *facet, realT dist);
void qh_setfacetplane(facetT *newfacets);
void qh_sethyperplane_det(int dim, coordT **rows, coordT *point0,
boolT toporient, coordT *normal, realT *offset, boolT *nearzero);
void qh_sethyperplane_gauss(int dim, coordT **rows, pointT *point0,
boolT toporient, coordT *normal, coordT *offset, boolT *nearzero);
boolT qh_sharpnewfacets(void);
/*========= infrequently used code in geom2.c =============*/
coordT *qh_copypoints(coordT *points, int numpoints, int dimension);
void qh_crossproduct(int dim, realT vecA[3], realT vecB[3], realT vecC[3]);
realT qh_determinant(realT **rows, int dim, boolT *nearzero);
realT qh_detjoggle(pointT *points, int numpoints, int dimension);
void qh_detroundoff(void);
realT qh_detsimplex(pointT *apex, setT *points, int dim, boolT *nearzero);
realT qh_distnorm(int dim, pointT *point, pointT *normal, realT *offsetp);
realT qh_distround(int dimension, realT maxabs, realT maxsumabs);
realT qh_divzero(realT numer, realT denom, realT mindenom1, boolT *zerodiv);
realT qh_facetarea(facetT *facet);
realT qh_facetarea_simplex(int dim, coordT *apex, setT *vertices,
vertexT *notvertex, boolT toporient, coordT *normal, realT *offset);
pointT *qh_facetcenter(setT *vertices);
facetT *qh_findgooddist(pointT *point, facetT *facetA, realT *distp, facetT **facetlist);
void qh_getarea(facetT *facetlist);
boolT qh_gram_schmidt(int dim, realT **rows);
boolT qh_inthresholds(coordT *normal, realT *angle);
void qh_joggleinput(void);
realT *qh_maxabsval(realT *normal, int dim);
setT *qh_maxmin(pointT *points, int numpoints, int dimension);
realT qh_maxouter(void);
void qh_maxsimplex(int dim, setT *maxpoints, pointT *points, int numpoints, setT **simplex);
realT qh_minabsval(realT *normal, int dim);
int qh_mindiff(realT *vecA, realT *vecB, int dim);
boolT qh_orientoutside(facetT *facet);
void qh_outerinner(facetT *facet, realT *outerplane, realT *innerplane);
coordT qh_pointdist(pointT *point1, pointT *point2, int dim);
void qh_printmatrix(FILE *fp, const char *string, realT **rows, int numrow, int numcol);
void qh_printpoints(FILE *fp, const char *string, setT *points);
void qh_projectinput(void);
void qh_projectpoints(signed char *project, int n, realT *points,
int numpoints, int dim, realT *newpoints, int newdim);
void qh_rotateinput(realT **rows);
void qh_rotatepoints(realT *points, int numpoints, int dim, realT **rows);
void qh_scaleinput(void);
void qh_scalelast(coordT *points, int numpoints, int dim, coordT low,
coordT high, coordT newhigh);
void qh_scalepoints(pointT *points, int numpoints, int dim,
realT *newlows, realT *newhighs);
boolT qh_sethalfspace(int dim, coordT *coords, coordT **nextp,
coordT *normal, coordT *offset, coordT *feasible);
coordT *qh_sethalfspace_all(int dim, int count, coordT *halfspaces, pointT *feasible);
pointT *qh_voronoi_center(int dim, setT *points);
#endif /* qhDEFgeom */
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
/*<html><pre> -<a href="qh-io.htm"
>-------------------------------</a><a name="TOP">-</a>
io.h
declarations of Input/Output functions
see README, libqhull.h and io.c
Copyright (c) 1993-2012 The Geometry Center.
$Id: //main/2011/qhull/src/libqhull/io.h#3 $$Change: 1464 $
$DateTime: 2012/01/25 22:58:41 $$Author: bbarber $
*/
#ifndef qhDEFio
#define qhDEFio 1
#include "libqhull.h"
/*============ constants and flags ==================*/
/*-<a href="qh-io.htm#TOC"