[Yt-svn] yt-commit r1279 - in trunk/yt/lagos: . fof hop
sskory at wrangler.dreamhost.com
sskory at wrangler.dreamhost.com
Mon Apr 27 09:21:44 PDT 2009
Author: sskory
Date: Mon Apr 27 09:21:43 2009
New Revision: 1279
URL: http://yt.spacepope.org/changeset/1279
Log:
Changed the FOF kdtree functional names to not conflict with HOP. Also removed a temporary file from hop_regroup.
Modified:
trunk/yt/lagos/HaloFinding.py
trunk/yt/lagos/fof/EnzoFOF.c
trunk/yt/lagos/fof/fof_main.c
trunk/yt/lagos/fof/kd.c
trunk/yt/lagos/fof/kd.h
trunk/yt/lagos/hop/hop_regroup.c
Modified: trunk/yt/lagos/HaloFinding.py
==============================================================================
--- trunk/yt/lagos/HaloFinding.py (original)
+++ trunk/yt/lagos/HaloFinding.py Mon Apr 27 09:21:43 2009
@@ -401,6 +401,7 @@
arr[arr > RE[i]+self.padding] -= dw[i]
def write_out(self, filename):
+ self.data_source.get_data(["particle_velocity_%s" % ax for ax in 'xyz'])
f = self._write_on_root(filename)
HaloList.write_out(self, f)
Modified: trunk/yt/lagos/fof/EnzoFOF.c
==============================================================================
--- trunk/yt/lagos/fof/EnzoFOF.c (original)
+++ trunk/yt/lagos/fof/EnzoFOF.c Mon Apr 27 09:21:43 2009
@@ -84,7 +84,7 @@
/* let's get started with the FOF stuff */
- KD kd;
+ KDFOF kd;
int nBucket,j;
float fPeriod[3],fEps;
int nMembers,nGroup,bVerbose=1;
@@ -101,15 +101,15 @@
/* initialize the kd FOF structure */
- kdInit(&kd,nBucket,fPeriod);
+ kdInitFoF(&kd,nBucket,fPeriod);
- /* kdReadTipsy(kd,stdin,bDark,bGas,bStar); */
+ /* kdReadTipsyFoF(kd,stdin,bDark,bGas,bStar); */
/* Copy positions into kd structure. */
fprintf(stdout, "Filling in %d particles\n", num_particles);
kd->nActive = num_particles;
- kd->p = (PARTICLE *)malloc(kd->nActive*sizeof(PARTICLE));
+ kd->p = (PARTICLEFOF *)malloc(kd->nActive*sizeof(PARTICLEFOF));
assert(kd->p != NULL);
for (i = 0; i < num_particles; i++) {
kd->p[i].iOrder = i;
@@ -118,19 +118,19 @@
kd->p[i].r[2] = (float)(*(npy_float64*) PyArray_GETPTR1(zpos, i));
}
- kdBuildTree(kd);
- kdTime(kd,&sec,&usec);
+ kdBuildTreeFoF(kd);
+ kdTimeFoF(kd,&sec,&usec);
nGroup = kdFoF(kd,fEps);
- kdTime(kd,&sec,&usec);
+ kdTimeFoF(kd,&sec,&usec);
if (bVerbose) printf("Number of initial groups:%d\n",nGroup);
- nGroup = kdTooSmall(kd,nMembers);
+ nGroup = kdTooSmallFoF(kd,nMembers);
if (bVerbose) {
printf("Number of groups:%d\n",nGroup);
printf("FOF CPU TIME: %d.%06d secs\n",sec,usec);
}
- kdOrder(kd);
+ kdOrderFoF(kd);
- /* kdOutGroup(kd,ach); */
+ /* kdOutGroupFoF(kd,ach); */
// Now we need to get the groupID, realID.
// This will give us the index into the original array.
@@ -149,7 +149,7 @@
(npy_int32) kd->p[i].iGroup;
}
- kdFinish(kd);
+ kdFinishFoF(kd);
PyArray_UpdateFlags(particle_group_id, NPY_OWNDATA | particle_group_id->flags);
PyObject *return_value = Py_BuildValue("N", particle_group_id);
Modified: trunk/yt/lagos/fof/fof_main.c
==============================================================================
--- trunk/yt/lagos/fof/fof_main.c (original)
+++ trunk/yt/lagos/fof/fof_main.c Mon Apr 27 09:21:43 2009
@@ -21,7 +21,7 @@
void main(int argc,char **argv)
{
- KD kd;
+ KDFOF kd;
int nBucket,i,j;
char ach[80];
float fPeriod[3],fEps;
@@ -110,22 +110,22 @@
}
else usage();
}
- kdInit(&kd,nBucket,fPeriod);
- kdReadTipsy(kd,stdin,bDark,bGas,bStar);
- kdBuildTree(kd);
- kdTime(kd,&sec,&usec);
+ kdInitFoF(&kd,nBucket,fPeriod);
+ kdReadTipsyFoF(kd,stdin,bDark,bGas,bStar);
+ kdBuildTreeFoF(kd);
+ kdTimeFoF(kd,&sec,&usec);
nGroup = kdFoF(kd,fEps);
- kdTime(kd,&sec,&usec);
+ kdTimeFoF(kd,&sec,&usec);
if (bVerbose) printf("Number of initial groups:%d\n",nGroup);
- nGroup = kdTooSmall(kd,nMembers);
+ nGroup = kdTooSmallFoF(kd,nMembers);
if (bVerbose) {
printf("Number of groups:%d\n",nGroup);
printf("FOF CPU TIME: %d.%06d secs\n",sec,usec);
}
- kdOrder(kd);
+ kdOrderFoF(kd);
strcat(ach,".grp");
- kdOutGroup(kd,ach);
- kdFinish(kd);
+ kdOutGroupFoF(kd,ach);
+ kdFinishFoF(kd);
}
Modified: trunk/yt/lagos/fof/kd.c
==============================================================================
--- trunk/yt/lagos/fof/kd.c (original)
+++ trunk/yt/lagos/fof/kd.c Mon Apr 27 09:21:43 2009
@@ -8,7 +8,7 @@
#include "tipsydefs.h"
-void kdTime(KD kd,int *puSecond,int *puMicro)
+void kdTimeFoF(KDFOF kd,int *puSecond,int *puMicro)
{
struct rusage ru;
@@ -24,12 +24,12 @@
}
-int kdInit(KD *pkd,int nBucket,float *fPeriod)
+int kdInitFoF(KDFOF *pkd,int nBucket,float *fPeriod)
{
- KD kd;
+ KDFOF kd;
int j;
- kd = (KD)malloc(sizeof(struct kdContext));
+ kd = (KDFOF)malloc(sizeof(struct kdContext));
assert(kd != NULL);
kd->nBucket = nBucket;
for (j=0;j<3;++j) kd->fPeriod[j] = fPeriod[j];
@@ -40,7 +40,7 @@
}
-void kdReadTipsy(KD kd,FILE *fp,int bDark,int bGas,int bStar)
+void kdReadTipsyFoF(KDFOF kd,FILE *fp,int bDark,int bGas,int bStar)
{
int i,j,nCnt;
struct dump h;
@@ -64,7 +64,7 @@
/*
** Allocate particles.
*/
- kd->p = (PARTICLE *)malloc(kd->nActive*sizeof(PARTICLE));
+ kd->p = (PARTICLEFOF *)malloc(kd->nActive*sizeof(PARTICLEFOF));
assert(kd->p != NULL);
/*
** Read Stuff!
@@ -97,9 +97,9 @@
}
-void kdSelect(KD kd,int d,int k,int l,int r)
+void kdSelectFoF(KDFOF kd,int d,int k,int l,int r)
{
- PARTICLE *p,t;
+ PARTICLEFOF *p,t;
double v;
int i,j;
@@ -128,7 +128,7 @@
}
-void kdCombine(KDN *p1,KDN *p2,KDN *pOut)
+void kdCombineFoF(KDFOFNFOF *p1,KDFOFNFOF *p2,KDFOFNFOF *pOut)
{
int j;
@@ -148,18 +148,18 @@
}
-void kdUpPass(KD kd,int iCell)
+void kdUpPassFoF(KDFOF kd,int iCell)
{
- KDN *c;
+ KDFOFNFOF *c;
int l,u,pj,j;
c = kd->kdNodes;
if (c[iCell].iDim != -1) {
- l = LOWER(iCell);
- u = UPPER(iCell);
- kdUpPass(kd,l);
- kdUpPass(kd,u);
- kdCombine(&c[l],&c[u],&c[iCell]);
+ l = LOWERFOF(iCell);
+ u = UPPERFOF(iCell);
+ kdUpPassFoF(kd,l);
+ kdUpPassFoF(kd,u);
+ kdCombineFoF(&c[l],&c[u],&c[iCell]);
}
else {
l = c[iCell].pLower;
@@ -179,11 +179,11 @@
}
}
-void kdBuildTree(KD kd)
+void kdBuildTreeFoF(KDFOF kd)
{
int l,n,i,d,m,j,diff;
- KDN *c;
- BND bnd;
+ KDFOFNFOF *c;
+ BNDFOF bnd;
n = kd->nActive;
kd->nLevels = 1;
@@ -196,7 +196,7 @@
kd->nSplit = l;
kd->nNodes = l<<1;
if (kd->kdNodes != NULL) free(kd->kdNodes);
- kd->kdNodes = (KDN *)malloc(kd->nNodes*sizeof(KDN));
+ kd->kdNodes = (KDFOFNFOF *)malloc(kd->nNodes*sizeof(KDFOFNFOF));
assert(kd->kdNodes != NULL);
/*
** Calculate Bounds.
@@ -214,13 +214,13 @@
}
}
/*
- ** Set up ROOT node
+ ** Set up ROOTFOF node
*/
c = kd->kdNodes;
- c[ROOT].pLower = 0;
- c[ROOT].pUpper = kd->nActive-1;
- c[ROOT].bnd = bnd;
- i = ROOT;
+ c[ROOTFOF].pLower = 0;
+ c[ROOTFOF].pUpper = kd->nActive-1;
+ c[ROOTFOF].bnd = bnd;
+ i = ROOTFOF;
while (1) {
assert(c[i].pUpper - c[i].pLower + 1 > 0);
if (i < kd->nSplit && (c[i].pUpper - c[i].pLower) > 0) {
@@ -232,35 +232,35 @@
c[i].iDim = d;
m = (c[i].pLower + c[i].pUpper)/2;
- kdSelect(kd,d,m,c[i].pLower,c[i].pUpper);
+ kdSelectFoF(kd,d,m,c[i].pLower,c[i].pUpper);
c[i].fSplit = kd->p[m].r[d];
- c[LOWER(i)].bnd = c[i].bnd;
- c[LOWER(i)].bnd.fMax[d] = c[i].fSplit;
- c[LOWER(i)].pLower = c[i].pLower;
- c[LOWER(i)].pUpper = m;
- c[UPPER(i)].bnd = c[i].bnd;
- c[UPPER(i)].bnd.fMin[d] = c[i].fSplit;
- c[UPPER(i)].pLower = m+1;
- c[UPPER(i)].pUpper = c[i].pUpper;
+ c[LOWERFOF(i)].bnd = c[i].bnd;
+ c[LOWERFOF(i)].bnd.fMax[d] = c[i].fSplit;
+ c[LOWERFOF(i)].pLower = c[i].pLower;
+ c[LOWERFOF(i)].pUpper = m;
+ c[UPPERFOF(i)].bnd = c[i].bnd;
+ c[UPPERFOF(i)].bnd.fMin[d] = c[i].fSplit;
+ c[UPPERFOF(i)].pLower = m+1;
+ c[UPPERFOF(i)].pUpper = c[i].pUpper;
diff = (m-c[i].pLower+1)-(c[i].pUpper-m);
assert(diff == 0 || diff == 1);
- i = LOWER(i);
+ i = LOWERFOF(i);
}
else {
c[i].iDim = -1;
- SETNEXT(i);
- if (i == ROOT) break;
+ SETNEXTFOF(i);
+ if (i == ROOTFOF) break;
}
}
- kdUpPass(kd,ROOT);
+ kdUpPassFoF(kd,ROOTFOF);
}
-int kdFoF(KD kd,float fEps)
+int kdFoF(KDFOF kd,float fEps)
{
- PARTICLE *p;
- KDN *c;
+ PARTICLEFOF *p;
+ KDFOFNFOF *c;
int pi,pj,pn,cp;
int iGroup;
@@ -300,14 +300,14 @@
x = p[pi].r[0];
y = p[pi].r[1];
z = p[pi].r[2];
- cp = ROOT;
+ cp = ROOTFOF;
while (1) {
- INTERSECT(c,cp,fEps2,lx,ly,lz,x,y,z,sx,sy,sz);
+ INTERSECTFOF(c,cp,fEps2,lx,ly,lz,x,y,z,sx,sy,sz);
/*
** We have an intersection to test.
*/
if (c[cp].iDim >= 0) {
- cp = LOWER(cp);
+ cp = LOWERFOF(cp);
continue;
}
else {
@@ -326,8 +326,8 @@
if (iTail == nFifo) iTail = 0;
}
}
- SETNEXT(cp);
- if (cp == ROOT) break;
+ SETNEXTFOF(cp);
+ if (cp == ROOTFOF) break;
continue;
}
ContainedCell:
@@ -341,8 +341,8 @@
if (iTail == nFifo) iTail = 0;
}
GetNextCell:
- SETNEXT(cp);
- if (cp == ROOT) break;
+ SETNEXTFOF(cp);
+ if (cp == ROOTFOF) break;
}
}
}
@@ -352,7 +352,7 @@
}
-int kdTooSmall(KD kd,int nMembers)
+int kdTooSmallFoF(KDFOF kd,int nMembers)
{
int *pnMembers,*pMap;
int i,pi,nGroup;
@@ -397,20 +397,20 @@
}
-int CmpParticles(const void *v1,const void *v2)
+int CmpParticlesFoF(const void *v1,const void *v2)
{
- PARTICLE *p1 = (PARTICLE *)v1;
- PARTICLE *p2 = (PARTICLE *)v2;
+ PARTICLEFOF *p1 = (PARTICLEFOF *)v1;
+ PARTICLEFOF *p2 = (PARTICLEFOF *)v2;
return(p1->iOrder - p2->iOrder);
}
-void kdOrder(KD kd)
+void kdOrderFoF(KDFOF kd)
{
- qsort(kd->p,kd->nActive,sizeof(PARTICLE),CmpParticles);
+ qsort(kd->p,kd->nActive,sizeof(PARTICLEFOF),CmpParticlesFoF);
}
-void kdOutGroup(KD kd,char *pszFile)
+void kdOutGroupFoF(KDFOF kd,char *pszFile)
{
FILE *fp;
int i,iCnt;
@@ -435,7 +435,7 @@
}
-void kdFinish(KD kd)
+void kdFinishFoF(KDFOF kd)
{
free(kd->p);
free(kd->kdNodes);
Modified: trunk/yt/lagos/fof/kd.h
==============================================================================
--- trunk/yt/lagos/fof/kd.h (original)
+++ trunk/yt/lagos/fof/kd.h Mon Apr 27 09:21:43 2009
@@ -1,41 +1,41 @@
-#ifndef KD_HINCLUDED
-#define KD_HINCLUDED
+#ifndef KDFOF_HINCLUDED
+#define KDFOF_HINCLUDED
-#define ROOT 1
-#define LOWER(i) (i<<1)
-#define UPPER(i) ((i<<1)+1)
-#define PARENT(i) (i>>1)
-#define SIBLING(i) ((i&1)?i-1:i+1)
-#define SETNEXT(i)\
+#define ROOTFOF 1
+#define LOWERFOF(i) (i<<1)
+#define UPPERFOF(i) ((i<<1)+1)
+#define PARENTFOF(i) (i>>1)
+#define SIBLINGFOF(i) ((i&1)?i-1:i+1)
+#define SETNEXTFOF(i)\
{\
while (i&1) i=i>>1;\
++i;\
}
-#define DARK 1
-#define GAS 2
-#define STAR 4
+#define DARKFOF 1
+#define GASFOF 2
+#define STARFOF 4
-#define KD_ORDERTEMP 256
+#define KDFOF_ORDERTEMP 256
typedef struct Particle {
float r[3];
int iGroup;
int iOrder;
- } PARTICLE;
+ } PARTICLEFOF;
typedef struct bndBound {
float fMin[3];
float fMax[3];
- } BND;
+ } BNDFOF;
typedef struct kdNode {
float fSplit;
- BND bnd;
+ BNDFOF bnd;
int iDim;
int pLower;
int pUpper;
- } KDN;
+ } KDNFOF;
typedef struct kdContext {
int nBucket;
@@ -52,15 +52,15 @@
int nLevels;
int nNodes;
int nSplit;
- PARTICLE *p;
- KDN *kdNodes;
+ PARTICLEFOF *p;
+ KDFOFNFOF *kdNodes;
int nGroup;
int uSecond;
int uMicro;
- } * KD;
+ } * KDFOF;
-#define INTERSECT(c,cp,fBall2,lx,ly,lz,x,y,z,sx,sy,sz)\
+#define INTERSECTFOF(c,cp,fBall2,lx,ly,lz,x,y,z,sx,sy,sz)\
{\
float dx,dy,dz,dx1,dy1,dz1,fDist2,fMax2;\
dx = c[cp].bnd.fMin[0]-x;\
@@ -179,15 +179,15 @@
}
-void kdTime(KD,int *,int *);
-int kdInit(KD *,int,float *);
-void kdReadTipsy(KD,FILE *,int,int,int);
-void kdBuildTree(KD);
-int kdFoF(KD,float);
-int kdTooSmall(KD,int);
-void kdOrder(KD);
-void kdOutGroup(KD,char *);
-void kdFinish(KD);
+void kdTimeFoF(KDFOF,int *,int *);
+int kdInitFoF(KDFOF *,int,float *);
+void kdReadTipsyFoF(KDFOF,FILE *,int,int,int);
+void kdBuildTreeFoF(KDFOF);
+int kdFoF(KDFOF,float);
+int kdTooSmallFoF(KDFOF,int);
+void kdOrderFoF(KDFOF);
+void kdOutGroupFoF(KDFOF,char *);
+void kdFinishFoF(KDFOF);
#endif
Modified: trunk/yt/lagos/hop/hop_regroup.c
==============================================================================
--- trunk/yt/lagos/hop/hop_regroup.c (original)
+++ trunk/yt/lagos/hop/hop_regroup.c Mon Apr 27 09:21:43 2009
@@ -441,24 +441,27 @@
below, and if the boundary density is higher than any seen previously
for the lower density group, then record this information */
/* If neither group is above peakdensthresh, skip the boundary */
-
- time_t t1;
- long rand_int;
- (void) time(&t1);
- srand48((long) t1); /* use time in seconds to set seed */
- rand_int = lrand48();
- sprintf(tempfilename,"%032d",rand_int);
- if ((boundfp=fopen(tempfilename,"w"))==NULL)
- myerror("Error opening scratch file");
-
+
+ /* make few arrays to eliminate the need to write a file to disk. The entries in
+ the arrays should be no larger than my_comm->nb.
+ Skory.
+ */
+ int g1temp[my_comm->nb], g2temp[my_comm->nb];
+ float denstemp[my_comm->nb];
+
+ int temppos = 0;
for(j=0;j<(my_comm->nb);j++) {
g1 = my_comm->g1vec[j];
g2 = my_comm->g2vec[j];
dens = my_comm->fdensity[j];
if (gdensity[g1]<peakdensthresh && gdensity[g2]<peakdensthresh) {
if (gdensity[g1]>densthresh && gdensity[g2]>densthresh &&
- dens>densthresh)
- fprintf(boundfp,"%"ISYM" %"ISYM" %"FSYM"\n", g1, g2, dens);
+ dens>densthresh) {
+ g1temp[temppos] = g1;
+ g2temp[temppos] = g2;
+ denstemp[temppos] = dens;
+ temppos += 1;
+ }
continue; /* group isn't dense enough */
}
if (gdensity[g1]>=peakdensthresh && gdensity[g2]>=peakdensthresh)
@@ -485,7 +488,6 @@
densestboundgroup[g2] = g1;
}
} /* Get the next boundary line */
- fclose(boundfp);
/* Now the fringe groups are connected to the proper group
@@ -495,12 +497,11 @@
/* Keep the density of the connection in densestbound, and the
proper group it leads to in densestboundgroup */
do {
- if ((boundfp=fopen(tempfilename,"r"))==NULL)
- myerror("Error opening scratch file for read");
changes = 0;
- while (fgets(line,80,boundfp)!=NULL) {
- if (sscanf(line,"%"ISYM" %"ISYM" %"FSYM, &g1, &g2, &dens)!=3)
- myerror("Error reading boundary.");
+ for (j=0;j<temppos;j++) {
+ g1 = g1temp[j];
+ g2 = g2temp[j];
+ dens = denstemp[j];
/* If the density of this boundary and the densestbound of
the other group is higher than a group's densestbound, then
replace it. */
@@ -515,7 +516,6 @@
densestboundgroup[g2] = densestboundgroup[g1];
}
}
- fclose(boundfp);
} while (changes);
/* Now connect the low-density groups to their densest boundaries */
More information about the yt-svn
mailing list