[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