[yt-svn] commit/yt: 3 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sat Mar 22 16:17:02 PDT 2014
3 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/1e8645ad2598/
Changeset: 1e8645ad2598
Branch: yt-3.0
User: MatthewTurk
Date: 2014-02-15 19:47:21
Summary: Adding ORIGAMI routines from Bridget Falck.
Affected #: 4 files
diff -r 13526c7b7b6d19787433919813cbf86a1df9c144 -r 1e8645ad2598b5c06bd2f6513332270104554cf9 yt/utilities/lib/origami.pyx
--- /dev/null
+++ b/yt/utilities/lib/origami.pyx
@@ -0,0 +1,45 @@
+"""
+This calls the ORIGAMI routines
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+from libc.stdlib cimport malloc, free
+
+cdef extern from "origami_tags.h":
+ int compute_tags(int ng, double boxsize, double **r, int npart,
+ unsigned char *m)
+
+def run_origami(np.ndarray[np.float64_t, ndim=1] pos_x,
+ np.ndarray[np.float64_t, ndim=1] pos_y,
+ np.ndarray[np.float64_t, ndim=1] pos_z,
+ double boxsize):
+ # We assume these have been passed in in the correct order and
+ # C-contiguous.
+ cdef int npart = pos_x.size
+ if npart == 1:
+ return np.zeros(1, dtype="uint8")
+ assert(sizeof(unsigned char) == sizeof(np.uint8_t))
+ assert(sizeof(double) == sizeof(np.float64_t))
+ cdef int ng = np.round(npart**(1./3))
+ assert(ng**3 == npart)
+ cdef double **r = <double **> malloc(sizeof(double *) * 3)
+ r[0] = <double *> pos_x.data
+ r[1] = <double *> pos_y.data
+ r[2] = <double *> pos_z.data
+ cdef np.ndarray[np.uint8_t, ndim=1] tags = np.zeros(npart, dtype="uint8")
+ cdef void *m = <void*> tags.data
+ compute_tags(ng, boxsize, r, npart, <unsigned char*> m)
+ free(r)
+ return tags
diff -r 13526c7b7b6d19787433919813cbf86a1df9c144 -r 1e8645ad2598b5c06bd2f6513332270104554cf9 yt/utilities/lib/origami_tags.c
--- /dev/null
+++ b/yt/utilities/lib/origami_tags.c
@@ -0,0 +1,220 @@
+// This code was originally written by Bridget Falck and Mark Neyrinck.
+// They have agreed to release it under the terms of the BSD license.
+#include "origami_tags.h"
+
+int isneg(int h) {
+ return (int)(h < 0);
+}
+
+int par(int i, int j, int k, int ng) {
+ return i + (j + k*ng)*ng;
+}
+
+int compute_tags(int ng, double boxsize, double **r, int np,
+ unsigned char *m) {
+ /* Note that the particles must be fed in according to the order specified in
+ * the README file */
+ double predict, xmin,xmax,ymin,ymax,zmin,zmax;
+
+ double negb2, b2;
+ int ng2,ng4, h, i,i2, x,y,z,nhalo,nhalo0,nhalo1,nhalo2,nhaloany;
+ unsigned char *m0,*m1,*m2, mn,m0n,m1n,m2n; /*Morphology tag */
+
+ double dx,d1,d2;
+
+ b2 = boxsize/2.;
+ negb2 = -boxsize/2.;
+ ng2=ng/2;
+ ng4=ng/4;
+
+ /* Boxsize should be the range in r, yielding a range 0-1 */
+ printf("%d particles\n",np);fflush(stdout);
+ xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+ m0 = (unsigned char *)malloc(np*sizeof(unsigned char)); /* for the diagonals */
+ m1 = (unsigned char *)malloc(np*sizeof(unsigned char));
+ m2 = (unsigned char *)malloc(np*sizeof(unsigned char));
+ for (i=0; i<np;i++) {
+ if (r[0][i]<xmin) xmin = r[0][i]; if (r[0][i]>xmax) xmax = r[0][i];
+ if (r[1][i]<ymin) ymin = r[1][i]; if (r[1][i]>ymax) ymax = r[1][i];
+ if (r[2][i]<zmin) zmin = r[2][i]; if (r[2][i]>zmax) zmax = r[2][i];
+
+ m[i] = 1;
+ m0[i] = 1;
+ m1[i] = 1;
+ m2[i] = 1;
+ }
+
+ if (m==NULL) {
+ printf("Morphology array cannot be allocated.\n");
+ return 1;
+ }
+ // printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
+
+ printf("Calculating ORIGAMI morphology.\n");
+
+ for (x=0; x<ng; x++){
+ // printf("%d\n",x);fflush(stdout);
+ for (y=0; y<ng; y++) {
+ for (z=0; z<ng; z++) {
+ i = par(x,y,z,ng);
+ /* First just along the Cartesian axes */
+ /* x-direction */
+ if (m[i] % 2 > 0) {
+ for (h=1; h<ng4; h++) {
+ i2 = par((x+h)%ng,y,z,ng);
+ dx = r[0][i2]-r[0][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("x:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ m[i] *= 2;
+ if (m[i2] % 2 > 0){
+ m[i2] *= 2;
+ }
+ }
+ }
+ }
+ if (m[i] % 3 > 0) {
+ for (h=1; h<ng4; h++) {
+ i2 = par(x,(y+h)%ng,z,ng);
+ dx = r[1][i2]-r[1][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("y:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ m[i] *= 3;
+ if (m[i2] % 3 > 0){
+ m[i2] *= 3;
+ }
+ break;
+ }
+ }
+ }
+ if (m[i] % 5 > 0) {
+ for (h=1; h<ng4; h++) {
+ i2 = par(x,y,(z+h)%ng,ng);
+ dx = r[2][i2]-r[2][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("z:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ m[i] *= 5;
+ if (m[i2] % 5 > 0){
+ m[i2] *= 5;
+ }
+ break;
+ }
+ }
+ }
+ // Now do diagonal directions
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(x,goodmod(y+h,ng),goodmod(z+h,ng),ng);
+ d1 = r[1][i2]-r[1][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m0[i] *= 2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(x,goodmod(y+h,ng),goodmod(z-h,ng),ng);
+ d1 = r[1][i2]-r[1][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m0[i] *= 3;
+ break;
+ }
+ }
+ // y
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),y,goodmod(z+h,ng),ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m1[i] *= 2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),y,goodmod(z-h,ng),ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m1[i] *= 3;
+ break;
+ }
+ }
+ // z
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),goodmod(y+h,ng),z,ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[1][i2]-r[1][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m2[i] *=2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),goodmod(y-h,ng),z,ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[1][i2]-r[1][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m2[i] *= 3;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ nhalo = 0;
+ nhalo0 = 0;
+ nhalo1 = 0;
+ nhalo2 = 0;
+ nhaloany = 0;
+ for (i=0;i<np;i++){
+ mn = (m[i]%2 == 0) + (m[i]%3 == 0) + (m[i]%5 == 0);
+ m0n = (unsigned char)(m[i]%2 == 0) + (unsigned char)(m0[i]%2 == 0) + (unsigned char)(m0[i]%3 == 0);
+ m1n = (unsigned char)(m[i]%3 == 0) + (unsigned char)(m1[i]%2 == 0) + (unsigned char)(m1[i]%3 == 0);
+ m2n = (unsigned char)(m[i]%5 == 0) + (unsigned char)(m2[i]%2 == 0) + (unsigned char)(m2[i]%3 == 0);
+ m[i] = max(mn,max(m0n,max(m1n,m2n)));
+ if (mn == 3) nhalo++;
+ if (m0n == 3) nhalo0++;
+ if (m1n == 3) nhalo1++;
+ if (m2n == 3) nhalo2++;
+ if (m[i] == 3) {
+ nhaloany++;
+ }
+ }
+ printf("nhalo=%d,%d,%d,%d,%d\n",nhalo,nhalo0,nhalo1,nhalo2,nhaloany);
+ free(m0);
+ free(m1);
+ free(m2);
+
+ /* Output m */
+ return 0;
+}
diff -r 13526c7b7b6d19787433919813cbf86a1df9c144 -r 1e8645ad2598b5c06bd2f6513332270104554cf9 yt/utilities/lib/origami_tags.h
--- /dev/null
+++ b/yt/utilities/lib/origami_tags.h
@@ -0,0 +1,14 @@
+#ifndef __ORIGAMI_TAGS_H__
+#define __ORIGAMI_TAGS_H__
+#include <stdio.h>
+#include <stdlib.h>
+
+#define BF 1e30
+#define max(A,B) (((A)>(B)) ? (A):(B))
+#define goodmod(A,B) (((A) >= (B)) ? (A-B):(((A) < 0) ? (A+B):(A)))
+
+int isneg(int h);
+int par(int i, int j, int k, int ng);
+int compute_tags(int ng, double boxsize, double **r, int npart,
+ unsigned char *m);
+#endif // __ORIGAMI_TAGS_H__
diff -r 13526c7b7b6d19787433919813cbf86a1df9c144 -r 1e8645ad2598b5c06bd2f6513332270104554cf9 yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -99,6 +99,10 @@
config.add_extension("Octree",
["yt/utilities/lib/Octree.pyx"],
libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ config.add_extension("origami",
+ ["yt/utilities/lib/origami.pyx",
+ "yt/utilities/lib/origami_tags.c"],
+ depends=["yt/utilities/lib/origami_tags.h"])
config.add_extension("image_utilities",
["yt/utilities/lib/image_utilities.pyx"],
libraries=["m"],
https://bitbucket.org/yt_analysis/yt/commits/bcfd8d77ee8d/
Changeset: bcfd8d77ee8d
Branch: yt-3.0
User: MatthewTurk
Date: 2014-03-19 13:53:02
Summary: Updating ORIGAMI with latest version from Bridget, adding citation.
Affected #: 3 files
diff -r 1e8645ad2598b5c06bd2f6513332270104554cf9 -r bcfd8d77ee8dc0c450c638eeddb71c4f7c1cc783 CITATION
--- a/CITATION
+++ b/CITATION
@@ -29,3 +29,28 @@
adsurl = {http://adsabs.harvard.edu/abs/2011ApJS..192....9T},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
+
+Using yt can also utilize other functionality. If you utilize ORIGAMI, we ask
+that you please cite the ORIGAMI paper:
+
+ at ARTICLE{2012ApJ...754..126F,
+ author = {{Falck}, B.~L. and {Neyrinck}, M.~C. and {Szalay}, A.~S.},
+ title = "{ORIGAMI: Delineating Halos Using Phase-space Folds}",
+ journal = {\apj},
+archivePrefix = "arXiv",
+ eprint = {1201.2353},
+ primaryClass = "astro-ph.CO",
+ keywords = {dark matter, galaxies: halos, large-scale structure of universe, methods: numerical},
+ year = 2012,
+ month = aug,
+ volume = 754,
+ eid = {126},
+ pages = {126},
+ doi = {10.1088/0004-637X/754/2/126},
+ adsurl = {http://adsabs.harvard.edu/abs/2012ApJ...754..126F},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+The main homepage for ORIGAMI can be found here:
+
+http://icg.port.ac.uk/~falckb/origami.html
diff -r 1e8645ad2598b5c06bd2f6513332270104554cf9 -r bcfd8d77ee8dc0c450c638eeddb71c4f7c1cc783 yt/utilities/lib/origami.pyx
--- a/yt/utilities/lib/origami.pyx
+++ b/yt/utilities/lib/origami.pyx
@@ -21,12 +21,19 @@
int compute_tags(int ng, double boxsize, double **r, int npart,
unsigned char *m)
+cdef int printed_citation = 0
+
def run_origami(np.ndarray[np.float64_t, ndim=1] pos_x,
np.ndarray[np.float64_t, ndim=1] pos_y,
np.ndarray[np.float64_t, ndim=1] pos_z,
double boxsize):
# We assume these have been passed in in the correct order and
# C-contiguous.
+ global printed_citation
+ if printed_citation == 0:
+ print "ORIGAMI was developed by Bridget Falck and Mark Neyrinck."
+ print "Please cite Falck, Neyrinck, & Szalay 2012, ApJ, 754, 2, 125."
+ printed_citation = 1
cdef int npart = pos_x.size
if npart == 1:
return np.zeros(1, dtype="uint8")
diff -r 1e8645ad2598b5c06bd2f6513332270104554cf9 -r bcfd8d77ee8dc0c450c638eeddb71c4f7c1cc783 yt/utilities/lib/origami_tags.c
--- a/yt/utilities/lib/origami_tags.c
+++ b/yt/utilities/lib/origami_tags.c
@@ -59,7 +59,6 @@
i = par(x,y,z,ng);
/* First just along the Cartesian axes */
/* x-direction */
- if (m[i] % 2 > 0) {
for (h=1; h<ng4; h++) {
i2 = par((x+h)%ng,y,z,ng);
dx = r[0][i2]-r[0][i];
@@ -67,14 +66,15 @@
if (dx > b2) dx -= boxsize;
if (dx < 0.) {
/*printf("x:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 2 > 0) {
m[i] *= 2;
+ }
if (m[i2] % 2 > 0){
m[i2] *= 2;
}
+ break;
}
}
- }
- if (m[i] % 3 > 0) {
for (h=1; h<ng4; h++) {
i2 = par(x,(y+h)%ng,z,ng);
dx = r[1][i2]-r[1][i];
@@ -82,15 +82,15 @@
if (dx > b2) dx -= boxsize;
if (dx < 0.) {
/*printf("y:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 3 > 0) {
m[i] *= 3;
+ }
if (m[i2] % 3 > 0){
m[i2] *= 3;
}
break;
}
}
- }
- if (m[i] % 5 > 0) {
for (h=1; h<ng4; h++) {
i2 = par(x,y,(z+h)%ng,ng);
dx = r[2][i2]-r[2][i];
@@ -98,14 +98,15 @@
if (dx > b2) dx -= boxsize;
if (dx < 0.) {
/*printf("z:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 5 > 0) {
m[i] *= 5;
+ }
if (m[i2] % 5 > 0){
m[i2] *= 5;
}
break;
}
}
- }
// Now do diagonal directions
for (h=1; h<ng4; h = -h + isneg(h)) {
i2 = par(x,goodmod(y+h,ng),goodmod(z+h,ng),ng);
https://bitbucket.org/yt_analysis/yt/commits/d3538d7d0240/
Changeset: d3538d7d0240
Branch: yt-3.0
User: samskillman
Date: 2014-03-23 00:16:56
Summary: Merged in MatthewTurk/yt/yt-3.0 (pull request #733)
Adding ORIGAMI routines from Bridget Falck.
Affected #: 5 files
diff -r 7305c9e6d46f05581859b3d95d80c912dd00b2d6 -r d3538d7d02407d920c0baa6090b5a3973196336d CITATION
--- a/CITATION
+++ b/CITATION
@@ -29,3 +29,28 @@
adsurl = {http://adsabs.harvard.edu/abs/2011ApJS..192....9T},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
+
+Using yt can also utilize other functionality. If you utilize ORIGAMI, we ask
+that you please cite the ORIGAMI paper:
+
+ at ARTICLE{2012ApJ...754..126F,
+ author = {{Falck}, B.~L. and {Neyrinck}, M.~C. and {Szalay}, A.~S.},
+ title = "{ORIGAMI: Delineating Halos Using Phase-space Folds}",
+ journal = {\apj},
+archivePrefix = "arXiv",
+ eprint = {1201.2353},
+ primaryClass = "astro-ph.CO",
+ keywords = {dark matter, galaxies: halos, large-scale structure of universe, methods: numerical},
+ year = 2012,
+ month = aug,
+ volume = 754,
+ eid = {126},
+ pages = {126},
+ doi = {10.1088/0004-637X/754/2/126},
+ adsurl = {http://adsabs.harvard.edu/abs/2012ApJ...754..126F},
+ adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+The main homepage for ORIGAMI can be found here:
+
+http://icg.port.ac.uk/~falckb/origami.html
diff -r 7305c9e6d46f05581859b3d95d80c912dd00b2d6 -r d3538d7d02407d920c0baa6090b5a3973196336d yt/utilities/lib/origami.pyx
--- /dev/null
+++ b/yt/utilities/lib/origami.pyx
@@ -0,0 +1,52 @@
+"""
+This calls the ORIGAMI routines
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+from libc.stdlib cimport malloc, free
+
+cdef extern from "origami_tags.h":
+ int compute_tags(int ng, double boxsize, double **r, int npart,
+ unsigned char *m)
+
+cdef int printed_citation = 0
+
+def run_origami(np.ndarray[np.float64_t, ndim=1] pos_x,
+ np.ndarray[np.float64_t, ndim=1] pos_y,
+ np.ndarray[np.float64_t, ndim=1] pos_z,
+ double boxsize):
+ # We assume these have been passed in in the correct order and
+ # C-contiguous.
+ global printed_citation
+ if printed_citation == 0:
+ print "ORIGAMI was developed by Bridget Falck and Mark Neyrinck."
+ print "Please cite Falck, Neyrinck, & Szalay 2012, ApJ, 754, 2, 125."
+ printed_citation = 1
+ cdef int npart = pos_x.size
+ if npart == 1:
+ return np.zeros(1, dtype="uint8")
+ assert(sizeof(unsigned char) == sizeof(np.uint8_t))
+ assert(sizeof(double) == sizeof(np.float64_t))
+ cdef int ng = np.round(npart**(1./3))
+ assert(ng**3 == npart)
+ cdef double **r = <double **> malloc(sizeof(double *) * 3)
+ r[0] = <double *> pos_x.data
+ r[1] = <double *> pos_y.data
+ r[2] = <double *> pos_z.data
+ cdef np.ndarray[np.uint8_t, ndim=1] tags = np.zeros(npart, dtype="uint8")
+ cdef void *m = <void*> tags.data
+ compute_tags(ng, boxsize, r, npart, <unsigned char*> m)
+ free(r)
+ return tags
diff -r 7305c9e6d46f05581859b3d95d80c912dd00b2d6 -r d3538d7d02407d920c0baa6090b5a3973196336d yt/utilities/lib/origami_tags.c
--- /dev/null
+++ b/yt/utilities/lib/origami_tags.c
@@ -0,0 +1,221 @@
+// This code was originally written by Bridget Falck and Mark Neyrinck.
+// They have agreed to release it under the terms of the BSD license.
+#include "origami_tags.h"
+
+int isneg(int h) {
+ return (int)(h < 0);
+}
+
+int par(int i, int j, int k, int ng) {
+ return i + (j + k*ng)*ng;
+}
+
+int compute_tags(int ng, double boxsize, double **r, int np,
+ unsigned char *m) {
+ /* Note that the particles must be fed in according to the order specified in
+ * the README file */
+ double predict, xmin,xmax,ymin,ymax,zmin,zmax;
+
+ double negb2, b2;
+ int ng2,ng4, h, i,i2, x,y,z,nhalo,nhalo0,nhalo1,nhalo2,nhaloany;
+ unsigned char *m0,*m1,*m2, mn,m0n,m1n,m2n; /*Morphology tag */
+
+ double dx,d1,d2;
+
+ b2 = boxsize/2.;
+ negb2 = -boxsize/2.;
+ ng2=ng/2;
+ ng4=ng/4;
+
+ /* Boxsize should be the range in r, yielding a range 0-1 */
+ printf("%d particles\n",np);fflush(stdout);
+ xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+ m0 = (unsigned char *)malloc(np*sizeof(unsigned char)); /* for the diagonals */
+ m1 = (unsigned char *)malloc(np*sizeof(unsigned char));
+ m2 = (unsigned char *)malloc(np*sizeof(unsigned char));
+ for (i=0; i<np;i++) {
+ if (r[0][i]<xmin) xmin = r[0][i]; if (r[0][i]>xmax) xmax = r[0][i];
+ if (r[1][i]<ymin) ymin = r[1][i]; if (r[1][i]>ymax) ymax = r[1][i];
+ if (r[2][i]<zmin) zmin = r[2][i]; if (r[2][i]>zmax) zmax = r[2][i];
+
+ m[i] = 1;
+ m0[i] = 1;
+ m1[i] = 1;
+ m2[i] = 1;
+ }
+
+ if (m==NULL) {
+ printf("Morphology array cannot be allocated.\n");
+ return 1;
+ }
+ // printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
+
+ printf("Calculating ORIGAMI morphology.\n");
+
+ for (x=0; x<ng; x++){
+ // printf("%d\n",x);fflush(stdout);
+ for (y=0; y<ng; y++) {
+ for (z=0; z<ng; z++) {
+ i = par(x,y,z,ng);
+ /* First just along the Cartesian axes */
+ /* x-direction */
+ for (h=1; h<ng4; h++) {
+ i2 = par((x+h)%ng,y,z,ng);
+ dx = r[0][i2]-r[0][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("x:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 2 > 0) {
+ m[i] *= 2;
+ }
+ if (m[i2] % 2 > 0){
+ m[i2] *= 2;
+ }
+ break;
+ }
+ }
+ for (h=1; h<ng4; h++) {
+ i2 = par(x,(y+h)%ng,z,ng);
+ dx = r[1][i2]-r[1][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("y:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 3 > 0) {
+ m[i] *= 3;
+ }
+ if (m[i2] % 3 > 0){
+ m[i2] *= 3;
+ }
+ break;
+ }
+ }
+ for (h=1; h<ng4; h++) {
+ i2 = par(x,y,(z+h)%ng,ng);
+ dx = r[2][i2]-r[2][i];
+ if (dx < negb2) dx += boxsize;
+ if (dx > b2) dx -= boxsize;
+ if (dx < 0.) {
+ /*printf("z:%d %d %d %d %f\n",x,y,z,h,dx);*/
+ if (m[i] % 5 > 0) {
+ m[i] *= 5;
+ }
+ if (m[i2] % 5 > 0){
+ m[i2] *= 5;
+ }
+ break;
+ }
+ }
+ // Now do diagonal directions
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(x,goodmod(y+h,ng),goodmod(z+h,ng),ng);
+ d1 = r[1][i2]-r[1][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m0[i] *= 2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(x,goodmod(y+h,ng),goodmod(z-h,ng),ng);
+ d1 = r[1][i2]-r[1][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m0[i] *= 3;
+ break;
+ }
+ }
+ // y
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),y,goodmod(z+h,ng),ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m1[i] *= 2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),y,goodmod(z-h,ng),ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[2][i2]-r[2][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m1[i] *= 3;
+ break;
+ }
+ }
+ // z
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),goodmod(y+h,ng),z,ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[1][i2]-r[1][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 + d2)*h < 0.) {
+ m2[i] *=2;
+ break;
+ }
+ }
+ for (h=1; h<ng4; h = -h + isneg(h)) {
+ i2 = par(goodmod(x+h,ng),goodmod(y-h,ng),z,ng);
+ d1 = r[0][i2]-r[0][i];
+ d2 = r[1][i2]-r[1][i];
+ if (d1 < negb2) d1 += boxsize;
+ if (d1 > b2) d1 -= boxsize;
+ if (d2 < negb2) d2 += boxsize;
+ if (d2 > b2) d2 -= boxsize;
+ if ((d1 - d2)*h < 0.) {
+ m2[i] *= 3;
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ nhalo = 0;
+ nhalo0 = 0;
+ nhalo1 = 0;
+ nhalo2 = 0;
+ nhaloany = 0;
+ for (i=0;i<np;i++){
+ mn = (m[i]%2 == 0) + (m[i]%3 == 0) + (m[i]%5 == 0);
+ m0n = (unsigned char)(m[i]%2 == 0) + (unsigned char)(m0[i]%2 == 0) + (unsigned char)(m0[i]%3 == 0);
+ m1n = (unsigned char)(m[i]%3 == 0) + (unsigned char)(m1[i]%2 == 0) + (unsigned char)(m1[i]%3 == 0);
+ m2n = (unsigned char)(m[i]%5 == 0) + (unsigned char)(m2[i]%2 == 0) + (unsigned char)(m2[i]%3 == 0);
+ m[i] = max(mn,max(m0n,max(m1n,m2n)));
+ if (mn == 3) nhalo++;
+ if (m0n == 3) nhalo0++;
+ if (m1n == 3) nhalo1++;
+ if (m2n == 3) nhalo2++;
+ if (m[i] == 3) {
+ nhaloany++;
+ }
+ }
+ printf("nhalo=%d,%d,%d,%d,%d\n",nhalo,nhalo0,nhalo1,nhalo2,nhaloany);
+ free(m0);
+ free(m1);
+ free(m2);
+
+ /* Output m */
+ return 0;
+}
diff -r 7305c9e6d46f05581859b3d95d80c912dd00b2d6 -r d3538d7d02407d920c0baa6090b5a3973196336d yt/utilities/lib/origami_tags.h
--- /dev/null
+++ b/yt/utilities/lib/origami_tags.h
@@ -0,0 +1,14 @@
+#ifndef __ORIGAMI_TAGS_H__
+#define __ORIGAMI_TAGS_H__
+#include <stdio.h>
+#include <stdlib.h>
+
+#define BF 1e30
+#define max(A,B) (((A)>(B)) ? (A):(B))
+#define goodmod(A,B) (((A) >= (B)) ? (A-B):(((A) < 0) ? (A+B):(A)))
+
+int isneg(int h);
+int par(int i, int j, int k, int ng);
+int compute_tags(int ng, double boxsize, double **r, int npart,
+ unsigned char *m);
+#endif // __ORIGAMI_TAGS_H__
diff -r 7305c9e6d46f05581859b3d95d80c912dd00b2d6 -r d3538d7d02407d920c0baa6090b5a3973196336d yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -99,6 +99,10 @@
config.add_extension("Octree",
["yt/utilities/lib/Octree.pyx"],
libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ config.add_extension("origami",
+ ["yt/utilities/lib/origami.pyx",
+ "yt/utilities/lib/origami_tags.c"],
+ depends=["yt/utilities/lib/origami_tags.h"])
config.add_extension("image_utilities",
["yt/utilities/lib/image_utilities.pyx"],
libraries=["m"],
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list