[yt-svn] commit/yt: 5 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Sep 9 05:44:52 PDT 2013
5 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/99a037526e5d/
Changeset: 99a037526e5d
Branch: yt
User: xarthisius
Date: 2013-08-21 23:46:46
Summary: Move bottle,peewee,pydot,rocket to extern module
Affected #: 14 files
diff -r 7e42ac327db55ab369e31419af6f40276ca58fc7 -r 99a037526e5d1b5d8ec29b24146344d33f6f154e yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
+++ b/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
@@ -55,7 +55,7 @@
from yt.funcs import *
from yt.utilities.pykdtree import KDTree
-import yt.utilities.pydot as pydot
+import yt.extern.pydot as pydot
# We don't currently use this, but we may again find a use for it in the
# future.
diff -r 7e42ac327db55ab369e31419af6f40276ca58fc7 -r 99a037526e5d1b5d8ec29b24146344d33f6f154e yt/analysis_modules/halo_merger_tree/merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/merger_tree.py
+++ b/yt/analysis_modules/halo_merger_tree/merger_tree.py
@@ -36,7 +36,7 @@
HaloProfiler
from yt.convenience import load
from yt.utilities.logger import ytLogger as mylog
-import yt.utilities.pydot as pydot
+import yt.extern.pydot as pydot
from yt.utilities.spatial import cKDTree
from yt.utilities.parallel_tools.parallel_analysis_interface import \
ParallelDummy, \
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/dd8cfa5ab087/
Changeset: dd8cfa5ab087
Branch: yt
User: xarthisius
Date: 2013-09-06 08:41:27
Summary: Drop bundled delaunay traiangulation and use the one shipped with matplotlib
Affected #: 20 files
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -955,7 +955,7 @@
interpolated using the nearest neighbor method, with *side* points on a
side.
"""
- import yt.utilities.delaunay as de
+ import matplotlib.delaunay.triangulate as de
if log_spacing:
zz = np.log10(self[field])
else:
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/utilities/delaunay/LICENSE.txt
--- a/yt/utilities/delaunay/LICENSE.txt
+++ /dev/null
@@ -1,29 +0,0 @@
-Copyright (c) 2001, 2002 Enthought, Inc.
-
-All rights reserved.
-
-Redistribution and use in source and binary forms, with or without
-modification, are permitted provided that the following conditions are met:
-
- a. Redistributions of source code must retain the above copyright notice,
- this list of conditions and the following disclaimer.
- b. Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
- c. Neither the name of the Enthought nor the names of its contributors
- may be used to endorse or promote products derived from this software
- without specific prior written permission.
-
-
-THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
-ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
-OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/utilities/delaunay/README.txt
--- a/yt/utilities/delaunay/README.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-The code in this directory is subject to the license as specified in
-LICENSE.txt. It was primarily written by Robert Kern, and is formally part of
-the SciKits project, a sub-project of SciPy.
-
-For more information, please see __init__.py, or the SciPy website at
-www.scipy.org .
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/utilities/delaunay/VoronoiDiagramGenerator.cpp
--- a/yt/utilities/delaunay/VoronoiDiagramGenerator.cpp
+++ /dev/null
@@ -1,1152 +0,0 @@
-/*
- * The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
- * Bell Laboratories.
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- */
-
-/*
- * This code was originally written by Stephan Fortune in C code. Shane O'Sullivan,
- * have since modified it, encapsulating it in a C++ class and, fixing memory leaks and
- * adding accessors to the Voronoi Edges.
- * Permission to use, copy, modify, and distribute this software for any
- * purpose without fee is hereby granted, provided that this entire notice
- * is included in all copies of any software which is or includes a copy
- * or modification of this software and in all copies of the supporting
- * documentation for such software.
- * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
- * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
- * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
- * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
- */
-
-/*
- * Subsequently, Robert Kern modified it to yield Python objects.
- * Copyright 2005 Robert Kern <robert.kern at gmail.com>
- * See LICENSE.txt in the scipy source directory.
- */
-
-#include "VoronoiDiagramGenerator.h"
-
-VoronoiDiagramGenerator::VoronoiDiagramGenerator()
-{
- siteidx = 0;
- sites = 0;
-
- allMemoryList = new FreeNodeArrayList;
- allMemoryList->memory = 0;
- allMemoryList->next = 0;
- currentMemoryBlock = allMemoryList;
- allEdges = 0;
- allEdgeList = 0;
- iteratorEdges = 0;
- iterEdgeList = 0;
- minDistanceBetweenSites = 0;
-}
-
-VoronoiDiagramGenerator::~VoronoiDiagramGenerator()
-{
- cleanupEdgeList();
- cleanup();
- cleanupEdges();
-
- if(allMemoryList != 0)
- delete allMemoryList;
-}
-
-
-
-bool VoronoiDiagramGenerator::generateVoronoi(double *xValues, double *yValues, int numPoints, double minX, double maxX, double minY, double maxY, double minDist)
-{
- cleanupEdgeList();
- cleanup();
- cleanupEdges();
- int i;
-
- minDistanceBetweenSites = minDist;
-
- nsites=numPoints;
- plot = 0;
- triangulate = 0;
- debug = 1;
- sorted = 0;
- freeinit(&sfl, sizeof (Site));
-
- sites = (struct Site *) myalloc(nsites*sizeof( *sites));
-
- if(sites == 0)
- return false;
-
- xmin = xValues[0];
- ymin = yValues[0];
- xmax = xValues[0];
- ymax = yValues[0];
-
- for(i = 0; i< nsites; i++)
- {
- sites[i].coord.x = xValues[i];
- sites[i].coord.y = yValues[i];
- sites[i].sitenbr = i;
- sites[i].refcnt = 0;
-
- if(xValues[i] < xmin)
- xmin = xValues[i];
- else if(xValues[i] > xmax)
- xmax = xValues[i];
-
- if(yValues[i] < ymin)
- ymin = yValues[i];
- else if(yValues[i] > ymax)
- ymax = yValues[i];
-
- //printf("\n%f %f\n",xValues[i],yValues[i]);
- }
-
- qsort(sites, nsites, sizeof (*sites), scomp);
-
- siteidx = 0;
- geominit();
- double temp = 0;
- if(minX > maxX)
- {
- temp = minX;
- minX = maxX;
- maxX = temp;
- }
- if(minY > maxY)
- {
- temp = minY;
- minY = maxY;
- maxY = temp;
- }
- borderMinX = minX;
- borderMinY = minY;
- borderMaxX = maxX;
- borderMaxY = maxY;
-
- siteidx = 0;
-
- voronoi(triangulate);
-
- return true;
-}
-
-bool VoronoiDiagramGenerator::ELinitialize()
-{
- int i;
- freeinit(&hfl, sizeof **ELhash);
- ELhashsize = 2 * sqrt_nsites;
- ELhash = (struct Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);
-
- if(ELhash == 0)
- return false;
-
- for(i=0; i<ELhashsize; i +=1) ELhash[i] = (struct Halfedge *)NULL;
- ELleftend = HEcreate( (struct Edge *)NULL, 0);
- ELrightend = HEcreate( (struct Edge *)NULL, 0);
- ELleftend -> ELleft = (struct Halfedge *)NULL;
- ELleftend -> ELright = ELrightend;
- ELrightend -> ELleft = ELleftend;
- ELrightend -> ELright = (struct Halfedge *)NULL;
- ELhash[0] = ELleftend;
- ELhash[ELhashsize-1] = ELrightend;
-
- return true;
-}
-
-
-struct Halfedge* VoronoiDiagramGenerator::HEcreate(struct Edge *e,int pm)
-{
- struct Halfedge *answer;
- answer = (struct Halfedge *) getfree(&hfl);
- answer -> ELedge = e;
- answer -> ELpm = pm;
- answer -> PQnext = (struct Halfedge *) NULL;
- answer -> vertex = (struct Site *) NULL;
- answer -> ELrefcnt = 0;
- return(answer);
-}
-
-
-void VoronoiDiagramGenerator::ELinsert(struct Halfedge *lb, struct Halfedge *newHe)
-{
- newHe -> ELleft = lb;
- newHe -> ELright = lb -> ELright;
- (lb -> ELright) -> ELleft = newHe;
- lb -> ELright = newHe;
-}
-
-/* Get entry from hash table, pruning any deleted nodes */
-struct Halfedge * VoronoiDiagramGenerator::ELgethash(int b)
-{
- struct Halfedge *he;
-
- if(b<0 || b>=ELhashsize)
- return((struct Halfedge *) NULL);
- he = ELhash[b];
- if (he == (struct Halfedge *) NULL || he->ELedge != (struct Edge *) DELETED )
- return (he);
-
- /* Hash table points to deleted half edge. Patch as necessary. */
- ELhash[b] = (struct Halfedge *) NULL;
- if ((he -> ELrefcnt -= 1) == 0)
- makefree((Freenode*)he, &hfl);
- return ((struct Halfedge *) NULL);
-}
-
-struct Halfedge * VoronoiDiagramGenerator::ELleftbnd(struct Point *p)
-{
- int i, bucket;
- struct Halfedge *he;
-
- /* Use hash table to get close to desired halfedge */
- bucket = (int)((p->x - xmin)/deltax * ELhashsize); //use the hash function to find the place in the hash map that this HalfEdge should be
-
- if(bucket<0) bucket =0; //make sure that the bucket position in within the range of the hash array
- if(bucket>=ELhashsize) bucket = ELhashsize - 1;
-
- he = ELgethash(bucket);
- if(he == (struct Halfedge *) NULL) //if the HE isn't found, search backwards and forwards in the hash map for the first non-null entry
- {
- for(i=1; 1 ; i += 1)
- {
- if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL)
- break;
- if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL)
- break;
- };
- totalsearch += i;
- };
- ntry += 1;
- /* Now search linear list of halfedges for the correct one */
- if (he==ELleftend || (he != ELrightend && right_of(he,p)))
- {
- do
- {
- he = he -> ELright;
- } while (he!=ELrightend && right_of(he,p)); //keep going right on the list until either the end is reached, or you find the 1st edge which the point
- he = he -> ELleft; //isn't to the right of
- }
- else //if the point is to the left of the HalfEdge, then search left for the HE just to the left of the point
- do
- {
- he = he -> ELleft;
- } while (he!=ELleftend && !right_of(he,p));
-
- /* Update hash table and reference counts */
- if(bucket > 0 && bucket <ELhashsize-1)
- {
- if(ELhash[bucket] != (struct Halfedge *) NULL)
- {
- ELhash[bucket] -> ELrefcnt -= 1;
- }
- ELhash[bucket] = he;
- ELhash[bucket] -> ELrefcnt += 1;
- };
- return (he);
-}
-
-
-/* This delete routine can't reclaim node, since pointers from hash
-table may be present. */
-void VoronoiDiagramGenerator::ELdelete(struct Halfedge *he)
-{
- (he -> ELleft) -> ELright = he -> ELright;
- (he -> ELright) -> ELleft = he -> ELleft;
- he -> ELedge = (struct Edge *)DELETED;
-}
-
-
-struct Halfedge * VoronoiDiagramGenerator::ELright(struct Halfedge *he)
-{
- return (he -> ELright);
-}
-
-struct Halfedge * VoronoiDiagramGenerator::ELleft(struct Halfedge *he)
-{
- return (he -> ELleft);
-}
-
-
-struct Site * VoronoiDiagramGenerator::leftreg(struct Halfedge *he)
-{
- if(he -> ELedge == (struct Edge *)NULL)
- return(bottomsite);
- return( he -> ELpm == le ?
- he -> ELedge -> reg[le] : he -> ELedge -> reg[re]);
-}
-
-struct Site * VoronoiDiagramGenerator::rightreg(struct Halfedge *he)
-{
- if(he -> ELedge == (struct Edge *)NULL) //if this halfedge has no edge, return the bottom site (whatever that is)
- return(bottomsite);
-
- //if the ELpm field is zero, return the site 0 that this edge bisects, otherwise return site number 1
- return( he -> ELpm == le ? he -> ELedge -> reg[re] : he -> ELedge -> reg[le]);
-}
-
-void VoronoiDiagramGenerator::geominit()
-{
- double sn;
-
- freeinit(&efl, sizeof(Edge));
- nvertices = 0;
- nedges = 0;
- sn = (double)nsites+4;
- sqrt_nsites = (int)sqrt(sn);
- deltay = ymax - ymin;
- deltax = xmax - xmin;
-}
-
-
-struct Edge * VoronoiDiagramGenerator::bisect(struct Site *s1, struct Site *s2)
-{
- double dx,dy,adx,ady;
- struct Edge *newedge;
-
- newedge = (struct Edge *) getfree(&efl);
-
- newedge -> reg[0] = s1; //store the sites that this edge is bisecting
- newedge -> reg[1] = s2;
- ref(s1);
- ref(s2);
- newedge -> ep[0] = (struct Site *) NULL; //to begin with, there are no endpoints on the bisector - it goes to infinity
- newedge -> ep[1] = (struct Site *) NULL;
-
- dx = s2->coord.x - s1->coord.x; //get the difference in x dist between the sites
- dy = s2->coord.y - s1->coord.y;
- adx = dx>0 ? dx : -dx; //make sure that the difference in positive
- ady = dy>0 ? dy : -dy;
- newedge -> c = (double)(s1->coord.x * dx + s1->coord.y * dy + (dx*dx + dy*dy)*0.5);//get the slope of the line
-
- if (adx>ady)
- {
- newedge -> a = 1.0; newedge -> b = dy/dx; newedge -> c /= dx;//set formula of line, with x fixed to 1
- }
- else
- {
- newedge -> b = 1.0; newedge -> a = dx/dy; newedge -> c /= dy;//set formula of line, with y fixed to 1
- };
-
- newedge -> edgenbr = nedges;
-
- //printf("\nbisect(%d) ((%f,%f) and (%f,%f)",nedges,s1->coord.x,s1->coord.y,s2->coord.x,s2->coord.y);
-
- nedges += 1;
- return(newedge);
-}
-
-//create a new site where the HalfEdges el1 and el2 intersect - note that the Point in the argument list is not used, don't know why it's there
-struct Site * VoronoiDiagramGenerator::intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p)
-{
- struct Edge *e1,*e2, *e;
- struct Halfedge *el;
- double d, xint, yint;
- int right_of_site;
- struct Site *v;
-
- e1 = el1 -> ELedge;
- e2 = el2 -> ELedge;
- if(e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL)
- return ((struct Site *) NULL);
-
- //if the two edges bisect the same parent, return null
- if (e1->reg[1] == e2->reg[1])
- return ((struct Site *) NULL);
-
- d = e1->a * e2->b - e1->b * e2->a;
- if (-1.0e-10<d && d<1.0e-10)
- return ((struct Site *) NULL);
-
- xint = (e1->c*e2->b - e2->c*e1->b)/d;
- yint = (e2->c*e1->a - e1->c*e2->a)/d;
-
- if( (e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
- (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
- e1->reg[1]->coord.x < e2->reg[1]->coord.x) )
- {
- el = el1;
- e = e1;
- }
- else
- {
- el = el2;
- e = e2;
- };
-
- right_of_site = xint >= e -> reg[1] -> coord.x;
- if ((right_of_site && el -> ELpm == le) || (!right_of_site && el -> ELpm == re))
- return ((struct Site *) NULL);
-
- //create a new site at the point of intersection - this is a new vector event waiting to happen
- v = (struct Site *) getfree(&sfl);
- v -> refcnt = 0;
- v -> coord.x = xint;
- v -> coord.y = yint;
- return(v);
-}
-
-/* returns 1 if p is to right of halfedge e */
-int VoronoiDiagramGenerator::right_of(struct Halfedge *el,struct Point *p)
-{
- struct Edge *e;
- struct Site *topsite;
- int right_of_site, above, fast;
- double dxp, dyp, dxs, t1, t2, t3, yl;
-
- e = el -> ELedge;
- topsite = e -> reg[1];
- right_of_site = p -> x > topsite -> coord.x;
- if(right_of_site && el -> ELpm == le) return(1);
- if(!right_of_site && el -> ELpm == re) return (0);
-
- if (e->a == 1.0)
- { dyp = p->y - topsite->coord.y;
- dxp = p->x - topsite->coord.x;
- fast = 0;
- if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
- { above = dyp>= e->b*dxp;
- fast = above;
- }
- else
- { above = p->x + p->y*e->b > e-> c;
- if(e->b<0.0) above = !above;
- if (!above) fast = 1;
- };
- if (!fast)
- { dxs = topsite->coord.x - (e->reg[0])->coord.x;
- above = e->b * (dxp*dxp - dyp*dyp) <
- dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
- if(e->b<0.0) above = !above;
- };
- }
- else /*e->b==1.0 */
- { yl = e->c - e->a*p->x;
- t1 = p->y - yl;
- t2 = p->x - topsite->coord.x;
- t3 = yl - topsite->coord.y;
- above = t1*t1 > t2*t2 + t3*t3;
- };
- return (el->ELpm==le ? above : !above);
-}
-
-
-void VoronoiDiagramGenerator::endpoint(struct Edge *e,int lr,struct Site * s)
-{
- e -> ep[lr] = s;
- ref(s);
- if(e -> ep[re-lr]== (struct Site *) NULL)
- return;
-
- clip_line(e);
-
- deref(e->reg[le]);
- deref(e->reg[re]);
- makefree((Freenode*)e, &efl);
-}
-
-
-double VoronoiDiagramGenerator::dist(struct Site *s,struct Site *t)
-{
- double dx,dy;
- dx = s->coord.x - t->coord.x;
- dy = s->coord.y - t->coord.y;
- return (double)(sqrt(dx*dx + dy*dy));
-}
-
-
-void VoronoiDiagramGenerator::makevertex(struct Site *v)
-{
- v -> sitenbr = nvertices;
- nvertices += 1;
- out_vertex(v);
-}
-
-
-void VoronoiDiagramGenerator::deref(struct Site *v)
-{
- v -> refcnt -= 1;
- if (v -> refcnt == 0 )
- makefree((Freenode*)v, &sfl);
-}
-
-void VoronoiDiagramGenerator::ref(struct Site *v)
-{
- v -> refcnt += 1;
-}
-
-//push the HalfEdge into the ordered linked list of vertices
-void VoronoiDiagramGenerator::PQinsert(struct Halfedge *he,struct Site * v, double offset)
-{
- struct Halfedge *last, *next;
-
- he -> vertex = v;
- ref(v);
- he -> ystar = (double)(v -> coord.y + offset);
- last = &PQhash[PQbucket(he)];
- while ((next = last -> PQnext) != (struct Halfedge *) NULL &&
- (he -> ystar > next -> ystar ||
- (he -> ystar == next -> ystar && v -> coord.x > next->vertex->coord.x)))
- {
- last = next;
- };
- he -> PQnext = last -> PQnext;
- last -> PQnext = he;
- PQcount += 1;
-}
-
-//remove the HalfEdge from the list of vertices
-void VoronoiDiagramGenerator::PQdelete(struct Halfedge *he)
-{
- struct Halfedge *last;
-
- if(he -> vertex != (struct Site *) NULL)
- {
- last = &PQhash[PQbucket(he)];
- while (last -> PQnext != he)
- last = last -> PQnext;
-
- last -> PQnext = he -> PQnext;
- PQcount -= 1;
- deref(he -> vertex);
- he -> vertex = (struct Site *) NULL;
- };
-}
-
-int VoronoiDiagramGenerator::PQbucket(struct Halfedge *he)
-{
- int bucket;
-
- bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);
- if (bucket<0) bucket = 0;
- if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
- if (bucket < PQmin) PQmin = bucket;
- return(bucket);
-}
-
-
-
-int VoronoiDiagramGenerator::PQempty()
-{
- return(PQcount==0);
-}
-
-
-struct Point VoronoiDiagramGenerator::PQ_min()
-{
- struct Point answer;
-
- while(PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {PQmin += 1;};
- answer.x = PQhash[PQmin].PQnext -> vertex -> coord.x;
- answer.y = PQhash[PQmin].PQnext -> ystar;
- return (answer);
-}
-
-struct Halfedge * VoronoiDiagramGenerator::PQextractmin()
-{
- struct Halfedge *curr;
-
- curr = PQhash[PQmin].PQnext;
- PQhash[PQmin].PQnext = curr -> PQnext;
- PQcount -= 1;
- return(curr);
-}
-
-
-bool VoronoiDiagramGenerator::PQinitialize()
-{
- int i;
-
- PQcount = 0;
- PQmin = 0;
- PQhashsize = 4 * sqrt_nsites;
- PQhash = (struct Halfedge *) myalloc(PQhashsize * sizeof *PQhash);
-
- if(PQhash == 0)
- return false;
-
- for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (struct Halfedge *)NULL;
-
- return true;
-}
-
-
-void VoronoiDiagramGenerator::freeinit(struct Freelist *fl,int size)
-{
- fl -> head = (struct Freenode *) NULL;
- fl -> nodesize = size;
-}
-
-char * VoronoiDiagramGenerator::getfree(struct Freelist *fl)
-{
- int i;
- struct Freenode *t;
-
- if(fl->head == (struct Freenode *) NULL)
- {
- t = (struct Freenode *) myalloc(sqrt_nsites * fl->nodesize);
-
- if(t == 0)
- return 0;
-
- currentMemoryBlock->next = new FreeNodeArrayList;
- currentMemoryBlock = currentMemoryBlock->next;
- currentMemoryBlock->memory = t;
- currentMemoryBlock->next = 0;
-
- for(i=0; i<sqrt_nsites; i+=1)
- makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);
- };
- t = fl -> head;
- fl -> head = (fl -> head) -> nextfree;
- return((char *)t);
-}
-
-
-
-void VoronoiDiagramGenerator::makefree(struct Freenode *curr,struct Freelist *fl)
-{
- curr -> nextfree = fl -> head;
- fl -> head = curr;
-}
-
-void VoronoiDiagramGenerator::cleanup()
-{
- if(sites != 0)
- {
- free(sites);
- sites = 0;
- }
-
- FreeNodeArrayList* current=0, *prev = 0;
-
- current = prev = allMemoryList;
-
- while(current->next != 0)
- {
- prev = current;
- current = current->next;
- free(prev->memory);
- delete prev;
- prev = 0;
- }
-
- if(current != 0 && current->memory != 0)
- {
- free(current->memory);
- delete current;
- }
-
- allMemoryList = new FreeNodeArrayList;
- allMemoryList->next = 0;
- allMemoryList->memory = 0;
- currentMemoryBlock = allMemoryList;
-}
-
-void VoronoiDiagramGenerator::cleanupEdges()
-{
- GraphEdge* geCurrent = 0, *gePrev = 0;
- geCurrent = gePrev = allEdges;
-
- while(geCurrent != 0 && geCurrent->next != 0)
- {
- gePrev = geCurrent;
- geCurrent = geCurrent->next;
- delete gePrev;
- }
-
- allEdges = 0;
-
-}
-
-void VoronoiDiagramGenerator::cleanupEdgeList()
-{
- EdgeList* elCurrent = 0, *elPrev = 0;
- elCurrent = elPrev = allEdgeList;
-
- while (elCurrent != 0 && elCurrent->next != 0)
- {
- elPrev = elCurrent;
- elCurrent = elCurrent->next;
- delete elPrev;
- }
-
- allEdgeList = 0;
-}
-
-void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2)
-{
- GraphEdge* newEdge = new GraphEdge;
- newEdge->next = allEdges;
- allEdges = newEdge;
- newEdge->x1 = x1;
- newEdge->y1 = y1;
- newEdge->x2 = x2;
- newEdge->y2 = y2;
-}
-
-void VoronoiDiagramGenerator::pushEdgeList(Edge *e)
-{
- EdgeList* newEdge = new EdgeList;
- newEdge->next = allEdgeList;
- allEdgeList = newEdge;
- newEdge->a = e->a;
- newEdge->b = e->b;
- newEdge->c = e->c;
- if (e->ep[0]) {
- newEdge->ep0nbr = e->ep[0]->sitenbr;
- newEdge->ep0x = e->ep[0]->coord.x;
- newEdge->ep0y = e->ep[0]->coord.y;
- } else {
- newEdge->ep0nbr = -1;
- }
- if (e->ep[1]) {
- newEdge->ep1nbr = e->ep[1]->sitenbr;
- newEdge->ep1x = e->ep[1]->coord.x;
- newEdge->ep1y = e->ep[1]->coord.y;
- } else {
- newEdge->ep1nbr = -1;
- }
- newEdge->reg0nbr = e->reg[0]->sitenbr;
- newEdge->reg1nbr = e->reg[1]->sitenbr;
- newEdge->edgenbr = e->edgenbr;
-}
-
-char * VoronoiDiagramGenerator::myalloc(unsigned n)
-{
- char *t=0;
- t=(char*)malloc(n);
- total_alloc += n;
- return(t);
-}
-
-
-/* for those who don't have Cherry's plot */
-/* #include <plot.h> */
-void VoronoiDiagramGenerator::openpl(){}
-void VoronoiDiagramGenerator::line(double x1, double y1, double x2, double y2)
-{
- pushGraphEdge(x1,y1,x2,y2);
-
-}
-void VoronoiDiagramGenerator::circle(double x, double y, double radius){}
-void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){}
-
-
-
-void VoronoiDiagramGenerator::out_bisector(struct Edge *e)
-{
-
-
-}
-
-
-void VoronoiDiagramGenerator::out_ep(struct Edge *e)
-{
-
-
-}
-
-void VoronoiDiagramGenerator::out_vertex(struct Site *v)
-{
-
-}
-
-
-void VoronoiDiagramGenerator::out_site(struct Site *s)
-{
- if(!triangulate & plot & !debug)
- circle (s->coord.x, s->coord.y, cradius);
-
-}
-
-
-void VoronoiDiagramGenerator::out_triple(struct Site *s1, struct Site *s2,struct Site * s3)
-{
-
-}
-
-
-
-void VoronoiDiagramGenerator::plotinit()
-{
-// double dx,dy,d;
-//
-// dy = ymax - ymin;
-// dx = xmax - xmin;
-// d = (double)(( dx > dy ? dx : dy) * 1.1);
-// pxmin = (double)(xmin - (d-dx)/2.0);
-// pxmax = (double)(xmax + (d-dx)/2.0);
-// pymin = (double)(ymin - (d-dy)/2.0);
-// pymax = (double)(ymax + (d-dy)/2.0);
-// cradius = (double)((pxmax - pxmin)/350.0);
-// openpl();
-// range(pxmin, pymin, pxmax, pymax);
-}
-
-
-void VoronoiDiagramGenerator::clip_line(struct Edge *e)
-{
-// struct Site *s1, *s2;
-// double x1=0,x2=0,y1=0,y2=0;
-
- pushEdgeList(e);
-
-// x1 = e->reg[0]->coord.x;
-// x2 = e->reg[1]->coord.x;
-// y1 = e->reg[0]->coord.y;
-// y2 = e->reg[1]->coord.y;
-//
-// //if the distance between the two points this line was created from is less than
-// //the square root of 2, then ignore it
-// if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
-// {
-// return;
-// }
-// pxmin = borderMinX;
-// pxmax = borderMaxX;
-// pymin = borderMinY;
-// pymax = borderMaxY;
-//
-// if(e -> a == 1.0 && e ->b >= 0.0)
-// {
-// s1 = e -> ep[1];
-// s2 = e -> ep[0];
-// }
-// else
-// {
-// s1 = e -> ep[0];
-// s2 = e -> ep[1];
-// };
-//
-// if(e -> a == 1.0)
-// {
-// y1 = pymin;
-// if (s1!=(struct Site *)NULL && s1->coord.y > pymin)
-// {
-// y1 = s1->coord.y;
-// }
-// if(y1>pymax)
-// {
-// // printf("\nClipped (1) y1 = %f to %f",y1,pymax);
-// y1 = pymax;
-// //return;
-// }
-// x1 = e -> c - e -> b * y1;
-// y2 = pymax;
-// if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
-// y2 = s2->coord.y;
-//
-// if(y2<pymin)
-// {
-// //printf("\nClipped (2) y2 = %f to %f",y2,pymin);
-// y2 = pymin;
-// //return;
-// }
-// x2 = (e->c) - (e->b) * y2;
-// if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))
-// {
-// //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
-// return;
-// }
-// if(x1> pxmax)
-// { x1 = pxmax; y1 = (e -> c - x1)/e -> b;};
-// if(x1<pxmin)
-// { x1 = pxmin; y1 = (e -> c - x1)/e -> b;};
-// if(x2>pxmax)
-// { x2 = pxmax; y2 = (e -> c - x2)/e -> b;};
-// if(x2<pxmin)
-// { x2 = pxmin; y2 = (e -> c - x2)/e -> b;};
-// }
-// else
-// {
-// x1 = pxmin;
-// if (s1!=(struct Site *)NULL && s1->coord.x > pxmin)
-// x1 = s1->coord.x;
-// if(x1>pxmax)
-// {
-// //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
-// //return;
-// x1 = pxmax;
-// }
-// y1 = e -> c - e -> a * x1;
-// x2 = pxmax;
-// if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
-// x2 = s2->coord.x;
-// if(x2<pxmin)
-// {
-// //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);
-// //return;
-// x2 = pxmin;
-// }
-// y2 = e -> c - e -> a * x2;
-// if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))
-// {
-// //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
-// return;
-// }
-// if(y1> pymax)
-// { y1 = pymax; x1 = (e -> c - y1)/e -> a;};
-// if(y1<pymin)
-// { y1 = pymin; x1 = (e -> c - y1)/e -> a;};
-// if(y2>pymax)
-// { y2 = pymax; x2 = (e -> c - y2)/e -> a;};
-// if(y2<pymin)
-// { y2 = pymin; x2 = (e -> c - y2)/e -> a;};
-// };
-//
-// //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2);
-// line(x1,y1,x2,y2);
-}
-
-
-/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
-deltax, deltay (can all be estimates).
-Performance suffers if they are wrong; better to make nsites,
-deltax, and deltay too big than too small. (?) */
-
-bool VoronoiDiagramGenerator::voronoi(int triangulate)
-{
- struct Site *newsite, *bot, *top, *temp, *p;
- struct Site *v;
- struct Point newintstar;
- int pm;
- struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
- struct Edge *e;
-
- PQinitialize();
- bottomsite = nextone();
- out_site(bottomsite);
- bool retval = ELinitialize();
-
- if(!retval)
- return false;
-
- newsite = nextone();
- while(1)
- {
-
- if(!PQempty())
- newintstar = PQ_min();
-
- //if the lowest site has a smaller y value than the lowest vector intersection, process the site
- //otherwise process the vector intersection
-
- if (newsite != (struct Site *)NULL && (PQempty() || newsite -> coord.y < newintstar.y
- || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
- {/* new site is smallest - this is a site event*/
- out_site(newsite); //output the site
- lbnd = ELleftbnd(&(newsite->coord)); //get the first HalfEdge to the LEFT of the new site
- rbnd = ELright(lbnd); //get the first HalfEdge to the RIGHT of the new site
- bot = rightreg(lbnd); //if this halfedge has no edge, , bot = bottom site (whatever that is)
- e = bisect(bot, newsite); //create a new edge that bisects
- bisector = HEcreate(e, le); //create a new HalfEdge, setting its ELpm field to 0
- ELinsert(lbnd, bisector); //insert this new bisector edge between the left and right vectors in a linked list
-
- if ((p = intersect(lbnd, bisector)) != (struct Site *) NULL) //if the new bisector intersects with the left edge, remove the left edge's vertex, and put in the new one
- {
- PQdelete(lbnd);
- PQinsert(lbnd, p, dist(p,newsite));
- };
- lbnd = bisector;
- bisector = HEcreate(e, re); //create a new HalfEdge, setting its ELpm field to 1
- ELinsert(lbnd, bisector); //insert the new HE to the right of the original bisector earlier in the IF stmt
-
- if ((p = intersect(bisector, rbnd)) != (struct Site *) NULL) //if this new bisector intersects with the
- {
- PQinsert(bisector, p, dist(p,newsite)); //push the HE into the ordered linked list of vertices
- };
- newsite = nextone();
- }
- else if (!PQempty()) /* intersection is smallest - this is a vector event */
- {
- lbnd = PQextractmin(); //pop the HalfEdge with the lowest vector off the ordered list of vectors
- llbnd = ELleft(lbnd); //get the HalfEdge to the left of the above HE
- rbnd = ELright(lbnd); //get the HalfEdge to the right of the above HE
- rrbnd = ELright(rbnd); //get the HalfEdge to the right of the HE to the right of the lowest HE
- bot = leftreg(lbnd); //get the Site to the left of the left HE which it bisects
- top = rightreg(rbnd); //get the Site to the right of the right HE which it bisects
-
- out_triple(bot, top, rightreg(lbnd)); //output the triple of sites, stating that a circle goes through them
-
- v = lbnd->vertex; //get the vertex that caused this event
- makevertex(v); //set the vertex number - couldn't do this earlier since we didn't know when it would be processed
- endpoint(lbnd->ELedge,lbnd->ELpm,v); //set the endpoint of the left HalfEdge to be this vector
- endpoint(rbnd->ELedge,rbnd->ELpm,v); //set the endpoint of the right HalfEdge to be this vector
- ELdelete(lbnd); //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map
- PQdelete(rbnd); //remove all vertex events to do with the right HE
- ELdelete(rbnd); //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map
- pm = le; //set the pm variable to zero
-
- if (bot->coord.y > top->coord.y) //if the site to the left of the event is higher than the Site
- { //to the right of it, then swap them and set the 'pm' variable to 1
- temp = bot;
- bot = top;
- top = temp;
- pm = re;
- }
- e = bisect(bot, top); //create an Edge (or line) that is between the two Sites. This creates
- //the formula of the line, and assigns a line number to it
- bisector = HEcreate(e, pm); //create a HE from the Edge 'e', and make it point to that edge with its ELedge field
- ELinsert(llbnd, bisector); //insert the new bisector to the right of the left HE
- endpoint(e, re-pm, v); //set one endpoint to the new edge to be the vector point 'v'.
- //If the site to the left of this bisector is higher than the right
- //Site, then this endpoint is put in position 0; otherwise in pos 1
- deref(v); //delete the vector 'v'
-
- //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it
- if((p = intersect(llbnd, bisector)) != (struct Site *) NULL)
- {
- PQdelete(llbnd);
- PQinsert(llbnd, p, dist(p,bot));
- };
-
- //if right HE and the new bisector don't intersect, then reinsert it
- if ((p = intersect(bisector, rrbnd)) != (struct Site *) NULL)
- {
- PQinsert(bisector, p, dist(p,bot));
- };
- }
- else break;
- };
-
-
-
-
- for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
- {
- e = lbnd -> ELedge;
-
- clip_line(e);
-
- };
-
- cleanup();
-
- return true;
-
-}
-
-
-int scomp(const void *p1,const void *p2)
-{
- struct Point *s1 = (Point*)p1, *s2=(Point*)p2;
- if(s1 -> y < s2 -> y) return(-1);
- if(s1 -> y > s2 -> y) return(1);
- if(s1 -> x < s2 -> x) return(-1);
- if(s1 -> x > s2 -> x) return(1);
- return(0);
-}
-
-/* return a single in-storage site */
-struct Site * VoronoiDiagramGenerator::nextone()
-{
- struct Site *s;
- if(siteidx < nsites)
- {
- s = &sites[siteidx];
- siteidx += 1;
- return(s);
- }
- else
- return( (struct Site *)NULL);
-}
-
-bool VoronoiDiagramGenerator::getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
- int& ep1, double& ep1x, double& ep1y,
- int& reg0, int& reg1)
-{
- if (iterEdgeList == 0)
- return false;
-
- ep0 = iterEdgeList->ep0nbr;
- ep0x = iterEdgeList->ep0x;
- ep0y = iterEdgeList->ep0y;
- ep1 = iterEdgeList->ep1nbr;
- ep1x = iterEdgeList->ep1x;
- ep1y = iterEdgeList->ep1y;
- reg0 = iterEdgeList->reg0nbr;
- reg1 = iterEdgeList->reg1nbr;
-
- iterEdgeList = iterEdgeList->next;
-
- return true;
-}
-
-//PyObject* VoronoiDiagramGenerator::_getMesh()
-//{
-// PyObject *vlist, *dlist, *tlist;
-// PyObject *temp, *faces, *face;
-// int tri0, tri1, reg0, reg1;
-// double tri0x, tri0y, tri1x, tri1y;
-// int length, numtri, i;
-//
-// length = nedges;
-// numtri = nvertices;
-//
-// dlist = PyList_New(length);
-// if (!dlist) goto fail;
-// vlist = PyList_New(numtri);
-// if (!vlist) goto fail;
-// tlist = PyList_New(numtri);
-// if (!tlist) goto fail;
-//
-// for (i=0; i<numtri; i++) {
-// faces = PyList_New(0);
-// if (!faces) goto fail;
-// PyList_SET_ITEM(tlist, i, faces);
-// }
-//
-// resetEdgeListIter();
-// i = -1;
-// while (getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
-// i++;
-// face = Py_BuildValue("(ii)", reg0, reg1);
-// if (!face) goto fail;
-// PyList_SET_ITEM(dlist, i, face);
-// if (tri0 > -1) {
-// temp = PyList_GET_ITEM(vlist, tri0);
-// if (!temp) {
-// temp = Py_BuildValue("(dd)", tri0x, tri0y);
-// PyList_SET_ITEM(vlist, tri0, temp);
-// }
-// faces = PyList_GET_ITEM(tlist, tri0);
-// if (PyList_Append(faces, face) < 0) goto fail;
-// }
-// if (tri1 > -1) {
-// temp = PyList_GET_ITEM(vlist, tri1);
-// if (!temp) {
-// temp = Py_BuildValue("(dd)", tri1x, tri1y);
-// PyList_SET_ITEM(vlist, tri1, temp);
-// }
-// faces = PyList_GET_ITEM(tlist, tri1);
-// if (PyList_Append(faces, face) < 0) goto fail;
-// }
-// }
-//
-// temp = PyTuple_Pack(3, vlist, dlist, tlist);
-// if (!temp) goto fail;
-//
-// Py_DECREF(vlist);
-// Py_DECREF(dlist);
-// Py_DECREF(tlist);
-//
-// return temp;
-//
-//fail:
-// Py_XDECREF(vlist);
-// Py_XDECREF(dlist);
-// Py_XDECREF(temp);
-// Py_XDECREF(faces);
-// Py_XDECREF(face);
-// return NULL;
-//}
-
-
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/utilities/delaunay/VoronoiDiagramGenerator.h
--- a/yt/utilities/delaunay/VoronoiDiagramGenerator.h
+++ /dev/null
@@ -1,283 +0,0 @@
-/*
-* The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
-* Bell Laboratories.
-* Permission to use, copy, modify, and distribute this software for any
-* purpose without fee is hereby granted, provided that this entire notice
-* is included in all copies of any software which is or includes a copy
-* or modification of this software and in all copies of the supporting
-* documentation for such software.
-* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
-* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
-* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
-* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
-*/
-
-/*
-* This code was originally written by Stephan Fortune in C code. I, Shane O'Sullivan,
-* have since modified it, encapsulating it in a C++ class and, fixing memory leaks and
-* adding accessors to the Voronoi Edges.
-* Permission to use, copy, modify, and distribute this software for any
-* purpose without fee is hereby granted, provided that this entire notice
-* is included in all copies of any software which is or includes a copy
-* or modification of this software and in all copies of the supporting
-* documentation for such software.
-* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
-* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
-* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
-* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
-*/
-
-#ifndef VORONOI_DIAGRAM_GENERATOR
-#define VORONOI_DIAGRAM_GENERATOR
-
-#include "Python.h"
-#include "numpy/arrayobject.h"
-
-#include <math.h>
-#include <stdlib.h>
-#include <string.h>
-#include <iostream>
-
-
-#ifndef NULL
-#define NULL 0
-#endif
-#define DELETED -2
-
-#define le 0
-#define re 1
-
-#ifndef MAX
-#define MAX(x, y) (x > y ? x: y)
-#endif
-
-struct Freenode
-{
- struct Freenode *nextfree;
-};
-
-struct FreeNodeArrayList
-{
- struct Freenode* memory;
- struct FreeNodeArrayList* next;
-
-};
-
-struct Freelist
-{
- struct Freenode *head;
- int nodesize;
-};
-
-struct Point
-{
- double x,y;
-};
-
-// structure used both for sites and for vertices
-struct Site
-{
- struct Point coord;
- int sitenbr;
- int refcnt;
-};
-
-
-
-struct Edge
-{
- double a,b,c;
- struct Site *ep[2];
- struct Site *reg[2];
- int edgenbr;
-};
-
-struct EdgeList
-{
- double a,b,c;
- int ep0nbr;
- double ep0x, ep0y;
- int ep1nbr;
- double ep1x, ep1y;
- int reg0nbr;
- int reg1nbr;
- int edgenbr;
- struct EdgeList *next;
-};
-
-struct GraphEdge
-{
- double x1,y1,x2,y2;
- struct GraphEdge* next;
-};
-
-
-
-
-struct Halfedge
-{
- struct Halfedge *ELleft, *ELright;
- struct Edge *ELedge;
- int ELrefcnt;
- char ELpm;
- struct Site *vertex;
- double ystar;
- struct Halfedge *PQnext;
-};
-
-
-
-
-class VoronoiDiagramGenerator
-{
-public:
- VoronoiDiagramGenerator();
- ~VoronoiDiagramGenerator();
-
- bool generateVoronoi(double *xValues, double *yValues, int numPoints, double minX, double maxX, double minY, double maxY, double minDist=0);
-
- void resetIterator()
- {
- iteratorEdges = allEdges;
- }
-
- bool getNext(double& x1, double& y1, double& x2, double& y2)
- {
- if(iteratorEdges == 0)
- return false;
-
- x1 = iteratorEdges->x1;
- x2 = iteratorEdges->x2;
- y1 = iteratorEdges->y1;
- y2 = iteratorEdges->y2;
-
- iteratorEdges = iteratorEdges->next;
-
- return true;
- }
-
- void resetEdgeListIter()
- {
- iterEdgeList = allEdgeList;
- }
-
- bool getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
- int& ep1, double& ep1x, double& ep1y,
- int& reg0, int& reg1);
-
- void getNumbers(int& edges, int& vertices) {
- edges = nedges;
- vertices = nvertices;
- }
-
-private:
- void cleanup();
- void cleanupEdgeList();
- void cleanupEdges();
- char *getfree(struct Freelist *fl);
- struct Halfedge *PQfind();
- int PQempty();
-
-
- struct Halfedge **ELhash;
- struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
- struct Halfedge *HEcreate(struct Edge *e,int pm);
-
-
- struct Point PQ_min();
- struct Halfedge *PQextractmin();
- void freeinit(struct Freelist *fl,int size);
- void makefree(struct Freenode *curr,struct Freelist *fl);
- void geominit();
- void plotinit();
- bool voronoi(int triangulate);
- void ref(struct Site *v);
- void deref(struct Site *v);
- void endpoint(struct Edge *e,int lr,struct Site * s);
-
- void ELdelete(struct Halfedge *he);
- struct Halfedge *ELleftbnd(struct Point *p);
- struct Halfedge *ELright(struct Halfedge *he);
- void makevertex(struct Site *v);
- void out_triple(struct Site *s1, struct Site *s2,struct Site * s3);
-
- void PQinsert(struct Halfedge *he,struct Site * v, double offset);
- void PQdelete(struct Halfedge *he);
- bool ELinitialize();
- void ELinsert(struct Halfedge *lb, struct Halfedge *newHe);
- struct Halfedge *ELgethash(int b);
- struct Halfedge *ELleft(struct Halfedge *he);
- struct Site *leftreg(struct Halfedge *he);
- void out_site(struct Site *s);
- bool PQinitialize();
- int PQbucket(struct Halfedge *he);
- void clip_line(struct Edge *e);
- char *myalloc(unsigned n);
- int right_of(struct Halfedge *el,struct Point *p);
-
- struct Site *rightreg(struct Halfedge *he);
- struct Edge *bisect(struct Site *s1, struct Site *s2);
- double dist(struct Site *s,struct Site *t);
- struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p=0);
-
- void out_bisector(struct Edge *e);
- void out_ep(struct Edge *e);
- void out_vertex(struct Site *v);
- struct Site *nextone();
-
- void pushGraphEdge(double x1, double y1, double x2, double y2);
- void pushEdgeList(Edge *e);
-
- void openpl();
- void line(double x1, double y1, double x2, double y2);
- void circle(double x, double y, double radius);
- void range(double minX, double minY, double maxX, double maxY);
-
-
- struct Freelist hfl;
- struct Halfedge *ELleftend, *ELrightend;
- int ELhashsize;
-
- int triangulate, sorted, plot, debug;
- double xmin, xmax, ymin, ymax, deltax, deltay;
-
- struct Site *sites;
- int nsites;
- int siteidx;
- int sqrt_nsites;
- int nvertices;
- struct Freelist sfl;
- struct Site *bottomsite;
-
- int nedges;
- struct Freelist efl;
- int PQhashsize;
- struct Halfedge *PQhash;
- int PQcount;
- int PQmin;
-
- int ntry, totalsearch;
- double pxmin, pxmax, pymin, pymax, cradius;
- int total_alloc;
-
- double borderMinX, borderMaxX, borderMinY, borderMaxY;
-
- FreeNodeArrayList* allMemoryList;
- FreeNodeArrayList* currentMemoryBlock;
-
- GraphEdge* allEdges;
- GraphEdge* iteratorEdges;
-
- EdgeList* allEdgeList;
- EdgeList* iterEdgeList;
-
- double minDistanceBetweenSites;
-
-};
-
-int scomp(const void *p1, const void *p2);
-
-
-#endif
-
-
diff -r 99a037526e5d1b5d8ec29b24146344d33f6f154e -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd yt/utilities/delaunay/__init__.py
--- a/yt/utilities/delaunay/__init__.py
+++ /dev/null
@@ -1,10 +0,0 @@
-"""Delaunay triangulation and interpolation tools.
-
-:Author: Robert Kern <robert.kern at gmail.com>
-:Copyright: Copyright 2005 Robert Kern.
-:License: BSD-style license. See LICENSE.txt
-"""
-
-from _delaunay import delaunay
-from triangulate import *
-from interpolate import *
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/8960f029b62f/
Changeset: 8960f029b62f
Branch: yt
User: xarthisius
Date: 2013-09-06 08:44:01
Summary: Move progressbar to extern
Affected #: 4 files
diff -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd -r 8960f029b62f4c703889a784449014ddaf7e0a40 yt/extern/progressbar.py
--- /dev/null
+++ b/yt/extern/progressbar.py
@@ -0,0 +1,368 @@
+#!/usr/bin/python
+# -*- coding: iso-8859-1 -*-
+#
+# progressbar - Text progressbar library for python.
+# Copyright (c) 2005 Nilton Volpato
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+"""Text progressbar library for python.
+
+This library provides a text mode progressbar. This is tipically used
+to display the progress of a long running operation, providing a
+visual clue that processing is underway.
+
+The ProgressBar class manages the progress, and the format of the line
+is given by a number of widgets. A widget is an object that may
+display diferently depending on the state of the progress. There are
+three types of widget:
+- a string, which always shows itself;
+- a ProgressBarWidget, which may return a diferent value every time
+it's update method is called; and
+- a ProgressBarWidgetHFill, which is like ProgressBarWidget, except it
+expands to fill the remaining width of the line.
+
+The progressbar module is very easy to use, yet very powerful. And
+automatically supports features like auto-resizing when available.
+"""
+
+__author__ = "Nilton Volpato"
+__author_email__ = "first-name dot last-name @ gmail.com"
+__date__ = "2006-05-07"
+__version__ = "2.2"
+
+# Changelog
+#
+# 2006-05-07: v2.2 fixed bug in windows
+# 2005-12-04: v2.1 autodetect terminal width, added start method
+# 2005-12-04: v2.0 everything is now a widget (wow!)
+# 2005-12-03: v1.0 rewrite using widgets
+# 2005-06-02: v0.5 rewrite
+# 2004-??-??: v0.1 first version
+
+
+import sys, time
+from array import array
+try:
+ from fcntl import ioctl
+ import termios
+except ImportError:
+ pass
+import signal
+
+class ProgressBarWidget(object):
+ """This is an element of ProgressBar formatting.
+
+ The ProgressBar object will call it's update value when an update
+ is needed. It's size may change between call, but the results will
+ not be good if the size changes drastically and repeatedly.
+ """
+ def update(self, pbar):
+ """Returns the string representing the widget.
+
+ The parameter pbar is a reference to the calling ProgressBar,
+ where one can access attributes of the class for knowing how
+ the update must be made.
+
+ At least this function must be overriden."""
+ pass
+
+class ProgressBarWidgetHFill(object):
+ """This is a variable width element of ProgressBar formatting.
+
+ The ProgressBar object will call it's update value, informing the
+ width this object must the made. This is like TeX \\hfill, it will
+ expand to fill the line. You can use more than one in the same
+ line, and they will all have the same width, and together will
+ fill the line.
+ """
+ def update(self, pbar, width):
+ """Returns the string representing the widget.
+
+ The parameter pbar is a reference to the calling ProgressBar,
+ where one can access attributes of the class for knowing how
+ the update must be made. The parameter width is the total
+ horizontal width the widget must have.
+
+ At least this function must be overriden."""
+ pass
+
+
+class ETA(ProgressBarWidget):
+ "Widget for the Estimated Time of Arrival"
+ def format_time(self, seconds):
+ return time.strftime('%H:%M:%S', time.gmtime(seconds))
+ def update(self, pbar):
+ if pbar.currval == 0:
+ return 'ETA: --:--:--'
+ elif pbar.finished:
+ return 'Time: %s' % self.format_time(pbar.seconds_elapsed)
+ else:
+ elapsed = pbar.seconds_elapsed
+ eta = elapsed * pbar.maxval / pbar.currval - elapsed
+ return 'ETA: %s' % self.format_time(eta)
+
+class FileTransferSpeed(ProgressBarWidget):
+ "Widget for showing the transfer speed (useful for file transfers)."
+ def __init__(self):
+ self.fmt = '%6.2f %s'
+ self.units = ['B','K','M','G','T','P']
+ def update(self, pbar):
+ if pbar.seconds_elapsed < 2e-6:#== 0:
+ bps = 0.0
+ else:
+ bps = float(pbar.currval) / pbar.seconds_elapsed
+ spd = bps
+ for u in self.units:
+ if spd < 1000:
+ break
+ spd /= 1000
+ return self.fmt % (spd, u+'/s')
+
+class RotatingMarker(ProgressBarWidget):
+ "A rotating marker for filling the bar of progress."
+ def __init__(self, markers='|/-\\'):
+ self.markers = markers
+ self.curmark = -1
+ def update(self, pbar):
+ if pbar.finished:
+ return self.markers[0]
+ self.curmark = (self.curmark + 1)%len(self.markers)
+ return self.markers[self.curmark]
+
+class Percentage(ProgressBarWidget):
+ "Just the percentage done."
+ def update(self, pbar):
+ return '%3d%%' % pbar.percentage()
+
+class Bar(ProgressBarWidgetHFill):
+ "The bar of progress. It will strech to fill the line."
+ def __init__(self, marker='#', left='|', right='|'):
+ self.marker = marker
+ self.left = left
+ self.right = right
+ def _format_marker(self, pbar):
+ if isinstance(self.marker, (str, unicode)):
+ return self.marker
+ else:
+ return self.marker.update(pbar)
+ def update(self, pbar, width):
+ percent = pbar.percentage()
+ cwidth = width - len(self.left) - len(self.right)
+ marked_width = int(percent * cwidth / 100)
+ m = self._format_marker(pbar)
+ bar = (self.left + (m*marked_width).ljust(cwidth) + self.right)
+ return bar
+
+class ReverseBar(Bar):
+ "The reverse bar of progress, or bar of regress. :)"
+ def update(self, pbar, width):
+ percent = pbar.percentage()
+ cwidth = width - len(self.left) - len(self.right)
+ marked_width = int(percent * cwidth / 100)
+ m = self._format_marker(pbar)
+ bar = (self.left + (m*marked_width).rjust(cwidth) + self.right)
+ return bar
+
+default_widgets = [Percentage(), ' ', Bar()]
+class ProgressBar(object):
+ """This is the ProgressBar class, it updates and prints the bar.
+
+ The term_width parameter may be an integer. Or None, in which case
+ it will try to guess it, if it fails it will default to 80 columns.
+
+ The simple use is like this:
+ >>> pbar = ProgressBar().start()
+ >>> for i in xrange(100):
+ ... # do something
+ ... pbar.update(i+1)
+ ...
+ >>> pbar.finish()
+
+ But anything you want to do is possible (well, almost anything).
+ You can supply different widgets of any type in any order. And you
+ can even write your own widgets! There are many widgets already
+ shipped and you should experiment with them.
+
+ When implementing a widget update method you may access any
+ attribute or function of the ProgressBar object calling the
+ widget's update method. The most important attributes you would
+ like to access are:
+ - currval: current value of the progress, 0 <= currval <= maxval
+ - maxval: maximum (and final) value of the progress
+ - finished: True if the bar is have finished (reached 100%), False o/w
+ - start_time: first time update() method of ProgressBar was called
+ - seconds_elapsed: seconds elapsed since start_time
+ - percentage(): percentage of the progress (this is a method)
+ """
+ def __init__(self, maxval=100, widgets=default_widgets, term_width=None,
+ fd=sys.stderr):
+ assert maxval > 0
+ self.maxval = maxval
+ self.widgets = widgets
+ self.fd = fd
+ self.signal_set = False
+ if term_width is None:
+ try:
+ self.handle_resize(None,None)
+ signal.signal(signal.SIGWINCH, self.handle_resize)
+ self.signal_set = True
+ except:
+ self.term_width = 79
+ else:
+ self.term_width = term_width
+
+ self.currval = 0
+ self.finished = False
+ self.prev_percentage = -1
+ self.start_time = None
+ self.seconds_elapsed = 0
+
+ def handle_resize(self, signum, frame):
+ h,w=array('h', ioctl(self.fd,termios.TIOCGWINSZ,'\0'*8))[:2]
+ self.term_width = w
+
+ def percentage(self):
+ "Returns the percentage of the progress."
+ return self.currval*100.0 / self.maxval
+
+ def _format_widgets(self):
+ r = []
+ hfill_inds = []
+ num_hfill = 0
+ currwidth = 0
+ for i, w in enumerate(self.widgets):
+ if isinstance(w, ProgressBarWidgetHFill):
+ r.append(w)
+ hfill_inds.append(i)
+ num_hfill += 1
+ elif isinstance(w, (str, unicode)):
+ r.append(w)
+ currwidth += len(w)
+ else:
+ weval = w.update(self)
+ currwidth += len(weval)
+ r.append(weval)
+ for iw in hfill_inds:
+ r[iw] = r[iw].update(self, (self.term_width-currwidth)/num_hfill)
+ return r
+
+ def _format_line(self):
+ return ''.join(self._format_widgets()).ljust(self.term_width)
+
+ def _need_update(self):
+ return int(self.percentage()) != int(self.prev_percentage)
+
+ def update(self, value):
+ "Updates the progress bar to a new value."
+ assert 0 <= value <= self.maxval
+ self.currval = value
+ if not self._need_update() or self.finished:
+ return
+ if not self.start_time:
+ self.start_time = time.time()
+ self.seconds_elapsed = time.time() - self.start_time
+ self.prev_percentage = self.percentage()
+ if value != self.maxval:
+ self.fd.write(self._format_line() + '\r')
+ else:
+ self.finished = True
+ self.fd.write(self._format_line() + '\n')
+
+ def start(self):
+ """Start measuring time, and prints the bar at 0%.
+
+ It returns self so you can use it like this:
+ >>> pbar = ProgressBar().start()
+ >>> for i in xrange(100):
+ ... # do something
+ ... pbar.update(i+1)
+ ...
+ >>> pbar.finish()
+ """
+ self.update(0)
+ return self
+
+ def finish(self):
+ """Used to tell the progress is finished."""
+ self.update(self.maxval)
+ if self.signal_set:
+ signal.signal(signal.SIGWINCH, signal.SIG_DFL)
+
+
+
+
+
+
+if __name__=='__main__':
+ import os
+
+ def example1():
+ widgets = ['Test: ', Percentage(), ' ', Bar(marker=RotatingMarker()),
+ ' ', ETA(), ' ', FileTransferSpeed()]
+ pbar = ProgressBar(widgets=widgets, maxval=10000000).start()
+ for i in range(1000000):
+ # do something
+ pbar.update(10*i+1)
+ pbar.finish()
+ print
+
+ def example2():
+ class CrazyFileTransferSpeed(FileTransferSpeed):
+ "It's bigger between 45 and 80 percent"
+ def update(self, pbar):
+ if 45 < pbar.percentage() < 80:
+ return 'Bigger Now ' + FileTransferSpeed.update(self,pbar)
+ else:
+ return FileTransferSpeed.update(self,pbar)
+
+ widgets = [CrazyFileTransferSpeed(),' <<<', Bar(), '>>> ', Percentage(),' ', ETA()]
+ pbar = ProgressBar(widgets=widgets, maxval=10000000)
+ # maybe do something
+ pbar.start()
+ for i in range(2000000):
+ # do something
+ pbar.update(5*i+1)
+ pbar.finish()
+ print
+
+ def example3():
+ widgets = [Bar('>'), ' ', ETA(), ' ', ReverseBar('<')]
+ pbar = ProgressBar(widgets=widgets, maxval=10000000).start()
+ for i in range(1000000):
+ # do something
+ pbar.update(10*i+1)
+ pbar.finish()
+ print
+
+ def example4():
+ widgets = ['Test: ', Percentage(), ' ',
+ Bar(marker='0',left='[',right=']'),
+ ' ', ETA(), ' ', FileTransferSpeed()]
+ pbar = ProgressBar(widgets=widgets, maxval=500)
+ pbar.start()
+ for i in range(100,500+1,50):
+ time.sleep(0.2)
+ pbar.update(i)
+ pbar.finish()
+ print
+
+
+ example1()
+ example2()
+ example3()
+ example4()
+
diff -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd -r 8960f029b62f4c703889a784449014ddaf7e0a40 yt/funcs.py
--- a/yt/funcs.py
+++ b/yt/funcs.py
@@ -34,7 +34,7 @@
from yt.utilities.exceptions import *
from yt.utilities.logger import ytLogger as mylog
from yt.utilities.definitions import inv_axis_names, axis_names, x_dict, y_dict
-import yt.utilities.progressbar as pb
+import yt.extern.progressbar as pb
import yt.utilities.rpdb as rpdb
from collections import defaultdict
from functools import wraps
diff -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd -r 8960f029b62f4c703889a784449014ddaf7e0a40 yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -46,7 +46,7 @@
import matplotlib.image as mpimg
import yt.visualization.plot_window as pw
-import yt.utilities.progressbar as progressbar
+import yt.extern.progressbar as progressbar
mylog = logging.getLogger('nose.plugins.answer-testing')
run_big_data = False
diff -r dd8cfa5ab087dc72b01a5c48537fe8b3e89c20bd -r 8960f029b62f4c703889a784449014ddaf7e0a40 yt/utilities/progressbar.py
--- a/yt/utilities/progressbar.py
+++ /dev/null
@@ -1,368 +0,0 @@
-#!/usr/bin/python
-# -*- coding: iso-8859-1 -*-
-#
-# progressbar - Text progressbar library for python.
-# Copyright (c) 2005 Nilton Volpato
-#
-# This library is free software; you can redistribute it and/or
-# modify it under the terms of the GNU Lesser General Public
-# License as published by the Free Software Foundation; either
-# version 2.1 of the License, or (at your option) any later version.
-#
-# This library is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-# Lesser General Public License for more details.
-#
-# You should have received a copy of the GNU Lesser General Public
-# License along with this library; if not, write to the Free Software
-# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-
-"""Text progressbar library for python.
-
-This library provides a text mode progressbar. This is tipically used
-to display the progress of a long running operation, providing a
-visual clue that processing is underway.
-
-The ProgressBar class manages the progress, and the format of the line
-is given by a number of widgets. A widget is an object that may
-display diferently depending on the state of the progress. There are
-three types of widget:
-- a string, which always shows itself;
-- a ProgressBarWidget, which may return a diferent value every time
-it's update method is called; and
-- a ProgressBarWidgetHFill, which is like ProgressBarWidget, except it
-expands to fill the remaining width of the line.
-
-The progressbar module is very easy to use, yet very powerful. And
-automatically supports features like auto-resizing when available.
-"""
-
-__author__ = "Nilton Volpato"
-__author_email__ = "first-name dot last-name @ gmail.com"
-__date__ = "2006-05-07"
-__version__ = "2.2"
-
-# Changelog
-#
-# 2006-05-07: v2.2 fixed bug in windows
-# 2005-12-04: v2.1 autodetect terminal width, added start method
-# 2005-12-04: v2.0 everything is now a widget (wow!)
-# 2005-12-03: v1.0 rewrite using widgets
-# 2005-06-02: v0.5 rewrite
-# 2004-??-??: v0.1 first version
-
-
-import sys, time
-from array import array
-try:
- from fcntl import ioctl
- import termios
-except ImportError:
- pass
-import signal
-
-class ProgressBarWidget(object):
- """This is an element of ProgressBar formatting.
-
- The ProgressBar object will call it's update value when an update
- is needed. It's size may change between call, but the results will
- not be good if the size changes drastically and repeatedly.
- """
- def update(self, pbar):
- """Returns the string representing the widget.
-
- The parameter pbar is a reference to the calling ProgressBar,
- where one can access attributes of the class for knowing how
- the update must be made.
-
- At least this function must be overriden."""
- pass
-
-class ProgressBarWidgetHFill(object):
- """This is a variable width element of ProgressBar formatting.
-
- The ProgressBar object will call it's update value, informing the
- width this object must the made. This is like TeX \\hfill, it will
- expand to fill the line. You can use more than one in the same
- line, and they will all have the same width, and together will
- fill the line.
- """
- def update(self, pbar, width):
- """Returns the string representing the widget.
-
- The parameter pbar is a reference to the calling ProgressBar,
- where one can access attributes of the class for knowing how
- the update must be made. The parameter width is the total
- horizontal width the widget must have.
-
- At least this function must be overriden."""
- pass
-
-
-class ETA(ProgressBarWidget):
- "Widget for the Estimated Time of Arrival"
- def format_time(self, seconds):
- return time.strftime('%H:%M:%S', time.gmtime(seconds))
- def update(self, pbar):
- if pbar.currval == 0:
- return 'ETA: --:--:--'
- elif pbar.finished:
- return 'Time: %s' % self.format_time(pbar.seconds_elapsed)
- else:
- elapsed = pbar.seconds_elapsed
- eta = elapsed * pbar.maxval / pbar.currval - elapsed
- return 'ETA: %s' % self.format_time(eta)
-
-class FileTransferSpeed(ProgressBarWidget):
- "Widget for showing the transfer speed (useful for file transfers)."
- def __init__(self):
- self.fmt = '%6.2f %s'
- self.units = ['B','K','M','G','T','P']
- def update(self, pbar):
- if pbar.seconds_elapsed < 2e-6:#== 0:
- bps = 0.0
- else:
- bps = float(pbar.currval) / pbar.seconds_elapsed
- spd = bps
- for u in self.units:
- if spd < 1000:
- break
- spd /= 1000
- return self.fmt % (spd, u+'/s')
-
-class RotatingMarker(ProgressBarWidget):
- "A rotating marker for filling the bar of progress."
- def __init__(self, markers='|/-\\'):
- self.markers = markers
- self.curmark = -1
- def update(self, pbar):
- if pbar.finished:
- return self.markers[0]
- self.curmark = (self.curmark + 1)%len(self.markers)
- return self.markers[self.curmark]
-
-class Percentage(ProgressBarWidget):
- "Just the percentage done."
- def update(self, pbar):
- return '%3d%%' % pbar.percentage()
-
-class Bar(ProgressBarWidgetHFill):
- "The bar of progress. It will strech to fill the line."
- def __init__(self, marker='#', left='|', right='|'):
- self.marker = marker
- self.left = left
- self.right = right
- def _format_marker(self, pbar):
- if isinstance(self.marker, (str, unicode)):
- return self.marker
- else:
- return self.marker.update(pbar)
- def update(self, pbar, width):
- percent = pbar.percentage()
- cwidth = width - len(self.left) - len(self.right)
- marked_width = int(percent * cwidth / 100)
- m = self._format_marker(pbar)
- bar = (self.left + (m*marked_width).ljust(cwidth) + self.right)
- return bar
-
-class ReverseBar(Bar):
- "The reverse bar of progress, or bar of regress. :)"
- def update(self, pbar, width):
- percent = pbar.percentage()
- cwidth = width - len(self.left) - len(self.right)
- marked_width = int(percent * cwidth / 100)
- m = self._format_marker(pbar)
- bar = (self.left + (m*marked_width).rjust(cwidth) + self.right)
- return bar
-
-default_widgets = [Percentage(), ' ', Bar()]
-class ProgressBar(object):
- """This is the ProgressBar class, it updates and prints the bar.
-
- The term_width parameter may be an integer. Or None, in which case
- it will try to guess it, if it fails it will default to 80 columns.
-
- The simple use is like this:
- >>> pbar = ProgressBar().start()
- >>> for i in xrange(100):
- ... # do something
- ... pbar.update(i+1)
- ...
- >>> pbar.finish()
-
- But anything you want to do is possible (well, almost anything).
- You can supply different widgets of any type in any order. And you
- can even write your own widgets! There are many widgets already
- shipped and you should experiment with them.
-
- When implementing a widget update method you may access any
- attribute or function of the ProgressBar object calling the
- widget's update method. The most important attributes you would
- like to access are:
- - currval: current value of the progress, 0 <= currval <= maxval
- - maxval: maximum (and final) value of the progress
- - finished: True if the bar is have finished (reached 100%), False o/w
- - start_time: first time update() method of ProgressBar was called
- - seconds_elapsed: seconds elapsed since start_time
- - percentage(): percentage of the progress (this is a method)
- """
- def __init__(self, maxval=100, widgets=default_widgets, term_width=None,
- fd=sys.stderr):
- assert maxval > 0
- self.maxval = maxval
- self.widgets = widgets
- self.fd = fd
- self.signal_set = False
- if term_width is None:
- try:
- self.handle_resize(None,None)
- signal.signal(signal.SIGWINCH, self.handle_resize)
- self.signal_set = True
- except:
- self.term_width = 79
- else:
- self.term_width = term_width
-
- self.currval = 0
- self.finished = False
- self.prev_percentage = -1
- self.start_time = None
- self.seconds_elapsed = 0
-
- def handle_resize(self, signum, frame):
- h,w=array('h', ioctl(self.fd,termios.TIOCGWINSZ,'\0'*8))[:2]
- self.term_width = w
-
- def percentage(self):
- "Returns the percentage of the progress."
- return self.currval*100.0 / self.maxval
-
- def _format_widgets(self):
- r = []
- hfill_inds = []
- num_hfill = 0
- currwidth = 0
- for i, w in enumerate(self.widgets):
- if isinstance(w, ProgressBarWidgetHFill):
- r.append(w)
- hfill_inds.append(i)
- num_hfill += 1
- elif isinstance(w, (str, unicode)):
- r.append(w)
- currwidth += len(w)
- else:
- weval = w.update(self)
- currwidth += len(weval)
- r.append(weval)
- for iw in hfill_inds:
- r[iw] = r[iw].update(self, (self.term_width-currwidth)/num_hfill)
- return r
-
- def _format_line(self):
- return ''.join(self._format_widgets()).ljust(self.term_width)
-
- def _need_update(self):
- return int(self.percentage()) != int(self.prev_percentage)
-
- def update(self, value):
- "Updates the progress bar to a new value."
- assert 0 <= value <= self.maxval
- self.currval = value
- if not self._need_update() or self.finished:
- return
- if not self.start_time:
- self.start_time = time.time()
- self.seconds_elapsed = time.time() - self.start_time
- self.prev_percentage = self.percentage()
- if value != self.maxval:
- self.fd.write(self._format_line() + '\r')
- else:
- self.finished = True
- self.fd.write(self._format_line() + '\n')
-
- def start(self):
- """Start measuring time, and prints the bar at 0%.
-
- It returns self so you can use it like this:
- >>> pbar = ProgressBar().start()
- >>> for i in xrange(100):
- ... # do something
- ... pbar.update(i+1)
- ...
- >>> pbar.finish()
- """
- self.update(0)
- return self
-
- def finish(self):
- """Used to tell the progress is finished."""
- self.update(self.maxval)
- if self.signal_set:
- signal.signal(signal.SIGWINCH, signal.SIG_DFL)
-
-
-
-
-
-
-if __name__=='__main__':
- import os
-
- def example1():
- widgets = ['Test: ', Percentage(), ' ', Bar(marker=RotatingMarker()),
- ' ', ETA(), ' ', FileTransferSpeed()]
- pbar = ProgressBar(widgets=widgets, maxval=10000000).start()
- for i in range(1000000):
- # do something
- pbar.update(10*i+1)
- pbar.finish()
- print
-
- def example2():
- class CrazyFileTransferSpeed(FileTransferSpeed):
- "It's bigger between 45 and 80 percent"
- def update(self, pbar):
- if 45 < pbar.percentage() < 80:
- return 'Bigger Now ' + FileTransferSpeed.update(self,pbar)
- else:
- return FileTransferSpeed.update(self,pbar)
-
- widgets = [CrazyFileTransferSpeed(),' <<<', Bar(), '>>> ', Percentage(),' ', ETA()]
- pbar = ProgressBar(widgets=widgets, maxval=10000000)
- # maybe do something
- pbar.start()
- for i in range(2000000):
- # do something
- pbar.update(5*i+1)
- pbar.finish()
- print
-
- def example3():
- widgets = [Bar('>'), ' ', ETA(), ' ', ReverseBar('<')]
- pbar = ProgressBar(widgets=widgets, maxval=10000000).start()
- for i in range(1000000):
- # do something
- pbar.update(10*i+1)
- pbar.finish()
- print
-
- def example4():
- widgets = ['Test: ', Percentage(), ' ',
- Bar(marker='0',left='[',right=']'),
- ' ', ETA(), ' ', FileTransferSpeed()]
- pbar = ProgressBar(widgets=widgets, maxval=500)
- pbar.start()
- for i in range(100,500+1,50):
- time.sleep(0.2)
- pbar.update(i)
- pbar.finish()
- print
-
-
- example1()
- example2()
- example3()
- example4()
-
https://bitbucket.org/yt_analysis/yt/commits/6938c5757a92/
Changeset: 6938c5757a92
Branch: yt
User: xarthisius
Date: 2013-09-06 08:47:41
Summary: Move pykdtree to extern
Affected #: 3 files
diff -r 8960f029b62f4c703889a784449014ddaf7e0a40 -r 6938c5757a9246082a90a7d6e0678d052a1e7489 yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
+++ b/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
@@ -54,7 +54,7 @@
import glob
from yt.funcs import *
-from yt.utilities.pykdtree import KDTree
+from yt.extern.pykdtree import KDTree
import yt.extern.pydot as pydot
# We don't currently use this, but we may again find a use for it in the
diff -r 8960f029b62f4c703889a784449014ddaf7e0a40 -r 6938c5757a9246082a90a7d6e0678d052a1e7489 yt/extern/pykdtree.py
--- /dev/null
+++ b/yt/extern/pykdtree.py
@@ -0,0 +1,823 @@
+# Copyright Anne M. Archibald 2008
+# Released under the scipy license
+import sys
+import numpy as np
+from heapq import heappush, heappop
+
+def minkowski_distance_p(x,y,p=2):
+ """Compute the pth power of the L**p distance between x and y
+
+ For efficiency, this function computes the L**p distance but does
+ not extract the pth root. If p is 1 or infinity, this is equal to
+ the actual L**p distance.
+ """
+ x = np.asarray(x)
+ y = np.asarray(y)
+ if p==np.inf:
+ return np.amax(np.abs(y-x),axis=-1)
+ elif p==1:
+ return np.sum(np.abs(y-x),axis=-1)
+ else:
+ return np.sum(np.abs(y-x)**p,axis=-1)
+def minkowski_distance(x,y,p=2):
+ """Compute the L**p distance between x and y"""
+ x = np.asarray(x)
+ y = np.asarray(y)
+ if p==np.inf or p==1:
+ return minkowski_distance_p(x,y,p)
+ else:
+ return minkowski_distance_p(x,y,p)**(1./p)
+
+class Rectangle(object):
+ """Hyperrectangle class.
+
+ Represents a Cartesian product of intervals.
+ """
+ def __init__(self, maxes, mins):
+ """Construct a hyperrectangle."""
+ self.maxes = np.maximum(maxes,mins).astype(np.float)
+ self.mins = np.minimum(maxes,mins).astype(np.float)
+ self.m, = self.maxes.shape
+
+ def __repr__(self):
+ return "<Rectangle %s>" % zip(self.mins, self.maxes)
+
+ def volume(self):
+ """Total volume."""
+ return np.prod(self.maxes-self.mins)
+
+ def split(self, d, split):
+ """Produce two hyperrectangles by splitting along axis d.
+
+ In general, if you need to compute maximum and minimum
+ distances to the children, it can be done more efficiently
+ by updating the maximum and minimum distances to the parent.
+ """ # FIXME: do this
+ mid = np.copy(self.maxes)
+ mid[d] = split
+ less = Rectangle(self.mins, mid)
+ mid = np.copy(self.mins)
+ mid[d] = split
+ greater = Rectangle(mid, self.maxes)
+ return less, greater
+
+ def min_distance_point(self, x, p=2.):
+ """Compute the minimum distance between x and a point in the hyperrectangle."""
+ return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-x,x-self.maxes)),p)
+
+ def max_distance_point(self, x, p=2.):
+ """Compute the maximum distance between x and a point in the hyperrectangle."""
+ return minkowski_distance(0, np.maximum(self.maxes-x,x-self.mins),p)
+
+ def min_distance_rectangle(self, other, p=2.):
+ """Compute the minimum distance between points in the two hyperrectangles."""
+ return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-other.maxes,other.mins-self.maxes)),p)
+
+ def max_distance_rectangle(self, other, p=2.):
+ """Compute the maximum distance between points in the two hyperrectangles."""
+ return minkowski_distance(0, np.maximum(self.maxes-other.mins,other.maxes-self.mins),p)
+
+
+class KDTree(object):
+ """
+ kd-tree for quick nearest-neighbor lookup
+
+ This class provides an index into a set of k-dimensional points
+ which can be used to rapidly look up the nearest neighbors of any
+ point.
+
+ The algorithm used is described in Maneewongvatana and Mount 1999.
+ The general idea is that the kd-tree is a binary tree, each of whose
+ nodes represents an axis-aligned hyperrectangle. Each node specifies
+ an axis and splits the set of points based on whether their coordinate
+ along that axis is greater than or less than a particular value.
+
+ During construction, the axis and splitting point are chosen by the
+ "sliding midpoint" rule, which ensures that the cells do not all
+ become long and thin.
+
+ The tree can be queried for the r closest neighbors of any given point
+ (optionally returning only those within some maximum distance of the
+ point). It can also be queried, with a substantial gain in efficiency,
+ for the r approximate closest neighbors.
+
+ For large dimensions (20 is already large) do not expect this to run
+ significantly faster than brute force. High-dimensional nearest-neighbor
+ queries are a substantial open problem in computer science.
+
+ The tree also supports all-neighbors queries, both with arrays of points
+ and with other kd-trees. These do use a reasonably efficient algorithm,
+ but the kd-tree is not necessarily the best data structure for this
+ sort of calculation.
+
+ """
+
+ def __init__(self, data, leafsize=10):
+ """Construct a kd-tree.
+
+ Parameters:
+ ===========
+
+ data : array-like, shape (n,k)
+ The data points to be indexed. This array is not copied, and
+ so modifying this data will result in bogus results.
+ leafsize : positive integer
+ The number of points at which the algorithm switches over to
+ brute-force.
+ """
+ self.data = np.asarray(data)
+ self.n, self.m = np.shape(self.data)
+ self.leafsize = int(leafsize)
+ if self.leafsize<1:
+ raise ValueError("leafsize must be at least 1")
+ self.maxes = np.amax(self.data,axis=0)
+ self.mins = np.amin(self.data,axis=0)
+
+ self.tree = self.__build(np.arange(self.n), self.maxes, self.mins)
+
+ class node(object):
+ if sys.version_info[0] >= 3:
+ def __lt__(self, other): id(self) < id(other)
+ def __gt__(self, other): id(self) > id(other)
+ def __le__(self, other): id(self) <= id(other)
+ def __ge__(self, other): id(self) >= id(other)
+ def __eq__(self, other): id(self) == id(other)
+ class leafnode(node):
+ def __init__(self, idx):
+ self.idx = idx
+ self.children = len(idx)
+ class innernode(node):
+ def __init__(self, split_dim, split, less, greater):
+ self.split_dim = split_dim
+ self.split = split
+ self.less = less
+ self.greater = greater
+ self.children = less.children+greater.children
+
+ def __build(self, idx, maxes, mins):
+ if len(idx)<=self.leafsize:
+ return KDTree.leafnode(idx)
+ else:
+ data = self.data[idx]
+ #maxes = np.amax(data,axis=0)
+ #mins = np.amin(data,axis=0)
+ d = np.argmax(maxes-mins)
+ maxval = maxes[d]
+ minval = mins[d]
+ if maxval==minval:
+ # all points are identical; warn user?
+ return KDTree.leafnode(idx)
+ data = data[:,d]
+
+ # sliding midpoint rule; see Maneewongvatana and Mount 1999
+ # for arguments that this is a good idea.
+ split = (maxval+minval)/2
+ less_idx = np.nonzero(data<=split)[0]
+ greater_idx = np.nonzero(data>split)[0]
+ if len(less_idx)==0:
+ split = np.amin(data)
+ less_idx = np.nonzero(data<=split)[0]
+ greater_idx = np.nonzero(data>split)[0]
+ if len(greater_idx)==0:
+ split = np.amax(data)
+ less_idx = np.nonzero(data<split)[0]
+ greater_idx = np.nonzero(data>=split)[0]
+ if len(less_idx)==0:
+ # _still_ zero? all must have the same value
+ assert np.all(data==data[0]), "Troublesome data array: %s" % data
+ split = data[0]
+ less_idx = np.arange(len(data)-1)
+ greater_idx = np.array([len(data)-1])
+
+ lessmaxes = np.copy(maxes)
+ lessmaxes[d] = split
+ greatermins = np.copy(mins)
+ greatermins[d] = split
+ return KDTree.innernode(d, split,
+ self.__build(idx[less_idx],lessmaxes,mins),
+ self.__build(idx[greater_idx],maxes,greatermins))
+
+ def __query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
+
+ side_distances = np.maximum(0,np.maximum(x-self.maxes,self.mins-x))
+ if p!=np.inf:
+ side_distances**=p
+ min_distance = np.sum(side_distances)
+ else:
+ min_distance = np.amax(side_distances)
+
+ # priority queue for chasing nodes
+ # entries are:
+ # minimum distance between the cell and the target
+ # distances between the nearest side of the cell and the target
+ # the head node of the cell
+ q = [(min_distance,
+ tuple(side_distances),
+ self.tree)]
+ # priority queue for the nearest neighbors
+ # furthest known neighbor first
+ # entries are (-distance**p, i)
+ neighbors = []
+
+ if eps==0:
+ epsfac=1
+ elif p==np.inf:
+ epsfac = 1/(1+eps)
+ else:
+ epsfac = 1/(1+eps)**p
+
+ if p!=np.inf and distance_upper_bound!=np.inf:
+ distance_upper_bound = distance_upper_bound**p
+
+ while q:
+ min_distance, side_distances, node = heappop(q)
+ if isinstance(node, KDTree.leafnode):
+ # brute-force
+ data = self.data[node.idx]
+ ds = minkowski_distance_p(data,x[np.newaxis,:],p)
+ for i in range(len(ds)):
+ if ds[i]<distance_upper_bound:
+ if len(neighbors)==k:
+ heappop(neighbors)
+ heappush(neighbors, (-ds[i], node.idx[i]))
+ if len(neighbors)==k:
+ distance_upper_bound = -neighbors[0][0]
+ else:
+ # we don't push cells that are too far onto the queue at all,
+ # but since the distance_upper_bound decreases, we might get
+ # here even if the cell's too far
+ if min_distance>distance_upper_bound*epsfac:
+ # since this is the nearest cell, we're done, bail out
+ break
+ # compute minimum distances to the children and push them on
+ if x[node.split_dim]<node.split:
+ near, far = node.less, node.greater
+ else:
+ near, far = node.greater, node.less
+
+ # near child is at the same distance as the current node
+ heappush(q,(min_distance, side_distances, near))
+
+ # far child is further by an amount depending only
+ # on the split value
+ sd = list(side_distances)
+ if p == np.inf:
+ min_distance = max(min_distance, abs(node.split-x[node.split_dim]))
+ elif p == 1:
+ sd[node.split_dim] = np.abs(node.split-x[node.split_dim])
+ min_distance = min_distance - side_distances[node.split_dim] + sd[node.split_dim]
+ else:
+ sd[node.split_dim] = np.abs(node.split-x[node.split_dim])**p
+ min_distance = min_distance - side_distances[node.split_dim] + sd[node.split_dim]
+
+ # far child might be too far, if so, don't bother pushing it
+ if min_distance<=distance_upper_bound*epsfac:
+ heappush(q,(min_distance, tuple(sd), far))
+
+ if p==np.inf:
+ return sorted([(-d,i) for (d,i) in neighbors])
+ else:
+ return sorted([((-d)**(1./p),i) for (d,i) in neighbors])
+
+ def query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
+ """
+ query the kd-tree for nearest neighbors
+
+ Parameters
+ ----------
+
+ x : array-like, last dimension self.m
+ An array of points to query.
+ k : integer
+ The number of nearest neighbors to return.
+ eps : nonnegative float
+ Return approximate nearest neighbors; the kth returned value
+ is guaranteed to be no further than (1+eps) times the
+ distance to the real kth nearest neighbor.
+ p : float, 1<=p<=infinity
+ Which Minkowski p-norm to use.
+ 1 is the sum-of-absolute-values "Manhattan" distance
+ 2 is the usual Euclidean distance
+ infinity is the maximum-coordinate-difference distance
+ distance_upper_bound : nonnegative float
+ Return only neighbors within this distance. This is used to prune
+ tree searches, so if you are doing a series of nearest-neighbor
+ queries, it may help to supply the distance to the nearest neighbor
+ of the most recent point.
+
+ Returns
+ -------
+
+ d : array of floats
+ The distances to the nearest neighbors.
+ If x has shape tuple+(self.m,), then d has shape tuple if
+ k is one, or tuple+(k,) if k is larger than one. Missing
+ neighbors are indicated with infinite distances. If k is None,
+ then d is an object array of shape tuple, containing lists
+ of distances. In either case the hits are sorted by distance
+ (nearest first).
+ i : array of integers
+ The locations of the neighbors in self.data. i is the same
+ shape as d.
+
+ Examples
+ --------
+
+ >>> from scipy.spatial import KDTree
+ >>> x, y = np.mgrid[0:5, 2:8]
+ >>> tree = KDTree(zip(x.ravel(), y.ravel()))
+ >>> tree.data
+ array([[0, 2],
+ [0, 3],
+ [0, 4],
+ [0, 5],
+ [0, 6],
+ [0, 7],
+ [1, 2],
+ [1, 3],
+ [1, 4],
+ [1, 5],
+ [1, 6],
+ [1, 7],
+ [2, 2],
+ [2, 3],
+ [2, 4],
+ [2, 5],
+ [2, 6],
+ [2, 7],
+ [3, 2],
+ [3, 3],
+ [3, 4],
+ [3, 5],
+ [3, 6],
+ [3, 7],
+ [4, 2],
+ [4, 3],
+ [4, 4],
+ [4, 5],
+ [4, 6],
+ [4, 7]])
+ >>> pts = np.array([[0, 0], [2.1, 2.9]])
+ >>> tree.query(pts)
+ (array([ 2. , 0.14142136]), array([ 0, 13]))
+
+ """
+ x = np.asarray(x)
+ if np.shape(x)[-1] != self.m:
+ raise ValueError("x must consist of vectors of length %d but has shape %s" % (self.m, np.shape(x)))
+ if p<1:
+ raise ValueError("Only p-norms with 1<=p<=infinity permitted")
+ retshape = np.shape(x)[:-1]
+ if retshape!=():
+ if k is None:
+ dd = np.empty(retshape,dtype=np.object)
+ ii = np.empty(retshape,dtype=np.object)
+ elif k>1:
+ dd = np.empty(retshape+(k,),dtype=np.float)
+ dd.fill(np.inf)
+ ii = np.empty(retshape+(k,),dtype=np.int)
+ ii.fill(self.n)
+ elif k==1:
+ dd = np.empty(retshape,dtype=np.float)
+ dd.fill(np.inf)
+ ii = np.empty(retshape,dtype=np.int)
+ ii.fill(self.n)
+ else:
+ raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
+ for c in np.ndindex(retshape):
+ hits = self.__query(x[c], k=k, p=p, distance_upper_bound=distance_upper_bound)
+ if k is None:
+ dd[c] = [d for (d,i) in hits]
+ ii[c] = [i for (d,i) in hits]
+ elif k>1:
+ for j in range(len(hits)):
+ dd[c+(j,)], ii[c+(j,)] = hits[j]
+ elif k==1:
+ if len(hits)>0:
+ dd[c], ii[c] = hits[0]
+ else:
+ dd[c] = np.inf
+ ii[c] = self.n
+ return dd, ii
+ else:
+ hits = self.__query(x, k=k, p=p, distance_upper_bound=distance_upper_bound)
+ if k is None:
+ return [d for (d,i) in hits], [i for (d,i) in hits]
+ elif k==1:
+ if len(hits)>0:
+ return hits[0]
+ else:
+ return np.inf, self.n
+ elif k>1:
+ dd = np.empty(k,dtype=np.float)
+ dd.fill(np.inf)
+ ii = np.empty(k,dtype=np.int)
+ ii.fill(self.n)
+ for j in range(len(hits)):
+ dd[j], ii[j] = hits[j]
+ return dd, ii
+ else:
+ raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
+
+
+ def __query_ball_point(self, x, r, p=2., eps=0):
+ R = Rectangle(self.maxes, self.mins)
+
+ def traverse_checking(node, rect):
+ if rect.min_distance_point(x,p)>=r/(1.+eps):
+ return []
+ elif rect.max_distance_point(x,p)<r*(1.+eps):
+ return traverse_no_checking(node)
+ elif isinstance(node, KDTree.leafnode):
+ d = self.data[node.idx]
+ return node.idx[minkowski_distance(d,x,p)<=r].tolist()
+ else:
+ less, greater = rect.split(node.split_dim, node.split)
+ return traverse_checking(node.less, less)+traverse_checking(node.greater, greater)
+ def traverse_no_checking(node):
+ if isinstance(node, KDTree.leafnode):
+
+ return node.idx.tolist()
+ else:
+ return traverse_no_checking(node.less)+traverse_no_checking(node.greater)
+
+ return traverse_checking(self.tree, R)
+
+ def query_ball_point(self, x, r, p=2., eps=0):
+ """Find all points within r of x
+
+ Parameters
+ ==========
+
+ x : array_like, shape tuple + (self.m,)
+ The point or points to search for neighbors of
+ r : positive float
+ The radius of points to return
+ p : float 1<=p<=infinity
+ Which Minkowski p-norm to use
+ eps : nonnegative float
+ Approximate search. Branches of the tree are not explored
+ if their nearest points are further than r/(1+eps), and branches
+ are added in bulk if their furthest points are nearer than r*(1+eps).
+
+ Returns
+ =======
+
+ results : list or array of lists
+ If x is a single point, returns a list of the indices of the neighbors
+ of x. If x is an array of points, returns an object array of shape tuple
+ containing lists of neighbors.
+
+
+ Note: if you have many points whose neighbors you want to find, you may save
+ substantial amounts of time by putting them in a KDTree and using query_ball_tree.
+ """
+ x = np.asarray(x)
+ if x.shape[-1]!=self.m:
+ raise ValueError("Searching for a %d-dimensional point in a %d-dimensional KDTree" % (x.shape[-1],self.m))
+ if len(x.shape)==1:
+ return self.__query_ball_point(x,r,p,eps)
+ else:
+ retshape = x.shape[:-1]
+ result = np.empty(retshape,dtype=np.object)
+ for c in np.ndindex(retshape):
+ result[c] = self.__query_ball_point(x[c], r, p=p, eps=eps)
+ return result
+
+ def query_ball_tree(self, other, r, p=2., eps=0):
+ """Find all pairs of points whose distance is at most r
+
+ Parameters
+ ==========
+
+ other : KDTree
+ The tree containing points to search against
+ r : positive float
+ The maximum distance
+ p : float 1<=p<=infinity
+ Which Minkowski norm to use
+ eps : nonnegative float
+ Approximate search. Branches of the tree are not explored
+ if their nearest points are further than r/(1+eps), and branches
+ are added in bulk if their furthest points are nearer than r*(1+eps).
+
+ Returns
+ =======
+
+ results : list of lists
+ For each element self.data[i] of this tree, results[i] is a list of the
+ indices of its neighbors in other.data.
+ """
+ results = [[] for i in range(self.n)]
+ def traverse_checking(node1, rect1, node2, rect2):
+ if rect1.min_distance_rectangle(rect2, p)>r/(1.+eps):
+ return
+ elif rect1.max_distance_rectangle(rect2, p)<r*(1.+eps):
+ traverse_no_checking(node1, node2)
+ elif isinstance(node1, KDTree.leafnode):
+ if isinstance(node2, KDTree.leafnode):
+ d = other.data[node2.idx]
+ for i in node1.idx:
+ results[i] += node2.idx[minkowski_distance(d,self.data[i],p)<=r].tolist()
+ else:
+ less, greater = rect2.split(node2.split_dim, node2.split)
+ traverse_checking(node1,rect1,node2.less,less)
+ traverse_checking(node1,rect1,node2.greater,greater)
+ elif isinstance(node2, KDTree.leafnode):
+ less, greater = rect1.split(node1.split_dim, node1.split)
+ traverse_checking(node1.less,less,node2,rect2)
+ traverse_checking(node1.greater,greater,node2,rect2)
+ else:
+ less1, greater1 = rect1.split(node1.split_dim, node1.split)
+ less2, greater2 = rect2.split(node2.split_dim, node2.split)
+ traverse_checking(node1.less,less1,node2.less,less2)
+ traverse_checking(node1.less,less1,node2.greater,greater2)
+ traverse_checking(node1.greater,greater1,node2.less,less2)
+ traverse_checking(node1.greater,greater1,node2.greater,greater2)
+
+ def traverse_no_checking(node1, node2):
+ if isinstance(node1, KDTree.leafnode):
+ if isinstance(node2, KDTree.leafnode):
+ for i in node1.idx:
+ results[i] += node2.idx.tolist()
+ else:
+ traverse_no_checking(node1, node2.less)
+ traverse_no_checking(node1, node2.greater)
+ else:
+ traverse_no_checking(node1.less, node2)
+ traverse_no_checking(node1.greater, node2)
+
+ traverse_checking(self.tree, Rectangle(self.maxes, self.mins),
+ other.tree, Rectangle(other.maxes, other.mins))
+ return results
+
+ def query_pairs(self, r, p=2., eps=0):
+ """Find all pairs of points whose distance is at most r
+
+ Parameters
+ ==========
+
+ r : positive float
+ The maximum distance
+ p : float 1<=p<=infinity
+ Which Minkowski norm to use
+ eps : nonnegative float
+ Approximate search. Branches of the tree are not explored
+ if their nearest points are further than r/(1+eps), and branches
+ are added in bulk if their furthest points are nearer than r*(1+eps).
+
+ Returns
+ =======
+
+ results : set
+ set of pairs (i,j), i<j, for which the corresponing positions are
+ close.
+
+ """
+ results = set()
+ visited = set()
+ def test_set_visited(node1, node2):
+ i, j = sorted((id(node1),id(node2)))
+ if (i,j) in visited:
+ return True
+ else:
+ visited.add((i,j))
+ return False
+ def traverse_checking(node1, rect1, node2, rect2):
+ if test_set_visited(node1, node2):
+ return
+
+ if id(node2)<id(node1):
+ # This node pair will be visited in the other order
+ #return
+ pass
+
+ if isinstance(node1, KDTree.leafnode):
+ if isinstance(node2, KDTree.leafnode):
+ d = self.data[node2.idx]
+ for i in node1.idx:
+ for j in node2.idx[minkowski_distance(d,self.data[i],p)<=r]:
+ if i<j:
+ results.add((i,j))
+ elif j<i:
+ results.add((j,i))
+ else:
+ less, greater = rect2.split(node2.split_dim, node2.split)
+ traverse_checking(node1,rect1,node2.less,less)
+ traverse_checking(node1,rect1,node2.greater,greater)
+ elif isinstance(node2, KDTree.leafnode):
+ less, greater = rect1.split(node1.split_dim, node1.split)
+ traverse_checking(node1.less,less,node2,rect2)
+ traverse_checking(node1.greater,greater,node2,rect2)
+ elif rect1.min_distance_rectangle(rect2, p)>r/(1.+eps):
+ return
+ elif rect1.max_distance_rectangle(rect2, p)<r*(1.+eps):
+ traverse_no_checking(node1.less, node2)
+ traverse_no_checking(node1.greater, node2)
+ else:
+ less1, greater1 = rect1.split(node1.split_dim, node1.split)
+ less2, greater2 = rect2.split(node2.split_dim, node2.split)
+ traverse_checking(node1.less,less1,node2.less,less2)
+ traverse_checking(node1.less,less1,node2.greater,greater2)
+ traverse_checking(node1.greater,greater1,node2.less,less2)
+ traverse_checking(node1.greater,greater1,node2.greater,greater2)
+
+ def traverse_no_checking(node1, node2):
+ if test_set_visited(node1, node2):
+ return
+
+ if id(node2)<id(node1):
+ # This node pair will be visited in the other order
+ #return
+ pass
+ if isinstance(node1, KDTree.leafnode):
+ if isinstance(node2, KDTree.leafnode):
+ for i in node1.idx:
+ for j in node2.idx:
+ if i<j:
+ results.add((i,j))
+ elif j<i:
+ results.add((j,i))
+ else:
+ traverse_no_checking(node1, node2.less)
+ traverse_no_checking(node1, node2.greater)
+ else:
+ traverse_no_checking(node1.less, node2)
+ traverse_no_checking(node1.greater, node2)
+
+ traverse_checking(self.tree, Rectangle(self.maxes, self.mins),
+ self.tree, Rectangle(self.maxes, self.mins))
+ return results
+
+
+ def count_neighbors(self, other, r, p=2.):
+ """Count how many nearby pairs can be formed.
+
+ Count the number of pairs (x1,x2) can be formed, with x1 drawn
+ from self and x2 drawn from other, and where distance(x1,x2,p)<=r.
+ This is the "two-point correlation" described in Gray and Moore 2000,
+ "N-body problems in statistical learning", and the code here is based
+ on their algorithm.
+
+ Parameters
+ ==========
+
+ other : KDTree
+
+ r : float or one-dimensional array of floats
+ The radius to produce a count for. Multiple radii are searched with a single
+ tree traversal.
+ p : float, 1<=p<=infinity
+ Which Minkowski p-norm to use
+
+ Returns
+ =======
+
+ result : integer or one-dimensional array of integers
+ The number of pairs. Note that this is internally stored in a numpy int,
+ and so may overflow if very large (two billion).
+ """
+
+ def traverse(node1, rect1, node2, rect2, idx):
+ min_r = rect1.min_distance_rectangle(rect2,p)
+ max_r = rect1.max_distance_rectangle(rect2,p)
+ c_greater = r[idx]>max_r
+ result[idx[c_greater]] += node1.children*node2.children
+ idx = idx[(min_r<=r[idx]) & (r[idx]<=max_r)]
+ if len(idx)==0:
+ return
+
+ if isinstance(node1,KDTree.leafnode):
+ if isinstance(node2,KDTree.leafnode):
+ ds = minkowski_distance(self.data[node1.idx][:,np.newaxis,:],
+ other.data[node2.idx][np.newaxis,:,:],
+ p).ravel()
+ ds.sort()
+ result[idx] += np.searchsorted(ds,r[idx],side='right')
+ else:
+ less, greater = rect2.split(node2.split_dim, node2.split)
+ traverse(node1, rect1, node2.less, less, idx)
+ traverse(node1, rect1, node2.greater, greater, idx)
+ else:
+ if isinstance(node2,KDTree.leafnode):
+ less, greater = rect1.split(node1.split_dim, node1.split)
+ traverse(node1.less, less, node2, rect2, idx)
+ traverse(node1.greater, greater, node2, rect2, idx)
+ else:
+ less1, greater1 = rect1.split(node1.split_dim, node1.split)
+ less2, greater2 = rect2.split(node2.split_dim, node2.split)
+ traverse(node1.less,less1,node2.less,less2,idx)
+ traverse(node1.less,less1,node2.greater,greater2,idx)
+ traverse(node1.greater,greater1,node2.less,less2,idx)
+ traverse(node1.greater,greater1,node2.greater,greater2,idx)
+ R1 = Rectangle(self.maxes, self.mins)
+ R2 = Rectangle(other.maxes, other.mins)
+ if np.shape(r) == ():
+ r = np.array([r])
+ result = np.zeros(1,dtype=int)
+ traverse(self.tree, R1, other.tree, R2, np.arange(1))
+ return result[0]
+ elif len(np.shape(r))==1:
+ r = np.asarray(r)
+ n, = r.shape
+ result = np.zeros(n,dtype=int)
+ traverse(self.tree, R1, other.tree, R2, np.arange(n))
+ return result
+ else:
+ raise ValueError("r must be either a single value or a one-dimensional array of values")
+
+ def sparse_distance_matrix(self, other, max_distance, p=2.):
+ """Compute a sparse distance matrix
+
+ Computes a distance matrix between two KDTrees, leaving as zero
+ any distance greater than max_distance.
+
+ Parameters
+ ==========
+
+ other : KDTree
+
+ max_distance : positive float
+
+ Returns
+ =======
+
+ result : dok_matrix
+ Sparse matrix representing the results in "dictionary of keys" format.
+ """
+ result = scipy.sparse.dok_matrix((self.n,other.n))
+
+ def traverse(node1, rect1, node2, rect2):
+ if rect1.min_distance_rectangle(rect2, p)>max_distance:
+ return
+ elif isinstance(node1, KDTree.leafnode):
+ if isinstance(node2, KDTree.leafnode):
+ for i in node1.idx:
+ for j in node2.idx:
+ d = minkowski_distance(self.data[i],other.data[j],p)
+ if d<=max_distance:
+ result[i,j] = d
+ else:
+ less, greater = rect2.split(node2.split_dim, node2.split)
+ traverse(node1,rect1,node2.less,less)
+ traverse(node1,rect1,node2.greater,greater)
+ elif isinstance(node2, KDTree.leafnode):
+ less, greater = rect1.split(node1.split_dim, node1.split)
+ traverse(node1.less,less,node2,rect2)
+ traverse(node1.greater,greater,node2,rect2)
+ else:
+ less1, greater1 = rect1.split(node1.split_dim, node1.split)
+ less2, greater2 = rect2.split(node2.split_dim, node2.split)
+ traverse(node1.less,less1,node2.less,less2)
+ traverse(node1.less,less1,node2.greater,greater2)
+ traverse(node1.greater,greater1,node2.less,less2)
+ traverse(node1.greater,greater1,node2.greater,greater2)
+ traverse(self.tree, Rectangle(self.maxes, self.mins),
+ other.tree, Rectangle(other.maxes, other.mins))
+
+ return result
+
+
+def distance_matrix(x,y,p=2,threshold=1000000):
+ """Compute the distance matrix.
+
+ Computes the matrix of all pairwise distances.
+
+ Parameters
+ ==========
+
+ x : array-like, m by k
+ y : array-like, n by k
+ p : float 1<=p<=infinity
+ Which Minkowski p-norm to use.
+ threshold : positive integer
+ If m*n*k>threshold use a python loop instead of creating
+ a very large temporary.
+
+ Returns
+ =======
+
+ result : array-like, m by n
+
+
+ """
+
+ x = np.asarray(x)
+ m, k = x.shape
+ y = np.asarray(y)
+ n, kk = y.shape
+
+ if k != kk:
+ raise ValueError("x contains %d-dimensional vectors but y contains %d-dimensional vectors" % (k, kk))
+
+ if m*n*k <= threshold:
+ return minkowski_distance(x[:,np.newaxis,:],y[np.newaxis,:,:],p)
+ else:
+ result = np.empty((m,n),dtype=np.float) #FIXME: figure out the best dtype
+ if m<n:
+ for i in range(m):
+ result[i,:] = minkowski_distance(x[i],y,p)
+ else:
+ for j in range(n):
+ result[:,j] = minkowski_distance(x,y[j],p)
+ return result
diff -r 8960f029b62f4c703889a784449014ddaf7e0a40 -r 6938c5757a9246082a90a7d6e0678d052a1e7489 yt/utilities/pykdtree.py
--- a/yt/utilities/pykdtree.py
+++ /dev/null
@@ -1,823 +0,0 @@
-# Copyright Anne M. Archibald 2008
-# Released under the scipy license
-import sys
-import numpy as np
-from heapq import heappush, heappop
-
-def minkowski_distance_p(x,y,p=2):
- """Compute the pth power of the L**p distance between x and y
-
- For efficiency, this function computes the L**p distance but does
- not extract the pth root. If p is 1 or infinity, this is equal to
- the actual L**p distance.
- """
- x = np.asarray(x)
- y = np.asarray(y)
- if p==np.inf:
- return np.amax(np.abs(y-x),axis=-1)
- elif p==1:
- return np.sum(np.abs(y-x),axis=-1)
- else:
- return np.sum(np.abs(y-x)**p,axis=-1)
-def minkowski_distance(x,y,p=2):
- """Compute the L**p distance between x and y"""
- x = np.asarray(x)
- y = np.asarray(y)
- if p==np.inf or p==1:
- return minkowski_distance_p(x,y,p)
- else:
- return minkowski_distance_p(x,y,p)**(1./p)
-
-class Rectangle(object):
- """Hyperrectangle class.
-
- Represents a Cartesian product of intervals.
- """
- def __init__(self, maxes, mins):
- """Construct a hyperrectangle."""
- self.maxes = np.maximum(maxes,mins).astype(np.float)
- self.mins = np.minimum(maxes,mins).astype(np.float)
- self.m, = self.maxes.shape
-
- def __repr__(self):
- return "<Rectangle %s>" % zip(self.mins, self.maxes)
-
- def volume(self):
- """Total volume."""
- return np.prod(self.maxes-self.mins)
-
- def split(self, d, split):
- """Produce two hyperrectangles by splitting along axis d.
-
- In general, if you need to compute maximum and minimum
- distances to the children, it can be done more efficiently
- by updating the maximum and minimum distances to the parent.
- """ # FIXME: do this
- mid = np.copy(self.maxes)
- mid[d] = split
- less = Rectangle(self.mins, mid)
- mid = np.copy(self.mins)
- mid[d] = split
- greater = Rectangle(mid, self.maxes)
- return less, greater
-
- def min_distance_point(self, x, p=2.):
- """Compute the minimum distance between x and a point in the hyperrectangle."""
- return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-x,x-self.maxes)),p)
-
- def max_distance_point(self, x, p=2.):
- """Compute the maximum distance between x and a point in the hyperrectangle."""
- return minkowski_distance(0, np.maximum(self.maxes-x,x-self.mins),p)
-
- def min_distance_rectangle(self, other, p=2.):
- """Compute the minimum distance between points in the two hyperrectangles."""
- return minkowski_distance(0, np.maximum(0,np.maximum(self.mins-other.maxes,other.mins-self.maxes)),p)
-
- def max_distance_rectangle(self, other, p=2.):
- """Compute the maximum distance between points in the two hyperrectangles."""
- return minkowski_distance(0, np.maximum(self.maxes-other.mins,other.maxes-self.mins),p)
-
-
-class KDTree(object):
- """
- kd-tree for quick nearest-neighbor lookup
-
- This class provides an index into a set of k-dimensional points
- which can be used to rapidly look up the nearest neighbors of any
- point.
-
- The algorithm used is described in Maneewongvatana and Mount 1999.
- The general idea is that the kd-tree is a binary tree, each of whose
- nodes represents an axis-aligned hyperrectangle. Each node specifies
- an axis and splits the set of points based on whether their coordinate
- along that axis is greater than or less than a particular value.
-
- During construction, the axis and splitting point are chosen by the
- "sliding midpoint" rule, which ensures that the cells do not all
- become long and thin.
-
- The tree can be queried for the r closest neighbors of any given point
- (optionally returning only those within some maximum distance of the
- point). It can also be queried, with a substantial gain in efficiency,
- for the r approximate closest neighbors.
-
- For large dimensions (20 is already large) do not expect this to run
- significantly faster than brute force. High-dimensional nearest-neighbor
- queries are a substantial open problem in computer science.
-
- The tree also supports all-neighbors queries, both with arrays of points
- and with other kd-trees. These do use a reasonably efficient algorithm,
- but the kd-tree is not necessarily the best data structure for this
- sort of calculation.
-
- """
-
- def __init__(self, data, leafsize=10):
- """Construct a kd-tree.
-
- Parameters:
- ===========
-
- data : array-like, shape (n,k)
- The data points to be indexed. This array is not copied, and
- so modifying this data will result in bogus results.
- leafsize : positive integer
- The number of points at which the algorithm switches over to
- brute-force.
- """
- self.data = np.asarray(data)
- self.n, self.m = np.shape(self.data)
- self.leafsize = int(leafsize)
- if self.leafsize<1:
- raise ValueError("leafsize must be at least 1")
- self.maxes = np.amax(self.data,axis=0)
- self.mins = np.amin(self.data,axis=0)
-
- self.tree = self.__build(np.arange(self.n), self.maxes, self.mins)
-
- class node(object):
- if sys.version_info[0] >= 3:
- def __lt__(self, other): id(self) < id(other)
- def __gt__(self, other): id(self) > id(other)
- def __le__(self, other): id(self) <= id(other)
- def __ge__(self, other): id(self) >= id(other)
- def __eq__(self, other): id(self) == id(other)
- class leafnode(node):
- def __init__(self, idx):
- self.idx = idx
- self.children = len(idx)
- class innernode(node):
- def __init__(self, split_dim, split, less, greater):
- self.split_dim = split_dim
- self.split = split
- self.less = less
- self.greater = greater
- self.children = less.children+greater.children
-
- def __build(self, idx, maxes, mins):
- if len(idx)<=self.leafsize:
- return KDTree.leafnode(idx)
- else:
- data = self.data[idx]
- #maxes = np.amax(data,axis=0)
- #mins = np.amin(data,axis=0)
- d = np.argmax(maxes-mins)
- maxval = maxes[d]
- minval = mins[d]
- if maxval==minval:
- # all points are identical; warn user?
- return KDTree.leafnode(idx)
- data = data[:,d]
-
- # sliding midpoint rule; see Maneewongvatana and Mount 1999
- # for arguments that this is a good idea.
- split = (maxval+minval)/2
- less_idx = np.nonzero(data<=split)[0]
- greater_idx = np.nonzero(data>split)[0]
- if len(less_idx)==0:
- split = np.amin(data)
- less_idx = np.nonzero(data<=split)[0]
- greater_idx = np.nonzero(data>split)[0]
- if len(greater_idx)==0:
- split = np.amax(data)
- less_idx = np.nonzero(data<split)[0]
- greater_idx = np.nonzero(data>=split)[0]
- if len(less_idx)==0:
- # _still_ zero? all must have the same value
- assert np.all(data==data[0]), "Troublesome data array: %s" % data
- split = data[0]
- less_idx = np.arange(len(data)-1)
- greater_idx = np.array([len(data)-1])
-
- lessmaxes = np.copy(maxes)
- lessmaxes[d] = split
- greatermins = np.copy(mins)
- greatermins[d] = split
- return KDTree.innernode(d, split,
- self.__build(idx[less_idx],lessmaxes,mins),
- self.__build(idx[greater_idx],maxes,greatermins))
-
- def __query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
-
- side_distances = np.maximum(0,np.maximum(x-self.maxes,self.mins-x))
- if p!=np.inf:
- side_distances**=p
- min_distance = np.sum(side_distances)
- else:
- min_distance = np.amax(side_distances)
-
- # priority queue for chasing nodes
- # entries are:
- # minimum distance between the cell and the target
- # distances between the nearest side of the cell and the target
- # the head node of the cell
- q = [(min_distance,
- tuple(side_distances),
- self.tree)]
- # priority queue for the nearest neighbors
- # furthest known neighbor first
- # entries are (-distance**p, i)
- neighbors = []
-
- if eps==0:
- epsfac=1
- elif p==np.inf:
- epsfac = 1/(1+eps)
- else:
- epsfac = 1/(1+eps)**p
-
- if p!=np.inf and distance_upper_bound!=np.inf:
- distance_upper_bound = distance_upper_bound**p
-
- while q:
- min_distance, side_distances, node = heappop(q)
- if isinstance(node, KDTree.leafnode):
- # brute-force
- data = self.data[node.idx]
- ds = minkowski_distance_p(data,x[np.newaxis,:],p)
- for i in range(len(ds)):
- if ds[i]<distance_upper_bound:
- if len(neighbors)==k:
- heappop(neighbors)
- heappush(neighbors, (-ds[i], node.idx[i]))
- if len(neighbors)==k:
- distance_upper_bound = -neighbors[0][0]
- else:
- # we don't push cells that are too far onto the queue at all,
- # but since the distance_upper_bound decreases, we might get
- # here even if the cell's too far
- if min_distance>distance_upper_bound*epsfac:
- # since this is the nearest cell, we're done, bail out
- break
- # compute minimum distances to the children and push them on
- if x[node.split_dim]<node.split:
- near, far = node.less, node.greater
- else:
- near, far = node.greater, node.less
-
- # near child is at the same distance as the current node
- heappush(q,(min_distance, side_distances, near))
-
- # far child is further by an amount depending only
- # on the split value
- sd = list(side_distances)
- if p == np.inf:
- min_distance = max(min_distance, abs(node.split-x[node.split_dim]))
- elif p == 1:
- sd[node.split_dim] = np.abs(node.split-x[node.split_dim])
- min_distance = min_distance - side_distances[node.split_dim] + sd[node.split_dim]
- else:
- sd[node.split_dim] = np.abs(node.split-x[node.split_dim])**p
- min_distance = min_distance - side_distances[node.split_dim] + sd[node.split_dim]
-
- # far child might be too far, if so, don't bother pushing it
- if min_distance<=distance_upper_bound*epsfac:
- heappush(q,(min_distance, tuple(sd), far))
-
- if p==np.inf:
- return sorted([(-d,i) for (d,i) in neighbors])
- else:
- return sorted([((-d)**(1./p),i) for (d,i) in neighbors])
-
- def query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):
- """
- query the kd-tree for nearest neighbors
-
- Parameters
- ----------
-
- x : array-like, last dimension self.m
- An array of points to query.
- k : integer
- The number of nearest neighbors to return.
- eps : nonnegative float
- Return approximate nearest neighbors; the kth returned value
- is guaranteed to be no further than (1+eps) times the
- distance to the real kth nearest neighbor.
- p : float, 1<=p<=infinity
- Which Minkowski p-norm to use.
- 1 is the sum-of-absolute-values "Manhattan" distance
- 2 is the usual Euclidean distance
- infinity is the maximum-coordinate-difference distance
- distance_upper_bound : nonnegative float
- Return only neighbors within this distance. This is used to prune
- tree searches, so if you are doing a series of nearest-neighbor
- queries, it may help to supply the distance to the nearest neighbor
- of the most recent point.
-
- Returns
- -------
-
- d : array of floats
- The distances to the nearest neighbors.
- If x has shape tuple+(self.m,), then d has shape tuple if
- k is one, or tuple+(k,) if k is larger than one. Missing
- neighbors are indicated with infinite distances. If k is None,
- then d is an object array of shape tuple, containing lists
- of distances. In either case the hits are sorted by distance
- (nearest first).
- i : array of integers
- The locations of the neighbors in self.data. i is the same
- shape as d.
-
- Examples
- --------
-
- >>> from scipy.spatial import KDTree
- >>> x, y = np.mgrid[0:5, 2:8]
- >>> tree = KDTree(zip(x.ravel(), y.ravel()))
- >>> tree.data
- array([[0, 2],
- [0, 3],
- [0, 4],
- [0, 5],
- [0, 6],
- [0, 7],
- [1, 2],
- [1, 3],
- [1, 4],
- [1, 5],
- [1, 6],
- [1, 7],
- [2, 2],
- [2, 3],
- [2, 4],
- [2, 5],
- [2, 6],
- [2, 7],
- [3, 2],
- [3, 3],
- [3, 4],
- [3, 5],
- [3, 6],
- [3, 7],
- [4, 2],
- [4, 3],
- [4, 4],
- [4, 5],
- [4, 6],
- [4, 7]])
- >>> pts = np.array([[0, 0], [2.1, 2.9]])
- >>> tree.query(pts)
- (array([ 2. , 0.14142136]), array([ 0, 13]))
-
- """
- x = np.asarray(x)
- if np.shape(x)[-1] != self.m:
- raise ValueError("x must consist of vectors of length %d but has shape %s" % (self.m, np.shape(x)))
- if p<1:
- raise ValueError("Only p-norms with 1<=p<=infinity permitted")
- retshape = np.shape(x)[:-1]
- if retshape!=():
- if k is None:
- dd = np.empty(retshape,dtype=np.object)
- ii = np.empty(retshape,dtype=np.object)
- elif k>1:
- dd = np.empty(retshape+(k,),dtype=np.float)
- dd.fill(np.inf)
- ii = np.empty(retshape+(k,),dtype=np.int)
- ii.fill(self.n)
- elif k==1:
- dd = np.empty(retshape,dtype=np.float)
- dd.fill(np.inf)
- ii = np.empty(retshape,dtype=np.int)
- ii.fill(self.n)
- else:
- raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
- for c in np.ndindex(retshape):
- hits = self.__query(x[c], k=k, p=p, distance_upper_bound=distance_upper_bound)
- if k is None:
- dd[c] = [d for (d,i) in hits]
- ii[c] = [i for (d,i) in hits]
- elif k>1:
- for j in range(len(hits)):
- dd[c+(j,)], ii[c+(j,)] = hits[j]
- elif k==1:
- if len(hits)>0:
- dd[c], ii[c] = hits[0]
- else:
- dd[c] = np.inf
- ii[c] = self.n
- return dd, ii
- else:
- hits = self.__query(x, k=k, p=p, distance_upper_bound=distance_upper_bound)
- if k is None:
- return [d for (d,i) in hits], [i for (d,i) in hits]
- elif k==1:
- if len(hits)>0:
- return hits[0]
- else:
- return np.inf, self.n
- elif k>1:
- dd = np.empty(k,dtype=np.float)
- dd.fill(np.inf)
- ii = np.empty(k,dtype=np.int)
- ii.fill(self.n)
- for j in range(len(hits)):
- dd[j], ii[j] = hits[j]
- return dd, ii
- else:
- raise ValueError("Requested %s nearest neighbors; acceptable numbers are integers greater than or equal to one, or None")
-
-
- def __query_ball_point(self, x, r, p=2., eps=0):
- R = Rectangle(self.maxes, self.mins)
-
- def traverse_checking(node, rect):
- if rect.min_distance_point(x,p)>=r/(1.+eps):
- return []
- elif rect.max_distance_point(x,p)<r*(1.+eps):
- return traverse_no_checking(node)
- elif isinstance(node, KDTree.leafnode):
- d = self.data[node.idx]
- return node.idx[minkowski_distance(d,x,p)<=r].tolist()
- else:
- less, greater = rect.split(node.split_dim, node.split)
- return traverse_checking(node.less, less)+traverse_checking(node.greater, greater)
- def traverse_no_checking(node):
- if isinstance(node, KDTree.leafnode):
-
- return node.idx.tolist()
- else:
- return traverse_no_checking(node.less)+traverse_no_checking(node.greater)
-
- return traverse_checking(self.tree, R)
-
- def query_ball_point(self, x, r, p=2., eps=0):
- """Find all points within r of x
-
- Parameters
- ==========
-
- x : array_like, shape tuple + (self.m,)
- The point or points to search for neighbors of
- r : positive float
- The radius of points to return
- p : float 1<=p<=infinity
- Which Minkowski p-norm to use
- eps : nonnegative float
- Approximate search. Branches of the tree are not explored
- if their nearest points are further than r/(1+eps), and branches
- are added in bulk if their furthest points are nearer than r*(1+eps).
-
- Returns
- =======
-
- results : list or array of lists
- If x is a single point, returns a list of the indices of the neighbors
- of x. If x is an array of points, returns an object array of shape tuple
- containing lists of neighbors.
-
-
- Note: if you have many points whose neighbors you want to find, you may save
- substantial amounts of time by putting them in a KDTree and using query_ball_tree.
- """
- x = np.asarray(x)
- if x.shape[-1]!=self.m:
- raise ValueError("Searching for a %d-dimensional point in a %d-dimensional KDTree" % (x.shape[-1],self.m))
- if len(x.shape)==1:
- return self.__query_ball_point(x,r,p,eps)
- else:
- retshape = x.shape[:-1]
- result = np.empty(retshape,dtype=np.object)
- for c in np.ndindex(retshape):
- result[c] = self.__query_ball_point(x[c], r, p=p, eps=eps)
- return result
-
- def query_ball_tree(self, other, r, p=2., eps=0):
- """Find all pairs of points whose distance is at most r
-
- Parameters
- ==========
-
- other : KDTree
- The tree containing points to search against
- r : positive float
- The maximum distance
- p : float 1<=p<=infinity
- Which Minkowski norm to use
- eps : nonnegative float
- Approximate search. Branches of the tree are not explored
- if their nearest points are further than r/(1+eps), and branches
- are added in bulk if their furthest points are nearer than r*(1+eps).
-
- Returns
- =======
-
- results : list of lists
- For each element self.data[i] of this tree, results[i] is a list of the
- indices of its neighbors in other.data.
- """
- results = [[] for i in range(self.n)]
- def traverse_checking(node1, rect1, node2, rect2):
- if rect1.min_distance_rectangle(rect2, p)>r/(1.+eps):
- return
- elif rect1.max_distance_rectangle(rect2, p)<r*(1.+eps):
- traverse_no_checking(node1, node2)
- elif isinstance(node1, KDTree.leafnode):
- if isinstance(node2, KDTree.leafnode):
- d = other.data[node2.idx]
- for i in node1.idx:
- results[i] += node2.idx[minkowski_distance(d,self.data[i],p)<=r].tolist()
- else:
- less, greater = rect2.split(node2.split_dim, node2.split)
- traverse_checking(node1,rect1,node2.less,less)
- traverse_checking(node1,rect1,node2.greater,greater)
- elif isinstance(node2, KDTree.leafnode):
- less, greater = rect1.split(node1.split_dim, node1.split)
- traverse_checking(node1.less,less,node2,rect2)
- traverse_checking(node1.greater,greater,node2,rect2)
- else:
- less1, greater1 = rect1.split(node1.split_dim, node1.split)
- less2, greater2 = rect2.split(node2.split_dim, node2.split)
- traverse_checking(node1.less,less1,node2.less,less2)
- traverse_checking(node1.less,less1,node2.greater,greater2)
- traverse_checking(node1.greater,greater1,node2.less,less2)
- traverse_checking(node1.greater,greater1,node2.greater,greater2)
-
- def traverse_no_checking(node1, node2):
- if isinstance(node1, KDTree.leafnode):
- if isinstance(node2, KDTree.leafnode):
- for i in node1.idx:
- results[i] += node2.idx.tolist()
- else:
- traverse_no_checking(node1, node2.less)
- traverse_no_checking(node1, node2.greater)
- else:
- traverse_no_checking(node1.less, node2)
- traverse_no_checking(node1.greater, node2)
-
- traverse_checking(self.tree, Rectangle(self.maxes, self.mins),
- other.tree, Rectangle(other.maxes, other.mins))
- return results
-
- def query_pairs(self, r, p=2., eps=0):
- """Find all pairs of points whose distance is at most r
-
- Parameters
- ==========
-
- r : positive float
- The maximum distance
- p : float 1<=p<=infinity
- Which Minkowski norm to use
- eps : nonnegative float
- Approximate search. Branches of the tree are not explored
- if their nearest points are further than r/(1+eps), and branches
- are added in bulk if their furthest points are nearer than r*(1+eps).
-
- Returns
- =======
-
- results : set
- set of pairs (i,j), i<j, for which the corresponing positions are
- close.
-
- """
- results = set()
- visited = set()
- def test_set_visited(node1, node2):
- i, j = sorted((id(node1),id(node2)))
- if (i,j) in visited:
- return True
- else:
- visited.add((i,j))
- return False
- def traverse_checking(node1, rect1, node2, rect2):
- if test_set_visited(node1, node2):
- return
-
- if id(node2)<id(node1):
- # This node pair will be visited in the other order
- #return
- pass
-
- if isinstance(node1, KDTree.leafnode):
- if isinstance(node2, KDTree.leafnode):
- d = self.data[node2.idx]
- for i in node1.idx:
- for j in node2.idx[minkowski_distance(d,self.data[i],p)<=r]:
- if i<j:
- results.add((i,j))
- elif j<i:
- results.add((j,i))
- else:
- less, greater = rect2.split(node2.split_dim, node2.split)
- traverse_checking(node1,rect1,node2.less,less)
- traverse_checking(node1,rect1,node2.greater,greater)
- elif isinstance(node2, KDTree.leafnode):
- less, greater = rect1.split(node1.split_dim, node1.split)
- traverse_checking(node1.less,less,node2,rect2)
- traverse_checking(node1.greater,greater,node2,rect2)
- elif rect1.min_distance_rectangle(rect2, p)>r/(1.+eps):
- return
- elif rect1.max_distance_rectangle(rect2, p)<r*(1.+eps):
- traverse_no_checking(node1.less, node2)
- traverse_no_checking(node1.greater, node2)
- else:
- less1, greater1 = rect1.split(node1.split_dim, node1.split)
- less2, greater2 = rect2.split(node2.split_dim, node2.split)
- traverse_checking(node1.less,less1,node2.less,less2)
- traverse_checking(node1.less,less1,node2.greater,greater2)
- traverse_checking(node1.greater,greater1,node2.less,less2)
- traverse_checking(node1.greater,greater1,node2.greater,greater2)
-
- def traverse_no_checking(node1, node2):
- if test_set_visited(node1, node2):
- return
-
- if id(node2)<id(node1):
- # This node pair will be visited in the other order
- #return
- pass
- if isinstance(node1, KDTree.leafnode):
- if isinstance(node2, KDTree.leafnode):
- for i in node1.idx:
- for j in node2.idx:
- if i<j:
- results.add((i,j))
- elif j<i:
- results.add((j,i))
- else:
- traverse_no_checking(node1, node2.less)
- traverse_no_checking(node1, node2.greater)
- else:
- traverse_no_checking(node1.less, node2)
- traverse_no_checking(node1.greater, node2)
-
- traverse_checking(self.tree, Rectangle(self.maxes, self.mins),
- self.tree, Rectangle(self.maxes, self.mins))
- return results
-
-
- def count_neighbors(self, other, r, p=2.):
- """Count how many nearby pairs can be formed.
-
- Count the number of pairs (x1,x2) can be formed, with x1 drawn
- from self and x2 drawn from other, and where distance(x1,x2,p)<=r.
- This is the "two-point correlation" described in Gray and Moore 2000,
- "N-body problems in statistical learning", and the code here is based
- on their algorithm.
-
- Parameters
- ==========
-
- other : KDTree
-
- r : float or one-dimensional array of floats
- The radius to produce a count for. Multiple radii are searched with a single
- tree traversal.
- p : float, 1<=p<=infinity
- Which Minkowski p-norm to use
-
- Returns
- =======
-
- result : integer or one-dimensional array of integers
- The number of pairs. Note that this is internally stored in a numpy int,
- and so may overflow if very large (two billion).
- """
-
- def traverse(node1, rect1, node2, rect2, idx):
- min_r = rect1.min_distance_rectangle(rect2,p)
- max_r = rect1.max_distance_rectangle(rect2,p)
- c_greater = r[idx]>max_r
- result[idx[c_greater]] += node1.children*node2.children
- idx = idx[(min_r<=r[idx]) & (r[idx]<=max_r)]
- if len(idx)==0:
- return
-
- if isinstance(node1,KDTree.leafnode):
- if isinstance(node2,KDTree.leafnode):
- ds = minkowski_distance(self.data[node1.idx][:,np.newaxis,:],
- other.data[node2.idx][np.newaxis,:,:],
- p).ravel()
- ds.sort()
- result[idx] += np.searchsorted(ds,r[idx],side='right')
- else:
- less, greater = rect2.split(node2.split_dim, node2.split)
- traverse(node1, rect1, node2.less, less, idx)
- traverse(node1, rect1, node2.greater, greater, idx)
- else:
- if isinstance(node2,KDTree.leafnode):
- less, greater = rect1.split(node1.split_dim, node1.split)
- traverse(node1.less, less, node2, rect2, idx)
- traverse(node1.greater, greater, node2, rect2, idx)
- else:
- less1, greater1 = rect1.split(node1.split_dim, node1.split)
- less2, greater2 = rect2.split(node2.split_dim, node2.split)
- traverse(node1.less,less1,node2.less,less2,idx)
- traverse(node1.less,less1,node2.greater,greater2,idx)
- traverse(node1.greater,greater1,node2.less,less2,idx)
- traverse(node1.greater,greater1,node2.greater,greater2,idx)
- R1 = Rectangle(self.maxes, self.mins)
- R2 = Rectangle(other.maxes, other.mins)
- if np.shape(r) == ():
- r = np.array([r])
- result = np.zeros(1,dtype=int)
- traverse(self.tree, R1, other.tree, R2, np.arange(1))
- return result[0]
- elif len(np.shape(r))==1:
- r = np.asarray(r)
- n, = r.shape
- result = np.zeros(n,dtype=int)
- traverse(self.tree, R1, other.tree, R2, np.arange(n))
- return result
- else:
- raise ValueError("r must be either a single value or a one-dimensional array of values")
-
- def sparse_distance_matrix(self, other, max_distance, p=2.):
- """Compute a sparse distance matrix
-
- Computes a distance matrix between two KDTrees, leaving as zero
- any distance greater than max_distance.
-
- Parameters
- ==========
-
- other : KDTree
-
- max_distance : positive float
-
- Returns
- =======
-
- result : dok_matrix
- Sparse matrix representing the results in "dictionary of keys" format.
- """
- result = scipy.sparse.dok_matrix((self.n,other.n))
-
- def traverse(node1, rect1, node2, rect2):
- if rect1.min_distance_rectangle(rect2, p)>max_distance:
- return
- elif isinstance(node1, KDTree.leafnode):
- if isinstance(node2, KDTree.leafnode):
- for i in node1.idx:
- for j in node2.idx:
- d = minkowski_distance(self.data[i],other.data[j],p)
- if d<=max_distance:
- result[i,j] = d
- else:
- less, greater = rect2.split(node2.split_dim, node2.split)
- traverse(node1,rect1,node2.less,less)
- traverse(node1,rect1,node2.greater,greater)
- elif isinstance(node2, KDTree.leafnode):
- less, greater = rect1.split(node1.split_dim, node1.split)
- traverse(node1.less,less,node2,rect2)
- traverse(node1.greater,greater,node2,rect2)
- else:
- less1, greater1 = rect1.split(node1.split_dim, node1.split)
- less2, greater2 = rect2.split(node2.split_dim, node2.split)
- traverse(node1.less,less1,node2.less,less2)
- traverse(node1.less,less1,node2.greater,greater2)
- traverse(node1.greater,greater1,node2.less,less2)
- traverse(node1.greater,greater1,node2.greater,greater2)
- traverse(self.tree, Rectangle(self.maxes, self.mins),
- other.tree, Rectangle(other.maxes, other.mins))
-
- return result
-
-
-def distance_matrix(x,y,p=2,threshold=1000000):
- """Compute the distance matrix.
-
- Computes the matrix of all pairwise distances.
-
- Parameters
- ==========
-
- x : array-like, m by k
- y : array-like, n by k
- p : float 1<=p<=infinity
- Which Minkowski p-norm to use.
- threshold : positive integer
- If m*n*k>threshold use a python loop instead of creating
- a very large temporary.
-
- Returns
- =======
-
- result : array-like, m by n
-
-
- """
-
- x = np.asarray(x)
- m, k = x.shape
- y = np.asarray(y)
- n, kk = y.shape
-
- if k != kk:
- raise ValueError("x contains %d-dimensional vectors but y contains %d-dimensional vectors" % (k, kk))
-
- if m*n*k <= threshold:
- return minkowski_distance(x[:,np.newaxis,:],y[np.newaxis,:,:],p)
- else:
- result = np.empty((m,n),dtype=np.float) #FIXME: figure out the best dtype
- if m<n:
- for i in range(m):
- result[i,:] = minkowski_distance(x[i],y,p)
- else:
- for j in range(n):
- result[:,j] = minkowski_distance(x,y[j],p)
- return result
https://bitbucket.org/yt_analysis/yt/commits/d4719f998852/
Changeset: d4719f998852
Branch: yt
User: MatthewTurk
Date: 2013-09-09 14:44:43
Summary: Merged in xarthisius/yt (pull request #588)
Move bundled software to extern directory
Affected #: 40 files
diff -r b3b28b859b0e9d6c7699fb70768dbe08bd8a4d5a -r d4719f9988525d3f4502f33b04786fd948e23d34 yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
+++ b/yt/analysis_modules/halo_merger_tree/enzofof_merger_tree.py
@@ -54,8 +54,8 @@
import glob
from yt.funcs import *
-from yt.utilities.pykdtree import KDTree
-import yt.utilities.pydot as pydot
+from yt.extern.pykdtree import KDTree
+import yt.extern.pydot as pydot
# We don't currently use this, but we may again find a use for it in the
# future.
diff -r b3b28b859b0e9d6c7699fb70768dbe08bd8a4d5a -r d4719f9988525d3f4502f33b04786fd948e23d34 yt/analysis_modules/halo_merger_tree/merger_tree.py
--- a/yt/analysis_modules/halo_merger_tree/merger_tree.py
+++ b/yt/analysis_modules/halo_merger_tree/merger_tree.py
@@ -36,7 +36,7 @@
HaloProfiler
from yt.convenience import load
from yt.utilities.logger import ytLogger as mylog
-import yt.utilities.pydot as pydot
+import yt.extern.pydot as pydot
from yt.utilities.spatial import cKDTree
from yt.utilities.parallel_tools.parallel_analysis_interface import \
ParallelDummy, \
diff -r b3b28b859b0e9d6c7699fb70768dbe08bd8a4d5a -r d4719f9988525d3f4502f33b04786fd948e23d34 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -955,7 +955,7 @@
interpolated using the nearest neighbor method, with *side* points on a
side.
"""
- import yt.utilities.delaunay as de
+ import matplotlib.delaunay.triangulate as de
if log_spacing:
zz = np.log10(self[field])
else:
This diff is so big that we needed to truncate the remainder.
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