[Yt-svn] yt-commit r548 - in trunk/yt/lagos: . hop
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Tue Jun 10 09:05:14 PDT 2008
Author: mturk
Date: Tue Jun 10 09:05:10 2008
New Revision: 548
URL: http://yt.spacepope.org/changeset/548
Log:
Imported HOP code from the original authors with modifications by Stephen Skory
and myself. Right now it only returns density and tags, but the plan is to
have it return recarrays and have output to SQLite be an option.
I've tested it for memory leaks and found nothing significant, but I need to do
more testing. Additionally, while I have performed very simple tests of
correctness, I have not yet completed the full tests. So, use at your own
risk.
Added:
trunk/yt/lagos/hop/
trunk/yt/lagos/hop/EnzoHop.c
trunk/yt/lagos/hop/README
trunk/yt/lagos/hop/__init__.py
trunk/yt/lagos/hop/hop.h
trunk/yt/lagos/hop/hop_hop.c
trunk/yt/lagos/hop/hop_kd.c
trunk/yt/lagos/hop/hop_regroup.c
trunk/yt/lagos/hop/hop_slice.c
trunk/yt/lagos/hop/hop_smooth.c
trunk/yt/lagos/hop/kd.h
trunk/yt/lagos/hop/setup.py
trunk/yt/lagos/hop/slice.h
trunk/yt/lagos/hop/smooth.h
Modified:
trunk/yt/lagos/setup.py
Added: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/EnzoHop.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,225 @@
+/************************************************************************
+* Copyright (C) 2008 Matthew Turk. All Rights Reserved.
+*
+* This file is part of yt.
+*
+* yt is free software; you can redistribute it and/or modify
+* it under the terms of the GNU General Public License as published by
+* the Free Software Foundation; either version 3 of the License, or
+* (at your option) any later version.
+*
+* This program 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 General Public License for more details.
+*
+* You should have received a copy of the GNU General Public License
+* along with this program. If not, see <http://www.gnu.org/licenses/>.
+*
+************************************************************************/
+
+//
+// EnzoHop
+// A module for running HOP halo finding on a set of particles
+//
+
+#include "Python.h"
+#include <stdio.h>
+#include <math.h>
+#include <signal.h>
+#include <ctype.h>
+#include "kd.h"
+#include "hop.h"
+
+#include "numpy/ndarrayobject.h"
+
+void initgrouplist(Grouplist *g);
+void hop_main(KD kd, HC *my_comm, float densthres);
+void regroup_main(float dens_outer, HC *my_comm);
+static PyObject *_HOPerror;
+
+static PyObject *
+Py_EnzoHop(PyObject *obj, PyObject *args)
+{
+ PyObject *oxpos, *oypos, *ozpos,
+ *omass, *oID;
+
+ PyArrayObject *xpos, *ypos, *zpos,
+ *mass, *ID;
+ xpos=ypos=zpos=mass=ID=NULL;
+ npy_float64 totalmass = 0;
+
+ int i;
+
+ if (!PyArg_ParseTuple(args, "OOOOO",
+ &oxpos, &oypos, &ozpos, &omass, &oID))
+ return PyErr_Format(_HOPerror,
+ "EnzoHop: Invalid parameters.");
+
+ /* First the regular source arrays */
+
+ xpos = (PyArrayObject *) PyArray_FromAny(oxpos,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if(!xpos){
+ PyErr_Format(_HOPerror,
+ "EnzoHop: xpos didn't work.");
+ goto _fail;
+ }
+ int num_particles = PyArray_SIZE(xpos);
+
+ ypos = (PyArrayObject *) PyArray_FromAny(oypos,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((!ypos)||(PyArray_SIZE(ypos) != num_particles)) {
+ PyErr_Format(_HOPerror,
+ "EnzoHop: xpos and ypos must be the same length.");
+ goto _fail;
+ }
+
+ zpos = (PyArrayObject *) PyArray_FromAny(ozpos,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((!zpos)||(PyArray_SIZE(zpos) != num_particles)) {
+ PyErr_Format(_HOPerror,
+ "EnzoHop: xpos and zpos must be the same length.");
+ goto _fail;
+ }
+
+ mass = (PyArrayObject *) PyArray_FromAny(omass,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((!zpos)||(PyArray_SIZE(mass) != num_particles)) {
+ PyErr_Format(_HOPerror,
+ "EnzoHop: xpos and mass must be the same length.");
+ goto _fail;
+ }
+
+ ID = (PyArrayObject *) PyArray_FromAny(oID,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((!ID)||(PyArray_SIZE(ID) != num_particles)) {
+ PyErr_Format(_HOPerror,
+ "EnzoHop: xpos and ID must be the same length.");
+ goto _fail;
+ }
+
+ for(i = 0; i < num_particles; i++)
+ totalmass+=*(npy_float64*)PyArray_GETPTR1(mass,i);
+
+ /* initialize the kd hop structure */
+
+ KD kd;
+ int nBucket = 16, kdcount = 0;
+ kdInit(&kd, nBucket);
+ kd->nActive = num_particles;
+ kd->p = malloc(sizeof(PARTICLE)*num_particles);
+ if (kd->p == NULL) {
+ fprintf(stderr, "failed allocating particles.\n");
+ goto _fail;
+ }
+
+ /* Copy positions into kd structure. */
+
+ for (i = 0; i < num_particles; i++) {
+ kd->p[i].r[0] = (float)*(npy_float64*) PyArray_GETPTR1(xpos, i);
+ kd->p[i].r[1] = (float)*(npy_float64*) PyArray_GETPTR1(ypos, i);
+ kd->p[i].r[2] = (float)*(npy_float64*) PyArray_GETPTR1(zpos, i);
+ kd->p[i].fMass = (float)(*(npy_float64*) PyArray_GETPTR1(mass, i)/totalmass);
+ kd->p[i].iID = (int)*(npy_int64*) PyArray_GETPTR1(ID,i); /* S. Skory */
+ }
+
+ HC my_comm;
+ my_comm.s = newslice();
+ my_comm.gl = (Grouplist*)malloc(sizeof(Grouplist));
+ initgrouplist(my_comm.gl);
+
+ fprintf(stderr, "Calling hop... %d\n",num_particles);
+ float thresh = 80.0;
+ hop_main(kd, &my_comm, thresh);
+
+ fprintf(stderr, "Calling regroup...\n");
+ regroup_main(thresh, &my_comm);
+
+ // Now we need to get the groupID, realID and the density.
+ // This will give us the index into the original array.
+ // Additionally, note that we don't really need to tie the index
+ // back to the ID in this code, as we can do that back in the python code.
+ // All we need to do is provide density and group information.
+
+ // Tags (as per writetagsf77) are in gl.s->ntag+1 and there are gl.s->numlist of them.
+ PyArrayObject *particle_density = (PyArrayObject *)
+ PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
+ PyArray_DescrFromType(NPY_FLOAT64));
+ PyArrayObject *particle_group_id = (PyArrayObject *)
+ PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
+ PyArray_DescrFromType(NPY_INT64));
+
+ for (i = 0; i < num_particles; i++) {
+ // Density is in kd->p[i].fDensity
+ *(npy_float64*)(PyArray_GETPTR1(particle_density, i)) =
+ (npy_float64) kd->p[i].fDensity;
+ // tag is in gl.s->ntag[i+1]
+ *(npy_int64*)(PyArray_GETPTR1(particle_group_id, i)) =
+ (npy_int64) my_comm.s->ntag[i+1];
+ }
+
+ kdFinish(kd);
+ free(my_comm.gl);
+ free_slice(my_comm.s);
+
+ PyArray_UpdateFlags(particle_density, NPY_OWNDATA | particle_density->flags);
+ PyArray_UpdateFlags(particle_group_id, NPY_OWNDATA | particle_group_id->flags);
+ PyObject *return_value = Py_BuildValue("NN", particle_density, particle_group_id);
+
+ Py_DECREF(xpos);
+ Py_DECREF(ypos);
+ Py_DECREF(zpos);
+ Py_DECREF(mass);
+ Py_DECREF(ID);
+
+ /* We don't need this, as it's done in kdFinish
+ if(kd->p!=NULL)free(kd->p);
+ */
+
+ return return_value;
+
+_fail:
+ Py_XDECREF(xpos);
+ Py_XDECREF(ypos);
+ Py_XDECREF(zpos);
+ Py_XDECREF(mass);
+ Py_XDECREF(ID);
+
+ if(kd->p!=NULL)free(kd->p);
+
+ return NULL;
+
+}
+
+static PyMethodDef _HOPMethods[] = {
+ {"RunHOP", Py_EnzoHop, METH_VARARGS},
+ {NULL, NULL} /* Sentinel */
+};
+
+/* platform independent*/
+#ifdef MS_WIN32
+__declspec(dllexport)
+#endif
+
+void initEnzoHop(void)
+{
+ PyObject *m, *d;
+ m = Py_InitModule("EnzoHop", _HOPMethods);
+ d = PyModule_GetDict(m);
+ _HOPerror = PyErr_NewException("EnzoHop.HOPerror", NULL, NULL);
+ PyDict_SetItemString(d, "error", _HOPerror);
+ import_array();
+}
+
+/*
+ * Local Variables:
+ * mode: C
+ * c-file-style: "python"
+ * End:
+ */
Added: trunk/yt/lagos/hop/README
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/README Tue Jun 10 09:05:10 2008
@@ -0,0 +1,22 @@
+Matthew Turk
+May 2008
+
+This code (described below) has been modified to be wrapped as a shared library
+callable from Python, as a part of the yt toolkit.
+
+Stephen Skory
+May/June 2007
+
+This is a new implementation of hop for enzo datasets, to replace the
+fragile 'enzohop.' enzohop uses the enzo grid functionality to extract
+the particle data from the HDF5 datasets. newhop uses plain HDF5 C++ calls
+to extract the data, which is then fed into the hop mechanism. As far as I
+know, this version is fine with 64 bit integers/floats, which enzohop isn't.
+
+There are a few versions of newhop which build on datastar just fine. I
+haven't tested it on other machines. The default build 'newhop' is for
+packed datasets and will include both stars and dm in the grouping.
+THe other versions are for non-packed datasets, or if you want to only
+consider dm particles for grouping. Hop doesn't like datasets with too
+many particles, (I've never done tests, but I know that 20 million
+particles give hop problems) so sometimes dm-only is the only way to go.
Added: trunk/yt/lagos/hop/__init__.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/__init__.py Tue Jun 10 09:05:10 2008
@@ -0,0 +1 @@
+from EnzoHop import *
Added: trunk/yt/lagos/hop/hop.h
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop.h Tue Jun 10 09:05:10 2008
@@ -0,0 +1,35 @@
+#include "slice.h"
+
+/* ----------------------------------------------------------------------- */
+/* The following structures track all the information about the groups */
+
+typedef struct groupstruct {
+ int npart; /* Number of particles in the group */
+ int npartcum; /* Cumulative number of particles */
+ int nread; /* Number read so far, also a utility field */
+ double compos[3], comvel[3];/* Lists of group CoM position and velocities */
+ double comtemp[3]; /* Temporary CoM position */
+ int idmerge; /* The new group ID # after merging. Not necessarily
+ unique! */
+ int rootgroup; /* The fully traced group id */
+} Group; /* Type Group is defined */
+
+typedef struct groupliststruct {
+ int npart; /* Number of particles in the simulation */
+ int ngroups; /* Number of groups in list */
+ int nnewgroups; /* Number of groups after relabeling */
+ int npartingroups; /* Number of particles in groups */
+ Group *list; /* List of groups, zero-offset */
+} Grouplist; /* Type Grouplist is defined */
+
+
+typedef struct hopComm {
+ int ngroups;
+ int nb;
+ float *gdensity;
+ float *g1vec;
+ float *g2vec;
+ float *fdensity;
+ Grouplist *gl;
+ Slice *s;
+} HC;
Added: trunk/yt/lagos/hop/hop_hop.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop_hop.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,963 @@
+/* HOP.C -- Daniel Eisenstein, 1997
+Based on a paper by Daniel Eisenstein & Piet Hut,
+"HOP: A New Group-Finding Algorithm for N-body Simulations."
+See the included documentation or view it at
+http://www.sns.ias.edu/~eisenste/hop/hop_doc.html */
+
+/* main() was customized from that of SMOOTH v2.0.1, by Joachim Stadel
+ and the NASA HPCC ESS at the University of Washington Dept. of Astronomy. */
+/* PrepareKD() is a code fragment from the same program. */
+/* ssort() is a C-translation of the Slatec FORTRAN routine of
+ R.E. Jones and J.A. Wisniewski (SNLA) */
+/* The other routines were written by DJE. */
+
+/* Version 1.0 (12/15/97) -- Original Release */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <ctype.h>
+#include <assert.h>
+#include "kd.h"
+#include "hop.h"
+#include "smooth.h"
+
+//#include "macros_and_parameters.h"
+#define ISYM "d"
+#define GSYM "g"
+#define FSYM "f"
+
+/* To give info to the user: INFORM("info"); */
+#define INFORM(string) printf(string); fflush(stdout)
+
+int SnapNumber;
+
+int ReadSimulationFile(KD, FILE *);
+
+void smDensityTH(SMX smx,int pi,int nSmooth,int *pList,float *fList);
+
+void smHop(SMX smx,int pi,int nSmooth,int *pList,float *fList);
+void FindGroups(SMX smx);
+void SortGroups(SMX smx);
+
+void MergeGroupsHash(SMX smx);
+void smMergeHash(SMX smx,int pi,int nSmooth,int *pList,float *fList);
+void ReSizeSMX(SMX smx, int nSmooth);
+
+void PrepareKD(KD kd);
+void binOutHop(SMX smx, HC *my_comm, float densthres);
+void binOutDensity(SMX smx, HC *my_comm, float densthres);
+void binInDensity(SMX smx, FILE *fp);
+void outGroupMerge(SMX smx, HC *my_comm);
+
+/* void main(int argc,char **argv) */
+void hop_main(KD kd, HC *my_comm, float densthres)
+{
+ /* KD kd; */
+ SMX smx;
+ int nBucket,nSmooth,i,j;
+ FILE *fp, *fpb;
+ char ach[80],achFile[80], *inputfile, *densfile;
+ float fPeriod[3];
+ int bDensity,bGroup,bSym,bMerge,nDens,nHop,nMerge,bTopHat;
+ float fDensThresh;
+
+ nBucket = 16;
+ nSmooth = 64;
+ nDens = 64;
+ nHop = -1;
+ fDensThresh = -1.0;
+ bDensity = 3;
+ bGroup = 3;
+ bMerge = 3;
+ bSym = 1;
+ bTopHat = 0;
+ strcpy(achFile,"output_hop");
+ inputfile = NULL;
+ i = 1;
+/* for (j=0;j<3;++j) fPeriod[j] = HUGE; */
+ for (j=0;j<3;++j) fPeriod[j] = 1.0;
+ nMerge = 4;
+
+
+ if (nHop<0) nHop=nDens;
+ if (bDensity==0) nSmooth = nHop+1;
+ else nSmooth = nDens+1;
+ /* When smSmooth() is asked for nSmooth particles, it seems to
+ generally return nSmooth-1 particles, including primary itself.
+ Hence, when we want nDens or nSmooth particles (including the
+ primary) we should ask for one more. By convention, I've chosen
+ to have nMerge reflect the number *not* including the primary,
+ so in this case we need to ask for two more! */
+
+ assert(!bMerge || nMerge<nHop-1);
+ /* This is generally the case of interest, and there's an
+ optimization to be had. I haven't the coded the other case */
+ assert(!bMerge || bGroup);
+ /* Can't run Merging without Grouping */
+ if (!bMerge) nMerge = 9999999;
+ /* If we're not Merging, we don't want smHop() to set
+ up for Merging. This kludge avoids that */
+
+ PrepareKD(kd);
+
+ smInit(&smx,kd,nSmooth,fPeriod);
+
+ smx->nHop = nHop;
+ smx->nDens = nDens;
+ smx->nMerge = nMerge;
+ smx->nGroups = 0;
+ smx->fDensThresh = fDensThresh;
+
+ INFORM("Building Tree...\n");
+ kdBuildTree(kd);
+
+ if (bDensity) {
+ INFORM("Finding Densities...\n");
+ if (bTopHat) smSmooth(smx,smDensityTH);
+ else if (bSym) smSmooth(smx,smDensitySym);
+ else smSmooth(smx,smDensity);
+ } /* Else, we've read them */
+ if (bGroup) {
+ INFORM("Finding Densest Neighbors...\n");
+ if (bDensity && nHop<nSmooth) smReSmooth(smx,smHop);
+ else {
+ if (nHop>=nSmooth) {
+ nSmooth = nHop+1;
+ ReSizeSMX(smx,nSmooth);
+ }
+ smSmooth(smx,smHop);
+ }
+ }
+
+ INFORM("Grouping...\n");
+ if (bGroup) FindGroups(smx);
+ if (bGroup) SortGroups(smx);
+
+ if (bMerge) {
+ INFORM("Merging Groups...\n");
+ MergeGroupsHash(smx);
+ }
+
+ kdOrder(kd);
+ INFORM("Writing Output...\n");
+
+ if (bMerge&2) {
+ smx->nSmooth=nSmooth; /* Restore this for output */
+ outGroupMerge(smx, my_comm);
+ }
+ if (bMerge) free(smx->hash);
+
+ if (bGroup&2) {
+ binOutHop(smx, my_comm, densthres);
+ }
+ if (bGroup) {free(smx->densestingroup); free(smx->nmembers);}
+ smFinish(smx);
+ //kdFinish(kd);
+ INFORM("All Done!");
+ return;
+}
+
+
+/* ============================================================= */
+/* ===================== New Density Routine =================== */
+/* ============================================================= */
+
+void smDensityTH(SMX smx,int pi,int nSmooth,int *pList,float *fList)
+/* Find density only using top-hat kernal. */
+{
+#ifdef DIFFERENT_MASSES
+ int j;
+ float totalmass;
+ for (j=0,totalmass=0.0; j<nSmooth; j++)
+ totalmass += smx->kd->p[pList[j]].fMass;
+ smx->kd->p[pi].fDensity = totalmass*0.75*M_1_PI/
+ smx->pfBall2[pi]/sqrt(smx->pfBall2[pi]);
+#else
+ /* This case is simple: the total mass is nSmooth times the mass
+ per particle */
+ smx->kd->p[pi].fDensity = (nSmooth)*smx->kd->fMass*0.75*M_1_PI/
+ smx->pfBall2[pi]/sqrt(smx->pfBall2[pi]);
+#endif
+ return;
+}
+
+/* ============================================================= */
+/* ================== Hop to Neighbors/Form Groups ============= */
+/* ============================================================= */
+
+void smHop(SMX smx,int pi,int nSmooth,int *pList,float *fList)
+/* Look at the nHop nearest particles and find the one with the
+highest density. Store its ID number in iHop as -1-ID (to make it negative) */
+/* nSmooth tends to be the expected value (smx->nSmooth-1) but can vary plus
+or minus 1 (at least). */
+/* If Merging is turned off, nMerge should be huge so as to avoid extra
+sorting below */
+{
+ int i,max, search, didsort;
+ float maxden;
+ void ssort(float X[], int Y[], int N, int KFLAG);
+
+ /* If the density is less than the threshold requirement, then assign 0 */
+ if (smx->kd->p[pi].fDensity<smx->fDensThresh) {
+ smx->kd->p[pi].iHop = 0;
+ return;
+ }
+
+ /* If we have exactly the right number, then it doesn't matter how
+ we search them. Otherwise, we need to sort first. */
+ /* We can destroy pList and fList if we want. fList holds the square radii*/
+
+ search = smx->nHop;
+ if (smx->nHop>nSmooth) search = nSmooth; /* Very rare exception */
+
+ if (smx->nHop<nSmooth || smx->nMerge+2<nSmooth) {
+ ssort(fList-1, pList-1, nSmooth, 2);
+ didsort = 1;
+ } else didsort = 0;
+
+ max = 0;
+ maxden = 0.0;
+ for (i=0;i<search;++i) {
+ if (smx->kd->p[pList[i]].fDensity>maxden) {
+ max = i;
+ maxden = smx->kd->p[pList[i]].fDensity;
+ }
+ }
+ smx->kd->p[pi].iHop = -1-pList[max];
+
+ /* If a sort was done, then we can save time in the Merge step by
+ recording the new Ball radius. */
+ /* Place the new radius in between the two boundary because occasionally
+ the floating point math comes out strange when comparing two floats */
+ if (didsort && smx->nMerge+2<nSmooth)
+ smx->pfBall2[pi] = 0.5*(fList[smx->nMerge+1]+fList[smx->nMerge]);
+ return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void FindGroups(SMX smx)
+/* Number the maxima. Trace each particle uphill to a maximum. */
+/* The local maxima were stored as in iHop as -1-ID (negative numbers);
+now we will store group number as positive numbers (1..n) in the same spot */
+/* Zero is left as an error condition */
+/* The particles MUST be left in the BuildTree order, for that is how the
+iHop tracing is done */
+/* Allocate space for densestingroup, from 0 to nGroups
+(inclusive) and store the particle number of the maxima, which is the
+densest particle in the group. Ditto for nmembers[], the number of
+particles in the group */
+{
+ int j, ng, current, localmax;
+ PARTICLE *p;
+
+ smx->nGroups = 0;
+ /* First look for maxima, where particle ID = iHop. Number the groups */
+ for (j=0, p=smx->kd->p;j<smx->kd->nActive;j++,p++)
+ if (p->iHop == -1-j) { /* Was p->iOrder */
+ /* Yes, it's a maximum */
+ smx->nGroups++;
+ /* p->iHop = smx->nGroups; */
+ }
+
+ /* Now catalog the maxima, before numbering the groups */
+ smx->densestingroup = (int *)malloc((size_t)((smx->nGroups+1)*sizeof(int)));
+ assert(smx->densestingroup!=NULL);
+ smx->nmembers = (int *)malloc((size_t)(sizeof(int)*(smx->nGroups+1)));
+ assert(smx->nmembers!=NULL);
+
+ ng = 0;
+ for (j=0,p=smx->kd->p;j<smx->kd->nActive;j++,p++)
+ if (p->iHop== -1-j) {
+ /* It's a maximum */
+ ng++;
+ smx->densestingroup[ng] = p->iOrder;
+ p->iHop = ng;
+ }
+
+ /* Now take the remaining particles and trace up to a maximum */
+ for (j=0,p=smx->kd->p;j<smx->kd->nActive;j++,p++) {
+ if (p->iHop>=0) continue; /* It's a maximum or an error */
+ localmax = -1-p->iHop;
+ while (smx->kd->p[localmax].iHop<0)
+ localmax = -1-smx->kd->p[localmax].iHop;
+ ng = smx->kd->p[localmax].iHop;
+ /* Now assign this group number to the whole lineage! */
+ /* Note that errors (iHop=0) will propagate */
+ current = -1-p->iHop;
+ p->iHop = ng;
+ while (smx->kd->p[current].iHop<0) {
+ localmax = -1-smx->kd->p[current].iHop;
+ smx->kd->p[current].iHop = ng;
+ current = localmax;
+ }
+ }
+ return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void SortGroups(SMX smx)
+/* Renumber the groups in order of number of members. */
+/* Move the group numbering from unit offset to zero offset */
+/* Move error condition (group=0) to (group=-1) */
+/* Store the number of members and most dense particle in each group */
+{
+ int j, *indx, *irank, *ip;
+ PARTICLE *p;
+ void make_rank_table(int n, int *ivect, int *rank);
+
+ indx = (int *)malloc((size_t)(sizeof(int)*(smx->nGroups+1)));
+ assert(indx!=NULL);
+ irank = (int *)malloc((size_t)(sizeof(int)*(smx->nGroups+1)));
+ assert(irank!=NULL);
+
+ /* Count the number of members in each group */
+ for (j=0;j<=smx->nGroups;j++) smx->nmembers[j]=0;
+ for (j=0,p=smx->kd->p;j<smx->kd->nActive;j++,p++)
+ smx->nmembers[p->iHop]++;
+
+ make_rank_table(smx->nGroups, smx->nmembers,irank);
+
+ for (j=1;j<=smx->nGroups;j++) irank[j] = smx->nGroups-irank[j];
+ irank[0] = -1; /* Move old 0's to new -1's */
+ /* irank[j] is the new group number of group j: zero-offset, ordered
+ large to small */
+
+ /* Relabel all the particles */
+ for (j=0,p=smx->kd->p;j<smx->kd->nActive;j++,p++)
+ p->iHop = irank[p->iHop];
+
+ /* Sort the nmembers and densestingroup lists to reflect the new ordering */
+ /* Use indx as a temp space */
+ for (j=1;j<=smx->nGroups;j++) indx[irank[j]]=smx->densestingroup[j];
+ ip = smx->densestingroup;
+ smx->densestingroup = indx;
+ indx = ip;
+ /* Number of error (old_group=0) is in nmembers[0]; move it to [nGroups] */
+ for (j=1;j<=smx->nGroups;j++) indx[irank[j]]=smx->nmembers[j];
+ indx[smx->nGroups]=smx->nmembers[0];
+ free(smx->nmembers);
+ smx->nmembers = indx;
+
+ free(irank);
+ /* Note that the memory allocated to indx is now used by smx->densestingroup
+ and so it should not be free'd. */
+ return;
+}
+
+/* ================================================================== */
+/* ========================== Group Merging ========================= */
+/* ================================================================== */
+
+void MergeGroupsHash(SMX smx)
+/* We're going to look through the particles looking for boundary particles,
+defined as particles with close neighbors of a different group, and store
+the most dense boundary point (average of the two points) */
+/* The matrix of boundaries is stored in a hash table */
+/* SortGroups() should be called previous to this, so that all the
+particles are in the assumed group numbering, i.e. 0 to ngroup-1, with
+-1 being unattached. The tags are not altered */
+/* In smHop, if nMerge+2 was smaller than nSmooth, we set the new radius
+for searching. If not, we left the old radius alone. Either way, we're
+ready to go. */
+{
+ int j, k, g, next, newgroup;
+
+ ReSizeSMX(smx, smx->nMerge+2); /* Alter the smoothing scale on smx */
+ smx->nHashLength = smx->nGroups*10+1;
+ smx->hash = (Boundary *)malloc(smx->nHashLength*sizeof(Boundary));
+ assert(smx->hash!=NULL);
+ for (j=0;j<smx->nHashLength;j++) {
+ smx->hash[j].nGroup1 = -1;
+ smx->hash[j].nGroup2 = -1;
+ smx->hash[j].fDensity = -1.0;
+ } /* Mark the slot as unused */
+
+ smReSmooth(smx,smMergeHash); /* Record all the boundary particles */
+ return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void smMergeHash(SMX smx,int pi,int nSmooth,int *pList,float *fList)
+/* Look at the list for groups which are not that of the particle */
+/* If found, and if density is high enough, then mark it as a boundary */
+/* by recording it in the hash table */
+{
+ int j,group;
+ float averdensity;
+ int g1,g2,count;
+ unsigned long hashpoint;
+ Boundary *hp;
+ int search;
+ void ssort(float X[], int Y[], int N, int KFLAG);
+
+ group = smx->kd->p[pi].iHop;
+ if (group==(-1)) return; /* This particle isn't in a group */
+
+ /* It seems that BallGather doesn't always get the right number of
+ particles....*/
+ search = nSmooth;
+ if (nSmooth>smx->nMerge+1) {
+ ssort(fList-1,pList-1,nSmooth,2);
+ search = smx->nMerge+1;
+ }
+ for (j=0;j<search;j++) {
+ g2=smx->kd->p[pList[j]].iHop;
+ if (g2==-1 || g2==group) continue; /* Same group or unassigned */
+ /* It's in a different group; we need to connect the two */
+ if (group<g2) g1=group;
+ else {g1=g2; g2=group;}
+ averdensity = 0.5*(smx->kd->p[pi].fDensity +
+ smx->kd->p[pList[j]].fDensity);
+ hashpoint = (g1+1)*g2; /* Avoid multiplying by 0 */
+ hashpoint = hashpoint % smx->nHashLength;
+ hp = smx->hash+hashpoint;
+ count = 0;
+ for (;;) {
+ if (hp->nGroup1==(-1)) { /* Empty slot */
+ hp->nGroup1 = g1;
+ hp->nGroup2 = g2;
+ hp->fDensity = averdensity;
+ break;
+ }
+ if (hp->nGroup1==g1 && hp->nGroup2==g2) {
+ /* We've seen this pair of groups before */
+ if (hp->fDensity > averdensity) break;
+ else {
+ hp->fDensity = averdensity;
+ break;
+ }
+ }
+ /* Else, this slot was full, go to the next one */
+ hp++;
+ if (hp>=smx->hash+smx->nHashLength) hp = smx->hash;
+ if (++count>1000) {
+ fprintf(stderr,"Hash Table is too full.\n");
+ exit(1);
+ }
+ }
+ /* Look at the next particle */
+ }
+ return;
+}
+
+
+/* ----------------------------------------------------------------- */
+
+void ReSizeSMX(SMX smx, int nSmooth)
+/* Set a new smoothing length, resizing the arrays which depend on this,
+but leaving the particle information intact. */
+/* However, because we won't always have resized pfBall2 (the search
+radius) correctly, we won't reduce the size of the fList and pList
+arrays */
+{
+ PQ_STATIC;
+ if (nSmooth>smx->nSmooth) { /* We're increasing the size */
+ smx->nListSize = nSmooth+RESMOOTH_SAFE;
+ free(smx->fList);
+ smx->fList = (float *)malloc(smx->nListSize*sizeof(float));
+ assert(smx->fList != NULL);
+ free(smx->pList);
+ smx->pList = (int *)malloc(smx->nListSize*sizeof(int));
+ assert(smx->pList != NULL);
+ }
+ smx->nSmooth=nSmooth;
+ free(smx->pq);
+ smx->pq = (PQ *)malloc(nSmooth*sizeof(PQ));
+ assert(smx->pq != NULL);
+ PQ_INIT(smx->pq,nSmooth);
+ return;
+}
+
+/* ===================================================================== */
+/* ===================== Input/Output, Binary and ASCII ================ */
+/* ===================================================================== */
+
+void PrepareKD(KD kd)
+/* This labels all the particles and finds the min/max of the positions */
+/* It used to appear in kd.c within kdReadTipsy(), but it seems so general
+that I'll spare the user the trouble of including it in any custom input
+routines */
+{
+ BND bnd;
+ int i, j;
+
+ /* Label the particles, so that we can restore the order at the end */
+ for (i=0;i<kd->nActive;i++) {
+ kd->p[i].iOrder=i;
+ }
+ /*
+ ** Calculate Bounds.
+ */
+ for (j=0;j<3;++j) {
+ bnd.fMin[j] = kd->p[0].r[j];
+ bnd.fMax[j] = kd->p[0].r[j];
+ }
+ for (i=1;i<kd->nActive;++i) {
+ for (j=0;j<3;++j) {
+ if (bnd.fMin[j] > kd->p[i].r[j])
+ bnd.fMin[j] = kd->p[i].r[j];
+ else if (bnd.fMax[j] < kd->p[i].r[j])
+ bnd.fMax[j] = kd->p[i].r[j];
+ }
+ }
+ kd->bnd = bnd;
+ return;
+}
+
+void binOutHop(SMX smx, HC *my_comm, float densthres)
+/* Write Group tag for each particle. Particles should be ordered. */
+/* Binary file: nActive, nGroups, list of Groups */
+{
+ int j,dummy;
+ /* FILE *blahfp = fopen("Part-PreMergeGroup.txt","w"); *//* S Skory */
+ Grouplist *g = my_comm->gl;
+ Slice *s = my_comm->s;
+
+ g->npart = s->numlist = s->numpart = smx->kd->nActive;
+ g->ngroups = smx->nGroups;
+ s->ntag = ivector(1,s->numlist);
+ s->ID = ivector(1,s->numlist);
+ for (j=0;j<smx->kd->nActive;j++) {
+ s->ID[1+j] = smx->kd->p[j].iID; /* S Skory's addition */
+ if (smx->kd->p[j].fDensity < densthres) s->ntag[j+1] = -1;
+ else smx->kd->p[j].iHop;
+
+ }
+
+ /* Here I'm going to add on the end of the file the real particle IDs for all the particles
+ added above, in the same order as above. S Skory */
+ return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void binOutDensity(SMX smx, HC *my_comm, float densthres)
+/* Write the density for each particle. Particles should be ordered. */
+/* Binary file: nActive, list of Densities */
+{
+ int j,dummy;
+ Slice *s = my_comm->s;
+ fprintf(stderr, "binOutDensity %d\n", smx->kd->nActive);
+ for (j=0;j<smx->kd->nActive;j++); return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void binInDensity(SMX smx, FILE *fp)
+/* Read in the densities for each particle. Particles should be ordered. */
+/* Binary file: nActive, list of Densities */
+{
+ int j, dummy;
+
+ if (fread(&dummy,sizeof(int),1,fp)!=1) {
+ fprintf(stderr,"Format of density file seems wrong.\n"); exit(1);
+ }
+ assert(dummy==smx->kd->nActive);
+ for (j=0;j<smx->kd->nActive;j++)
+ if (fread(&(smx->kd->p[j].fDensity),sizeof(float),1,fp)!=1) {
+ fprintf(stderr,"Error reading density file.\n");
+ exit(1);
+ }
+ return;
+}
+
+/* ----------------------------------------------------------------- */
+
+void outGroupMerge(SMX smx, HC *my_comm)
+/* Write an ASCII file with information on the groups and group merging */
+/* Start the boundary list with the only ### line */
+/* Groups should be ordered before calling this (else densities will be wrong)*/
+{
+ int j, den;
+ Boundary *hp;
+
+ my_comm->gdensity = vector(0,smx->nGroups-1);
+ for (j=0;j<smx->nGroups;j++) {
+ my_comm->gdensity[j]=smx->kd->p[den].fDensity;
+ }
+ int nb = 0;
+ for (j=0, hp=smx->hash;j<smx->nHashLength; j++,hp++)
+ if (hp->nGroup1>=0)nb++;
+ my_comm->ngroups = smx->nGroups;
+ my_comm->nb = nb;
+ my_comm->g1vec = vector(0,nb);
+ my_comm->g2vec = vector(0,nb);
+ my_comm->fdensity = vector(0,smx->nHashLength);
+ nb = 0;
+ for (j=0, hp=smx->hash;j<smx->nHashLength; j++,hp++)
+ if (hp->nGroup1>=0){
+ my_comm->g1vec[nb] = hp->nGroup1;
+ my_comm->g2vec[nb] = hp->nGroup2;
+ my_comm->fdensity[nb++] = hp->fDensity;
+ }
+ return;
+}
+
+
+/* ================================================================== */
+/* ======================= Sorting ================================== */
+/* ================================================================== */
+
+typedef struct index_struct {
+ float value;
+ int index;
+} *ptrindex;
+
+int cmp_index(const void *a, const void *b)
+{
+ if ( ((ptrindex)a)->value<((ptrindex)b)->value) return -1;
+ else if ( ((ptrindex)a)->value>((ptrindex)b)->value) return 1;
+ else return 0;
+}
+
+void make_rank_table(int n, int *ivect, int *rank)
+/* Given a vector of integers ivect[1..n], construct a rank table rank[1..n]
+so that rank[j] contains the ordering of element j, with rank[j]=n indicating
+that the jth element was the highest, and rank[j]=1 indicating that it
+was the lowest. Storage for rank[] should be declared externally */
+/* I don't think this routine is particularly fast, but it's a
+miniscule fraction of the runtime */
+{
+ int j;
+ ptrindex sortvect;
+
+ sortvect = (ptrindex)malloc(n*sizeof(struct index_struct));
+ for (j=0;j<n;j++) sortvect[j].value = (float)ivect[j+1];
+ for (j=0;j<n;j++) sortvect[j].index = j+1; /* Label them prior to sort */
+ qsort(sortvect,n,sizeof(struct index_struct),cmp_index);
+ /* Now sortvect is in order (smallest to largest) */
+ for (j=0;j<n;j++) rank[sortvect[j].index]=j+1;
+ free(sortvect);
+ return;
+}
+
+/* DJE -- This is a C-translation of the Slatec FORTRAN routine ssort(). */
+/* I have kept the variable names and program flow unchanged; hence
+all the goto's and capitalized names.... */
+/* The input vectors must be UNIT-OFFSET */
+
+/* Here is the original Slatec header:
+
+***BEGIN PROLOGUE SSORT
+***PURPOSE Sort an array and optionally make the same interchanges in
+ an auxiliary array. The array may be sorted in increasing
+ or decreasing order. A slightly modified QUICKSORT
+ algorithm is used.
+***LIBRARY SLATEC
+***CATEGORY N6A2B
+***TYPE SINGLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
+***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING
+***AUTHOR Jones, R. E., (SNLA)
+ Wisniewski, J. A., (SNLA)
+***DESCRIPTION
+
+ SSORT sorts array X and optionally makes the same interchanges in
+ array Y. The array X may be sorted in increasing order or
+ decreasing order. A slightly modified quicksort algorithm is used.
+
+ Description of Parameters
+ X - array of values to be sorted (usually abscissas)
+ Y - array to be (optionally) carried along
+ N - number of values in array X to be sorted
+ KFLAG - control parameter
+ = 2 means sort X in increasing order and carry Y along.
+ = 1 means sort X in increasing order (ignoring Y)
+ = -1 means sort X in decreasing order (ignoring Y)
+ = -2 means sort X in decreasing order and carry Y along.
+
+***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm
+ for sorting with minimal storage, Communications of
+ the ACM, 12, 3 (1969), pp. 185-187.
+***ROUTINES CALLED XERMSG
+***REVISION HISTORY (YYMMDD)
+ 761101 DATE WRITTEN
+ 761118 Modified to use the Singleton quicksort algorithm. (JAW)
+ 890531 Changed all specific intrinsics to generic. (WRB)
+ 890831 Modified array declarations. (WRB)
+ 891009 Removed unreferenced statement labels. (WRB)
+ 891024 Changed category. (WRB)
+ 891024 REVISION DATE from Version 3.2
+ 891214 Prologue converted to Version 4.0 format. (BAB)
+ 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
+ 901012 Declared all variables; changed X,Y to SX,SY. (M. McClain)
+ 920501 Reformatted the REFERENCES section. (DWL, WRB)
+ 920519 Clarified error messages. (DWL)
+ 920801 Declarations section rebuilt and code restructured to use
+ IF-THEN-ELSE-ENDIF. (RWC, WRB)
+***END PROLOGUE SSORT
+*/
+
+#define TYPEOFY int /* DJE--To make the variable type of Y customizable */
+ /* because it has to be changed in two places...*/
+
+void ssort(float X[], TYPEOFY Y[], int N, int KFLAG)
+/* Note that the second array is an int array. If you want otherwise,
+alter the type above */
+{
+ /* .. Local Scalars .. */
+ float R, T, TT;
+ TYPEOFY TTY, TY;
+ int I, IJ, J, K, KK, L, M, NN;
+ /* .. Local Arrays .. */
+ int IL[31], IU[31]; /* DJE--These were 21, but I suspect that this
+ limits N to 2^21. Since memory is cheap, I'll set this a
+ little higher */
+
+ /* ***FIRST EXECUTABLE STATEMENT SSORT */
+ NN = N;
+ if (NN < 1) {
+ fprintf(stderr,"The number of values to be sorted is not positive.");
+ abort();
+ }
+
+ KK = abs(KFLAG);
+ if (KK != 1 && KK != 2) {
+ fprintf(stderr,"The sort control parameter, K, is not 2, 1, -1, or -2.");
+ abort();
+ }
+
+ /* Alter array X to get decreasing order if needed */
+
+ if (KFLAG <= -1)
+ for (I=1; I<=NN; I++)
+ X[I] = -X[I];
+
+
+ if (KK == 2) goto line100;
+
+ /* Sort X only */
+
+ M = 1;
+ I = 1;
+ J = NN;
+ R = 0.375E0;
+
+line20: if (I == J) goto line60;
+ if (R <= 0.5898437E0)
+ R = R+3.90625E-2;
+ else R = R-0.21875E0;
+
+
+line30: K = I;
+
+ /* Select a central element of the array and save it in location T */
+
+ IJ = I + (int)((J-I)*R);
+ T = X[IJ];
+
+ /* If first element of array is greater than T, interchange with T */
+
+ if (X[I] > T) {
+ X[IJ] = X[I];
+ X[I] = T;
+ T = X[IJ];
+ }
+ L = J;
+
+ /* If last element of array is less than than T, interchange with T */
+
+ if (X[J] < T) {
+ X[IJ] = X[J];
+ X[J] = T;
+ T = X[IJ];
+
+ /* If first element of array is greater than T, interchange with T */
+
+ if (X[I] > T) {
+ X[IJ] = X[I];
+ X[I] = T;
+ T = X[IJ];
+ }
+ }
+
+ /* Find an element in the second half of the array which is smaller */
+ /* than T */
+
+line40: L = L-1;
+ if (X[L] > T) goto line40;
+
+ /* Find an element in the first half of the array which is greater */
+ /* than T */
+
+line50: K = K+1;
+ if (X[K] < T) goto line50;
+
+ /* Interchange these elements */
+
+ if (K <= L) {
+ TT = X[L];
+ X[L] = X[K];
+ X[K] = TT;
+ goto line40;
+ }
+
+ /* Save upper and lower subscripts of the array yet to be sorted */
+
+ if (L-I > J-K) {
+ IL[M] = I;
+ IU[M] = L;
+ I = K;
+ M = M+1;
+ } else {
+ IL[M] = K;
+ IU[M] = J;
+ J = L;
+ M = M+1;
+ }
+ goto line70;
+
+ /* Begin again on another portion of the unsorted array */
+
+line60: M = M-1;
+ if (M == 0) goto line190;
+ I = IL[M];
+ J = IU[M];
+
+line70: if (J-I >= 1) goto line30;
+ if (I == 1) goto line20;
+ I = I-1;
+
+line80: I = I+1;
+ if (I == J) goto line60;
+ T = X[I+1];
+ if (X[I] <= T) goto line80;
+ K = I;
+
+line90: X[K+1] = X[K];
+ K = K-1;
+ if (T < X[K]) goto line90;
+ X[K+1] = T;
+ goto line80;
+
+ /* Sort X and carry Y along */
+
+line100: M = 1;
+ I = 1;
+ J = NN;
+ R = 0.375E0;
+
+line110: if (I == J) goto line150;
+ if (R <= 0.5898437E0)
+ R = R+3.90625E-2;
+ else R = R-0.21875E0;
+
+line120: K = I;
+
+ /* Select a central element of the array and save it in location T */
+
+ IJ = I + (int)((J-I)*R);
+ T = X[IJ];
+ TY = Y[IJ];
+
+ /* If first element of array is greater than T, interchange with T */
+
+ if (X[I] > T) {
+ X[IJ] = X[I];
+ X[I] = T;
+ T = X[IJ];
+ Y[IJ] = Y[I];
+ Y[I] = TY;
+ TY = Y[IJ];
+ }
+ L = J;
+;
+ /* If last element of array is less than T, interchange with T */
+
+ if (X[J] < T) {
+ X[IJ] = X[J];
+ X[J] = T;
+ T = X[IJ];
+ Y[IJ] = Y[J];
+ Y[J] = TY;
+ TY = Y[IJ];
+
+ /* If first element of array is greater than T, interchange with T */
+
+ if (X[I] > T) {
+ X[IJ] = X[I];
+ X[I] = T;
+ T = X[IJ];
+ Y[IJ] = Y[I];
+ Y[I] = TY;
+ TY = Y[IJ];
+ }
+ }
+
+ /* Find an element in the second half of the array which is smaller */
+ /* than T */
+
+line130: L = L-1;
+ if (X[L] > T) goto line130;
+
+ /* Find an element in the first half of the array which is greater */
+ /* than T */
+
+line140: K = K+1;
+ if (X[K] < T) goto line140;
+
+ /* Interchange these elements */
+
+ if (K <= L) {
+ TT = X[L];
+ X[L] = X[K];
+ X[K] = TT;
+ TTY = Y[L];
+ Y[L] = Y[K];
+ Y[K] = TTY;
+ goto line130;
+ }
+
+ /* Save upper and lower subscripts of the array yet to be sorted */
+
+ if (L-I > J-K) {
+ IL[M] = I;
+ IU[M] = L;
+ I = K;
+ M = M+1;
+ } else {
+ IL[M] = K;
+ IU[M] = J;
+ J = L;
+ M = M+1;
+ }
+ goto line160;
+
+ /* Begin again on another portion of the unsorted array */
+
+line150: M = M-1;
+ if (M == 0) goto line190;
+ I = IL[M];
+ J = IU[M];
+
+line160: if (J-I >= 1) goto line120;
+ if (I == 1) goto line110;
+ I = I-1;
+
+line170: I = I+1;
+ if (I == J) goto line150;
+ T = X[I+1];
+ TY = Y[I+1];
+ if (X[I] <= T) goto line170;
+ K = I;
+
+line180: X[K+1] = X[K];
+ Y[K+1] = Y[K];
+ K = K-1;
+ if (T < X[K]) goto line180;
+ X[K+1] = T;
+ Y[K+1] = TY;
+ goto line170;
+
+ /* Clean up */
+
+line190: if (KFLAG <= -1)
+ for (I=1; I<=NN; I++)
+ X[I] = -X[I];
+
+ return;
+}
Added: trunk/yt/lagos/hop/hop_kd.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop_kd.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,223 @@
+/* KD.C */
+/* This was written by Joachim Stadel and the NASA HPCC ESS at
+the University of Washington Department of Astronomy as part of
+the SMOOTH program, v2.0.1.
+URL: http://www-hpcc.astro.washington.edu/tools/SMOOTH */
+
+/* DJE--I have removed all the subroutines not used by HOP, notably
+the input and output routines. */
+
+/* HOP Version 1.0 (12/15/97) -- Original Release */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <sys/time.h>
+#include <sys/resource.h>
+#include <assert.h>
+#include "kd.h"
+//#include "macros_and_parameters.h"
+/* #include "tipsydefs.h" */ /* Don't need this, since I removed kdReadTipsy()*/
+
+
+#define MAX_ROOT_ITTR 32
+
+
+void kdTime(KD kd,int *puSecond,int *puMicro)
+{
+ struct rusage ru;
+
+ getrusage(0,&ru);
+ *puMicro = ru.ru_utime.tv_usec - kd->uMicro;
+ *puSecond = ru.ru_utime.tv_sec - kd->uSecond;
+ if (*puMicro < 0) {
+ *puMicro += 1000000;
+ *puSecond -= 1;
+ }
+ kd->uSecond = ru.ru_utime.tv_sec;
+ kd->uMicro = ru.ru_utime.tv_usec;
+ }
+
+
+int kdInit(KD *pkd,int nBucket)
+{
+ KD kd;
+
+ kd = (KD)malloc(sizeof(struct kdContext));
+ assert(kd != NULL);
+ kd->nBucket = nBucket;
+ *pkd = kd;
+ return(1);
+ }
+
+/*
+ ** JST's Median Algorithm
+ */
+int kdMedianJst(KD kd,int d,int l,int u)
+{
+ float fm;
+ int i,k,m;
+ PARTICLE *p,t;
+
+ p = kd->p;
+ k = (l+u)/2;
+ m = k;
+ while (l < u) {
+ m = (l+u)/2;
+ fm = p[m].r[d];
+ t = p[m];
+ p[m] = p[u];
+ p[u] = t;
+ i = u-1;
+ m = l;
+ while (p[m].r[d] < fm) ++m;
+ while (m < i) {
+ while (p[i].r[d] >= fm) if (--i == m) break;
+ /*
+ ** Swap
+ */
+ t = p[m];
+ p[m] = p[i];
+ p[i] = t;
+ --i;
+ while (p[m].r[d] < fm) ++m;
+ }
+ t = p[m];
+ p[m] = p[u];
+ p[u] = t;
+ if (k <= m) u = m-1;
+ if (k >= m) l = m+1;
+ }
+ return(m);
+ }
+
+
+void kdCombine(KDN *p1,KDN *p2,KDN *pOut)
+{
+ int j;
+
+ /*
+ ** Combine the bounds.
+ */
+ for (j=0;j<3;++j) {
+ if (p2->bnd.fMin[j] < p1->bnd.fMin[j])
+ pOut->bnd.fMin[j] = p2->bnd.fMin[j];
+ else
+ pOut->bnd.fMin[j] = p1->bnd.fMin[j];
+ if (p2->bnd.fMax[j] > p1->bnd.fMax[j])
+ pOut->bnd.fMax[j] = p2->bnd.fMax[j];
+ else
+ pOut->bnd.fMax[j] = p1->bnd.fMax[j];
+ }
+ }
+
+
+void kdUpPass(KD kd,int iCell)
+{
+ KDN *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]);
+ }
+ else {
+ l = c[iCell].pLower;
+ u = c[iCell].pUpper;
+ for (j=0;j<3;++j) {
+ c[iCell].bnd.fMin[j] = kd->p[u].r[j];
+ c[iCell].bnd.fMax[j] = kd->p[u].r[j];
+ }
+ for (pj=l;pj<u;++pj) {
+ for (j=0;j<3;++j) {
+ if (kd->p[pj].r[j] < c[iCell].bnd.fMin[j])
+ c[iCell].bnd.fMin[j] = kd->p[pj].r[j];
+ if (kd->p[pj].r[j] > c[iCell].bnd.fMax[j])
+ c[iCell].bnd.fMax[j] = kd->p[pj].r[j];
+ }
+ }
+ }
+ }
+
+int kdBuildTree(KD kd)
+{
+ int l,n,i,d,m,j,ct;
+ KDN *c;
+
+ n = kd->nActive;
+ kd->nLevels = 1;
+ l = 1;
+ while (n > kd->nBucket) {
+ n = n>>1;
+ l = l<<1;
+ ++kd->nLevels;
+ }
+ kd->nSplit = l;
+ kd->nNodes = l<<1;
+ kd->kdNodes = (KDN *)malloc(kd->nNodes*sizeof(KDN));
+ assert(kd->kdNodes != NULL);
+ /*
+ ** Set up ROOT node
+ */
+ c = kd->kdNodes;
+ c[ROOT].pLower = 0;
+ c[ROOT].pUpper = kd->nActive-1;
+ c[ROOT].bnd = kd->bnd;
+ i = ROOT;
+ ct = ROOT;
+ SETNEXT(ct);
+ while (1) {
+ if (i < kd->nSplit) {
+ d = 0;
+ for (j=1;j<3;++j) {
+ if (c[i].bnd.fMax[j]-c[i].bnd.fMin[j] >
+ c[i].bnd.fMax[d]-c[i].bnd.fMin[d]) d = j;
+ }
+ c[i].iDim = d;
+ m = kdMedianJst(kd,d,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-1;
+ c[UPPER(i)].bnd = c[i].bnd;
+ c[UPPER(i)].bnd.fMin[d] = c[i].fSplit;
+ c[UPPER(i)].pLower = m;
+ c[UPPER(i)].pUpper = c[i].pUpper;
+ i = LOWER(i);
+ }
+ else {
+ c[i].iDim = -1;
+ SETNEXT(i);
+ if (i == ct) break;
+ }
+ }
+ kdUpPass(kd,ROOT);
+ return(1);
+ }
+
+
+int cmpParticles(const void *v1,const void *v2)
+{
+ PARTICLE *p1=(PARTICLE *)v1,*p2=(PARTICLE *)v2;
+
+ return(p1->iOrder - p2->iOrder);
+ }
+
+
+void kdOrder(KD kd)
+{
+ qsort(kd->p,kd->nActive,sizeof(PARTICLE),cmpParticles);
+ }
+
+void kdFinish(KD kd)
+{
+ free(kd->p);
+ free(kd->kdNodes);
+ free(kd);
+ }
+
Added: trunk/yt/lagos/hop/hop_regroup.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop_regroup.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,714 @@
+/* REGROUP.C, Daniel Eisenstein, 1997 */
+/* Based on a paper by Daniel Eisenstein & Piet Hut,
+"HOP: A New Group-Finding Algorithm for N-body Simulations."
+See the included documentation or view it at
+http://www.sns.ias.edu/~eisenste/hop/hop_doc.html */
+
+/* Version 1.0 (12/15/97) -- Original Release */
+
+#include "slice.h"
+#include <string.h>
+#include <stdio.h>
+#include <limits.h>
+#include <float.h>
+//#include "macros_and_parameters.h"
+#include "hop.h"
+
+#define ISYM "d"
+#define GSYM "g"
+#define FSYM "f"
+
+/* #define MINDENS (-FLT_MAX/3.0) */
+#define MINDENS (-1.e+30/3.0)
+/* This is the most negative density that can be accomodated. Note
+that MINDENS*2.0 is referenced in the code and so must be properly
+represented by the machine. There's no reason for this to be close to
+the actual minimum of the density. */
+
+#define INFORM(pstr) printf(pstr); fflush(stdout)
+/* Used for messages, e.g. INFORM("Doing this"); */
+
+#define UNBOUND -2 /* The tag marker for unbound particles */
+
+/* ----------------------------------------------------------------------- */
+/* Prototypes */
+void initgrouplist(Grouplist *g);
+void readtags(Slice *s, Grouplist *g, char *fname);
+void densitycut(Slice *s, char *fname, float densthresh);
+void writegmerge(Slice *s, Grouplist *gl, char *fname, float pt, float mt);
+void readgmerge(Slice *s, Grouplist *gl, char *fname);
+void merge_groups_boundaries(Slice *s, Grouplist *gl, char *fname,
+ float peakdensthresh, float saddledensthresh, float densthresh, HC *my_comm);
+void translatetags(Slice *s, Grouplist *gl);
+void writetags(Slice *s, Grouplist *gl, char *fname);
+void writetagsf77(Slice *s, Grouplist *gl, char *fname);
+void count_membership(Slice *s, Grouplist *g);
+void sort_groups(Slice *s, Grouplist *gl, int mingroupsize, char *fname);
+
+/* ----------------------------------------------------------------------- */
+/* We use the following structure to handle the user interface: */
+
+typedef struct controlstruct {
+ char *tagname; /* Input file for group tags */
+ char *densname; /* Input file for density file */
+ char *gmergename; /* Input file for group boundary specifications, OR
+ input file for group merging data */
+ char *outsizename; /* Output file for size output*/
+ char *outtagname; /* Output file for group tags */
+ char *outgmergename; /* Output file for group merging */
+
+ int qdenscut; /* =1 if we're making a density cut, =0 otherwise */
+ float densthresh; /* The outer density threshold (delta_outer)*/
+
+ int qgbound; /* =1 if we are to read the boundaries file and
+ determine the merging.*/
+ float peak_thresh; /* Density threshold for peak (delta_peak) */
+ float saddle_thresh; /* Density threshold for merging (delta_saddle) */
+ int qgmerge_given; /* =1 if we are to use a group translation from file */
+
+ int mingroupsize; /* The minimum group size we follow */
+ int qoutput; /* =1 if we are to write the tags */
+ int qf77; /* =1 if binary output if in f77 format */
+ int qpipe; /* =1 if we are to write the output tags to stdout */
+ int qsort; /* =1 if we are to sort */
+
+ /* The following aren't used in the present version, but I included
+ them in case the user wants to customize the program: */
+ char *dataname; /* Input file for particle data */
+ int qunbind; /* =1 if we are to unbind at all */
+} Controls; /* Type Controls is defined */
+
+/* ====================================================================== */
+/* ===================== User Interface ================================= */
+/* ====================================================================== */
+
+void parsecommandline(float dens_outer, Controls *c)
+{
+ int narg, qmerge;
+ char *outname, *rootname;
+ narg = 1;
+ rootname = c->dataname = c->densname = c->gmergename = c->tagname =
+ outname = c->outsizename = c->outtagname = c->outgmergename = NULL;
+ c->qdenscut = -1;
+ qmerge = 1;
+ c->qgmerge_given = 0;
+
+ c->qunbind = 0;
+ c->qoutput = 1;
+ c->qsort = 1;
+ c->qpipe = 0;
+ c->qf77 = 0;
+
+ c->mingroupsize = -1;
+ if (2.0*MINDENS>=MINDENS || MINDENS>=0)
+ myerror("MINDENS seems to be illegal.");
+ /* Need MINDENS<0 and 2*MINDENS to be machine-representable */
+ c->densthresh = 2.0*MINDENS;
+ c->saddle_thresh = 2.0*MINDENS;
+ c->peak_thresh = 2.0*MINDENS;
+
+ /* GLB: hard-code some parameters. */
+
+ c->peak_thresh = 3.0*dens_outer;
+ c->saddle_thresh = 2.5*dens_outer;
+ c->densthresh = dens_outer;
+ c->qdenscut = 1;
+ rootname = "output_hop";
+
+ /* Get the input files ready */
+ if (c->qdenscut==-1) {
+ /* Neither -douter nor -nodens was chosen. */
+ mywarn("Outer density threshold left unspecified. Skipping this cut.");
+ c->qdenscut = 0;
+ } else if (c->qdenscut==1) {
+ /* We have a chosen density. Need to figure out the density file. */
+ if (c->densname==NULL) {
+ if (rootname==NULL)
+ myerror("No density file name or root has been specified.");
+ c->densname = (char *)malloc(80);
+ strcpy(c->densname,rootname); strcat(c->densname, ".den");
+ }
+ } else c->densname = NULL; /* We have no reason to read it */
+
+ if (c->tagname==NULL) {
+ if (rootname==NULL)
+ myerror("No .hop file name or root has been specified.");
+ c->tagname = (char *)malloc(80);
+ strcpy(c->tagname,rootname); strcat(c->tagname, ".hop");
+ }
+
+ if (qmerge==1) {
+ if (c->qgmerge_given==0) {
+ /* We need to have a .gbound file */
+ c->qgbound = 1;
+ if (c->saddle_thresh<MINDENS || c->peak_thresh<MINDENS)
+ myerror("-dsaddle and -dpeak need to be specified.");
+ if (c->gmergename==NULL) {
+ if (rootname==NULL)
+ myerror("No .gbound file name or root has been specified.");
+ c->gmergename = (char *)malloc(80);
+ strcpy(c->gmergename,rootname);
+ strcat(c->gmergename, ".gbound");
+ }
+ } else c->qgbound = 0; /* We know c->mergename is ready to go */
+ } else c->gmergename = NULL; /* No reason to read it */
+
+ /* Get the output files ready */
+ /* If a default name wasn't given, we'll assume zregroup */
+ if (outname==NULL) {
+ outname = (char *)malloc(20);
+ strcpy(outname,"zregroup");
+ }
+ /* Primary tag output: */
+ if (c->qoutput) { /* Need to figure out where we're sending the output */
+ if (c->qpipe&&c->outtagname)
+ myerror("Conflicting instructions--gave specific output name and told to pipe.");
+ if (c->qpipe>0) mywarn("Writing tags to stdout.");
+ if (c->qpipe) c->outtagname = NULL; /* Our signal to send to stdout */
+ else if (c->outtagname==NULL) {
+ c->outtagname = (char *)malloc(80);
+ strcpy(c->outtagname, outname);
+ strcat(c->outtagname, ".tag");
+ } /* Otherwise the name was set by the user */
+ } else {
+ /* We're not outputing tags */
+ if (c->qpipe) myerror("Conflicting instructions--told to pipe and not to output.");
+ }
+
+ if (c->qsort) {
+ if (c->qpipe>=0) { /* The user didn't specify quiet */
+ c->outsizename = (char *)malloc(80);
+ strcpy(c->outsizename, outname);
+ strcat(c->outsizename, ".size");
+ }
+ }
+
+ if (c->qpipe>=0) { /* The user didn't specify quiet */
+ c->outgmergename = (char *)malloc(80);
+ strcpy(c->outgmergename, outname);
+ strcat(c->outgmergename, ".gmerge");
+ }
+
+ if (c->mingroupsize >= 0 && !c->qsort)
+ myerror("Imposition of a certain group size occurs within the sort routine.");
+ if (c->qsort && c->mingroupsize < 0) {
+ mywarn("No minimum group size specified. Assuming 10 particles.");
+ c->mingroupsize = 10;
+ }
+
+ if (c->densthresh<MINDENS) c->densthresh=MINDENS;
+ /* This is our default--a very negative number */
+
+ return;
+}
+
+/* ====================================================================== */
+/* ============================== MAIN() ================================ */
+/* ====================================================================== */
+
+/* void main(int argc, char *argv[]) */
+void regroup_main(float dens_outer, HC *my_comm)
+{
+ Grouplist *gl = my_comm->gl;
+ Slice *s = my_comm->s;
+ FILE *f;
+ Controls c;
+
+ /* parsecommandline(argc, argv, &c); */
+ parsecommandline(dens_outer, &c);
+
+ //initgrouplist(gl);
+ //s=newslice();
+
+ /* We need to read the tag file and perhaps perform a density cut */
+ // We don't read anymore (mjt)
+ //readtags(s,gl,c.tagname);
+
+ // We cut in advance now (mjt)
+ //if (c.qdenscut) densitycut(s,c.densname,c.densthresh);
+
+ /* Next do the merging of the input groups */
+ if (c.qgbound) {
+ /* We're going to read a .gbound file and merge groups */
+ merge_groups_boundaries(s,gl,c.gmergename,
+ c.peak_thresh, c.saddle_thresh, c.densthresh, my_comm);
+ /* Renumber the groups from large to small; remove any tiny ones */
+ if (c.qsort) sort_groups(s, gl, c.mingroupsize, c.outsizename);
+ writegmerge(s, gl, c.outgmergename, c.peak_thresh, c.saddle_thresh);
+ translatetags(s,gl);
+ }
+ else if (c.qgmerge_given) {
+ /* We're going to read a .gmerge file and merge groups as it says */
+ readgmerge(s, gl, c.gmergename);
+ translatetags(s, gl);
+ } /* Else we'll use the tags as given by the original .hop file */
+
+ /* If one wants to manipulate the groups any more, this is a good
+ place to do it. For example, you might want to remove unbound particles:
+ if (c.qunbind) {
+ get_particle_data(s, gl, c.dataname);
+ unbind_particles(s, gl, c.mingroupsize);
+ }
+ */
+
+ /* Write the output */
+ /*if (c.qoutput) {
+ if (c.qf77) writetagsf77(s, gl, c.outtagname);
+ else writetags(s, gl, c.outtagname);
+ }*/
+
+ //free_slice(s);
+ return;
+}
+
+/* ================================================================= */
+/* =================== Initialization Routines ===================== */
+/* ================================================================= */
+
+void initgrouplist(Grouplist *g)
+/* Just make sure this stuff is zero */
+{
+ g->list = NULL;
+ g->npartingroups = g->npart = g->ngroups = 0; g->nnewgroups = 0;
+ return;
+}
+
+void readtags(Slice *s, Grouplist *g, char *fname)
+/* Read the tag file named fname into s->ntag[] */
+/* Groups need not be sorted, but must be numbered from 0 to ngroups-1 */
+{
+ FILE *f;
+
+ if ((f=fopen(fname,"r"))==NULL) myerror("Input tag file not found.");
+ if (fread(&(g->npart),sizeof(int),1,f)!=1) myerror("Tag file read error.");
+ if (fread(&(g->ngroups),sizeof(int),1,f)!=1) myerror("Tag file read error.");
+ fprintf(stderr,"Number of particles: %"ISYM". Number of groups: %"ISYM".\n",
+ g->npart, g->ngroups);
+
+ s->numpart = g->npart;
+ s->numlist = g->npart;
+ s->ntag = ivector(1,s->numlist);
+ s->ID = ivector(1,s->numlist);
+ fread(s->ntag+1, sizeof(int), s->numlist, f); /* Read in all the tags */
+ fread(s->ID+1, sizeof(int), s->numlist,f); /* Read in the real particle IDs. S Skory */
+ fclose(f);
+
+ return;
+}
+
+/* ========================== Density Cut ======================== */
+
+#define MAXBLOCK 65536 /* Read the file 256k at a time */
+
+void densitycut(Slice *s, char *fname, float densthresh)
+/* Read the density file and change the tag on any particle with density
+less than densthresh to -1, thus removing them from groups */
+/* This will leave some groups with no particles, which is fine */
+/* We read the file in segments, so as to reduce memory consumption */
+{
+ FILE *f;
+ int j, numread, npart, block; /* block was a float by mistake */
+ float density[MAXBLOCK];
+
+ if ((f=fopen(fname,"r"))==NULL)
+ myerror("Density file not found.");
+ npart = 0; fread(&npart,sizeof(int),1,f);
+ if (npart!=s->numpart)
+ mywarn("Density file doesn't match slice description.");
+
+ numread = 0;
+ block = MAXBLOCK; /* Start off big */
+ while (numread<npart) {
+ if (npart-numread<block) block = npart-numread;
+ if (fread(density,sizeof(float),block,f)!=block)
+ myerror("Read error in density file.");
+ for (j=1;j<=block;j++)
+ if (density[j-1]<densthresh) /* density is zero-offset */
+ s->ntag[numread+j]=(-1); /* s->ntag is unit-offset */
+ numread+=block;
+ }
+ fclose(f);
+ return;
+}
+
+/* ====================== Read/Write .gmerge files ======================= */
+/* The gmerge file is just a map from the old (pre-regroup) group numbers
+to the new (post-regroup) group numbers. Of course, there are more "old"
+groups than "new" groups, since the point of regroup() is to merge groups. */
+
+void writegmerge(Slice *s, Grouplist *gl, char *fname, float pt, float mt)
+/* Write the translation between old groups and new groups, ASCII */
+{
+ FILE *f;
+ int j;
+ Group *gr;
+
+ if (fname==NULL) return; /* We've been told not to write anything */
+
+ if ((f=fopen(fname,"w"))==NULL) myerror("Can't open gmerge file for write.");
+ fprintf(f,"%"ISYM"\n%"ISYM"\n%"ISYM"\n", gl->npart, gl->ngroups, gl->nnewgroups);
+ fprintf(f,"%"FSYM"\n%"FSYM"\n", pt, mt);
+ for (j=0,gr=gl->list;j<gl->ngroups;j++,gr++)
+ fprintf(f,"%"ISYM" %"ISYM"\n", j, gr->idmerge);
+ fclose(f);
+ return;
+}
+
+void readgmerge(Slice *s, Grouplist *gl, char *fname)
+/* Read the translation between old groups and new groups, ASCII */
+/* Also, set up gl->list for translation */
+{
+ FILE *f;
+ int j, dummy;
+ Group *gr;
+ float pt, mt;
+
+ if ((f=fopen(fname,"r"))==NULL) myerror("Can't open gmerge read file.");
+ if (fscanf(f,"%"ISYM"\n%"ISYM"\n%"ISYM"\n", &(gl->npart), &(gl->ngroups),
+ &(gl->nnewgroups))!=3) myerror("Error in header of gmerge file.");
+ if (gl->npart!=s->numpart) myerror("Number of particles in gmerge file doesn't match that of tags file.");
+ fscanf(f,"%"FSYM" %"FSYM"\n", &pt, &mt);
+
+ if (gl->list!=NULL) free(gl->list);
+ gl->list = (Group *)malloc((size_t)(gl->ngroups *sizeof(Group)));
+ if (gl->list==NULL) myerror("Error in allocating gl->list.");
+
+ for (j=0,gr=gl->list; j<gl->ngroups; j++,gr++) {
+ if (fscanf(f,"%"ISYM" %"ISYM"\n", &dummy, &(gr->idmerge))!=2 || dummy!=j)
+ myerror("Error in reading gmerge file.");
+ gr->npart = -1; /* We're not setting this */
+ }
+ fclose(f);
+ return;
+}
+
+/* ====================== GROUP MERGING BY BOUNDARIES ================ */
+
+void merge_groups_boundaries(Slice *s, Grouplist *gl, char *mergename,
+ float peakdensthresh, float saddledensthresh, float densthresh,
+ HC *my_comm)
+/* Read in the gmerge file and decide which groups are to be merged.
+Groups are numbered 0 to ngroups-1. Groups with boundaries greater
+than saddledensthresh are merged. Groups with maximum densities
+less than peakdensthresh are merged to the group with
+maxdensity above peakdensthresh with which it shares the highest
+density border. */
+/* Only groups with maximum densities above peakdensthresh can be group
+centers. */
+/* Allocate space for the grouplist and store the merging results in
+the idmerge field. */
+/* I think this will work even if saddledensthresh<densthresh */
+{
+ int j, k, g1, g2, ngroups, dummy[3];
+ Group *gr;
+ float *densestbound, fdum[3], dens;
+ int *densestboundgroup, changes;
+ char line[80], *tempfilename;
+ FILE *fp;
+ FILE *boundfp;
+ float *gdensity = my_comm->gdensity;
+ ngroups = my_comm->ngroups;
+
+ tempfilename = tmpnam(NULL);
+ if (densthresh<MINDENS) densthresh=MINDENS;
+ /* Have a 2*MINDENS condition below... */
+ densestbound = vector(0,ngroups-1);
+ densestboundgroup = ivector(0,ngroups-1);
+
+ /* Now allocate the grouplist */
+ gl->ngroups = ngroups;
+ if (gl->list!=NULL) free(gl->list);
+ gl->list = (Group *)malloc((size_t)(gl->ngroups *sizeof(Group)));
+ fprintf(stderr,"ngroups = %d\n",ngroups);
+ if (gl->list==NULL) myerror("Error in allocating gl->list.");
+ for (j=0,gr=gl->list;j<gl->ngroups;j++,gr++) {
+ /* If group is too underdense, it cannot be a group center */
+ if (gdensity[j]<peakdensthresh) gr->idmerge=(-1);
+ else gr->idmerge = j;
+ gr->npart = -1; /* Not doing anything with this */
+ densestbound[j] = 2.0*MINDENS; /* Initialize */
+ densestboundgroup[j] = -1; /* Initialize */
+ }
+
+ /* Now step through the list of boundaries */
+ /* If a boundary is between two groups with max densities above
+ peakdensthresh and if the boundary is above saddledensthresh, then
+ merge the groups (keeping the lower number of the two). */
+ /* If one of the groups is above peakdensthresh and the other is
+ 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 */
+
+ if ((boundfp=fopen(tempfilename,"w"))==NULL)
+ myerror("Error opening scratch file");
+
+ 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);
+ continue; /* group isn't dense enough */
+ }
+ if (gdensity[g1]>=peakdensthresh && gdensity[g2]>=peakdensthresh)
+ if (dens<saddledensthresh) continue;
+ /* Boundary is not dense enough to merge */
+ else { /* Groups should be merged */
+ /* Trace each group to its root */
+ while (g1!=gl->list[g1].idmerge)
+ g1=gl->list[g1].idmerge;
+ while (g2!=gl->list[g2].idmerge)
+ g2=gl->list[g2].idmerge;
+ if (g1<g2) gl->list[g2].idmerge=g1;
+ else gl->list[g1].idmerge=g2;
+ continue; /* Go to the next boundary */
+ }
+ /* Else one is above peakdensthresh, the other below. */
+ /* Make the high one g1 */
+ if (gdensity[g1]<gdensity[g2]) {
+ dummy[0] = g2; g2=g1; g1=dummy[0];
+ }
+ if (dens>densestbound[g2]) {
+ /* It's the densest boundary yet */
+ densestbound[g2] = dens;
+ densestboundgroup[g2] = g1;
+ }
+ } /* Get the next boundary line */
+ fclose(boundfp);
+
+
+ /* Now the fringe groups are connected to the proper group
+ (>peakdensthresh) with the largest boundary. But we want to look
+ through the boundaries between fringe groups to propagate this
+ along. Connections are only as good as their smallest boundary */
+ /* 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.");
+ /* If the density of this boundary and the densestbound of
+ the other group is higher than a group's densestbound, then
+ replace it. */
+ /* Make g1 have the higher densestbound */
+ if (densestbound[g2]>densestbound[g1]) {
+ dummy[0] = g2; g2=g1; g1=dummy[0];
+ }
+ if (dens>densestbound[g2]&&densestbound[g1]>densestbound[g2]) {
+ changes++;
+ if (dens<densestbound[g1]) densestbound[g2]=dens;
+ else densestbound[g2]=densestbound[g1];
+ densestboundgroup[g2] = densestboundgroup[g1];
+ }
+ }
+ fclose(boundfp);
+ } while (changes);
+
+ /* Now connect the low-density groups to their densest boundaries */
+ /* But only if the boundary exceeds densthresh! */
+ for (j=0;j<gl->ngroups;j++) {
+ if (densestbound[j]>=densthresh)
+ gl->list[j].idmerge = densestboundgroup[j];
+ }
+ /* Now we want to number the newly merged groups */
+ /* The center groups are given negative numbers <-1 */
+ for (j=0,gl->nnewgroups=0; j<gl->ngroups; j++)
+ if (gl->list[j].idmerge==j) {
+ gl->list[j].idmerge = -2-(gl->nnewgroups++);
+ }
+
+ /* Now trace each group through until a negative number is reached */
+ for (j=0; j<gl->ngroups; j++) {
+ if (gl->list[j].idmerge<0) continue;
+ g1 = j;
+ while ((g1=gl->list[g1].idmerge)>=0);
+ g2 = j;
+ do gl->list[g2].idmerge = g1;
+ while ((g2=gl->list[g2].idmerge)>=0);
+ }
+
+ /* Finally, renumber the groups 0..N-1 */
+ for (j=0,gr=gl->list;j<gl->ngroups;j++,gr++)
+ gr->idmerge = -2-gr->idmerge; /* Keep -1 -> -1 */
+
+
+ /* And delete the tempfile */
+ remove(tempfilename);
+ free_vector(gdensity,0,ngroups-1);
+ free_vector(densestbound,0,ngroups-1);
+ free_ivector(densestboundgroup,0,ngroups-1);
+ return;
+}
+
+/* ======================================================================= */
+/* =============== Update the tags and write them out ==================== */
+/* ======================================================================= */
+
+void translatetags(Slice *s, Grouplist *gl)
+/* Alter s->ntag to have the new groups. Reset gl so as to reflect the
+new number of groups. */
+{
+ int j;
+
+
+ for (j=1;j<=s->numlist;j++)
+ if (s->ntag[j]>=0) {
+ s->ntag[j] = gl->list[s->ntag[j]].idmerge;
+ }
+ /* Otherwise, translate the unbound particles */
+ else if (s->ntag[j]<-1)
+ s->ntag[j] = UNBOUND - gl->list[UNBOUND-s->ntag[j]].idmerge;
+ free(gl->list);
+ gl->list = NULL;
+ gl->ngroups = gl->nnewgroups;
+ return;
+}
+
+void writetags(Slice *s, Grouplist *gl, char *fname)
+/* Write s->ntag to file */
+/* If fname==NULL, write to stdout */
+{
+ FILE *f;
+
+
+ if (fname!=NULL) {
+ if ((f=fopen(fname,"w"))==NULL) myerror("Error opening new tag file.");
+ } else f=stdout;
+ fwrite(&(s->numpart),sizeof(int),1,f);
+ printf("writetags: s->numpart = %d gl->ngroups = %d\n",
+ s->numpart, gl->ngroups);
+ fwrite(&(gl->ngroups),sizeof(int),1,f);
+ fwrite(s->ntag+1,sizeof(int),s->numlist,f);
+ fwrite(s->ID+1,sizeof(int),s->numlist,f); /* S Skory */
+ fclose(f);
+
+ return;
+}
+
+void writetagsf77(Slice *s, Grouplist *gl, char *fname)
+/* Write s->ntag to file */
+/* If fname==NULL, write to stdout */
+/* Use a format readable for FORTRAN unformatted read commands */
+{
+ FILE *f;
+ int dummy;
+ if (fname!=NULL) {
+ if ((f=fopen(fname,"w"))==NULL) myerror("Error opening new tag file.");
+ } else f=stdout;
+ dummy = 8; fwrite(&dummy,sizeof(int),1,f);
+ fwrite(&(s->numpart),sizeof(int),1,f);
+ fwrite(&(gl->ngroups),sizeof(int),1,f);
+ fwrite(&dummy,sizeof(int),1,f);
+ dummy = s->numlist*sizeof(int); fwrite(&dummy,sizeof(int),1,f);
+ fwrite(s->ntag+1,sizeof(int),s->numlist,f);
+ fwrite(&dummy,sizeof(int),1,f);
+ fclose(f);
+ return;
+}
+
+/* ====================================================================== */
+/* ========================== Sorting the Groups ======================== */
+/* ====================================================================== */
+
+void sort_groups(Slice *s, Grouplist *gl, int mingroupsize, char *fname)
+/* Sort the groups, as labeled by the idmerge field not their original
+number, from largest to smallest. Alter the idmerge field to this new
+numbering, setting any below mingroupsize to -1. */
+/* If fname!=NULL, write a little output file listing the group sizes */
+{
+ FILE *f;
+ int j,k, *order, partingroup, igr, *newnum, nmergedgroups;
+ float *gsize;
+ Group *c;
+ void make_index_table(int n, float *fvect, int *index);
+
+ nmergedgroups = gl->nnewgroups;
+ gsize = vector(0,nmergedgroups-1);
+ order = ivector(1,nmergedgroups);
+ newnum = ivector(0,nmergedgroups-1);
+
+ /* First we need to find the number of particles in each group */
+ for (j=0,c=gl->list;j<gl->ngroups;j++,c++) c->npart=0;
+
+ for (j=1;j<=s->numlist;j++) { /* Look through all the particles */
+ igr = s->ntag[j];
+ if (igr>=0)
+ if (igr<gl->ngroups) gl->list[igr].npart++;
+ else myerror("Group tag is out of bounds.");
+ }
+ /* Now combine these to find the number in the new groups */
+ for (j=0;j<nmergedgroups;j++) gsize[j]=0;
+ for (j=0,c=gl->list;j<gl->ngroups;j++,c++)
+ if (c->idmerge>=0 && c->idmerge<nmergedgroups)
+ gsize[c->idmerge]+=c->npart;
+ else if (c->idmerge>=nmergedgroups)
+ myerror("Group idmerge is out of bounds.");
+
+ make_index_table(nmergedgroups, gsize-1, order);
+ /* But remember that order[] thinks that gsize is unit-offset */
+ for (j=nmergedgroups,k=0;j>0; j--,k++)
+ if (gsize[order[j]-1]>mingroupsize-0.5) newnum[order[j]-1]=k;
+ else break; /* All of the rest are too small */
+
+ gl->nnewgroups = k;
+ for (;j>0;j--) newnum[order[j]-1]=(-1);
+ /* Newnum[] holds the new sorted number for merged group j */
+
+ /* Now assign sorted group numbers to idmerge */
+ partingroup = 0;
+ for (j=0,c=gl->list;j<gl->ngroups;j++,c++)
+ if (c->idmerge>=0)
+ if ((c->idmerge = newnum[c->idmerge])>=0)
+ partingroup+=c->npart;
+
+ /* Output the .size file, if inputed name isn't NULL */
+ if (fname!=NULL) {
+ f = fopen(fname,"w");
+ fprintf(f,"%"ISYM"\n%"ISYM"\n%"ISYM"\n", s->numpart, partingroup, gl->nnewgroups);
+ for (j=0;j<gl->nnewgroups;j++)
+ fprintf(f,"%"ISYM" %"ISYM"\n", j, (int)gsize[order[nmergedgroups-j]-1]);
+ }
+ fclose(f);
+ free_ivector(order,1,nmergedgroups);
+ free_vector(gsize,0,nmergedgroups-1);
+ free_ivector(newnum,0,nmergedgroups-1);
+ return;
+}
+
+/* ======================== Sorting ============================ */
+
+typedef struct index_struct {
+ float value;
+ int index;
+} *ptrindex;
+
+int cmp_index_regroup(const void *a, const void *b)
+{
+ if ( ((ptrindex)a)->value<((ptrindex)b)->value) return -1;
+ else if ( ((ptrindex)a)->value>((ptrindex)b)->value) return 1;
+ else return 0;
+}
+
+void make_index_table(int n, float *fvect, int *index)
+/* Given a vector of floats fvect[1..n], construct a index table index[1..n]
+so that index[j] contains the ID number of the jth lowest element.
+Storage for index[] should be declared externally */
+/* This isn't fast, but it takes a tiny fraction of the runtime */
+{
+ int j;
+ ptrindex sortvect;
+
+ sortvect = (ptrindex)malloc(n*sizeof(struct index_struct));
+ for (j=0;j<n;j++) sortvect[j].value = fvect[j+1];
+ for (j=0;j<n;j++) sortvect[j].index = j+1; /* Label them prior to sort */
+ qsort(sortvect,n,sizeof(struct index_struct),cmp_index_regroup);
+ /* Now sortvect is in order (smallest to largest) */
+ for (j=0;j<n;j++) index[j+1] = sortvect[j].index;
+ free(sortvect);
+ return;
+}
Added: trunk/yt/lagos/hop/hop_slice.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop_slice.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,393 @@
+/* SLICE.C, Daniel Eisenstein, 1997 */
+/* Based on a paper by Daniel Eisenstein & Piet Hut,
+"HOP: A New Group-Finding Algorithm for N-body Simulations."
+See the included documentation or view it at
+http://www.sns.ias.edu/~eisenste/hop/hop_doc.html */
+
+/* Version 1.0 (12/15/97) -- Original Release */
+
+#include "slice.h"
+//#include "macros_and_parameters.h"
+
+/* ===================== Creating and Destroying Slices ================ */
+
+Slice *newslice()
+/* Returns a pointer to a newly allocated Slice variable, with internal
+variables initialized */
+{
+ Slice *s;
+ s = (Slice *)malloc(sizeof(Slice));
+ s->pid = NULL; s->offset = 0;
+ s->px = s->py = s->pz = s->vx = s->vy = s->vz = NULL;
+ s->ntag = NULL;
+ s->ID = NULL; /* S Skory */
+ s->numpart = s->numlist = 0;
+ return s;
+}
+
+void free_tags(Slice *s)
+/* Free the tag vector */
+{
+ if (s->ntag!=NULL) {
+ free_ivector(s->ntag, 1, s->numlist);
+ s->ntag=NULL;
+ free_ivector(s->ID, 1, s->numlist); /* S Skory */
+ s->ID=NULL;
+ }
+ return;
+}
+
+void free_data(Slice *s)
+/* Free all the data vectors */
+{
+ if (s->pid!=NULL) {free(s->pid); s->pid=NULL;}
+ if (s->px!=NULL) {free_vector(s->px,1,s->numlist); s->px=NULL;}
+ if (s->py!=NULL) {free_vector(s->py,1,s->numlist); s->py=NULL;}
+ if (s->pz!=NULL) {free_vector(s->pz,1,s->numlist); s->pz=NULL;}
+ if (s->vx!=NULL) {free_vector(s->vx,1,s->numlist); s->vx=NULL;}
+ if (s->vy!=NULL) {free_vector(s->vy,1,s->numlist); s->vy=NULL;}
+ if (s->vz!=NULL) {free_vector(s->vz,1,s->numlist); s->vz=NULL;}
+ return;
+}
+
+void free_slice(Slice *s)
+/* Free the space associated with the vectors in the given Slice */
+/* Then free the Slice variable itself */
+{
+ free_tags(s);
+ free_data(s);
+ free(s);
+ return;
+}
+
+/* =================================================================== */
+
+int f77write(FILE *f, void *p, int len)
+/* len is number of bytes to be written from p[0..len-1] */
+/* Return 0 if successful, 1 if not */
+{
+ if (fwrite(&len,sizeof(int),1,f)!=1) return 1;
+ if (fwrite(p,1,len,f)!=len) return 1;
+ if (fwrite(&len,sizeof(int),1,f)!=1) return 1;
+ return 0;
+}
+
+int f77read(FILE *f, void *p, int maxbytes)
+/* Read a FORTRAN style block from the given file */
+/* maxbytes is the amount of space the pointer p points to */
+/* Space must be allocated to read the whole block into p */
+/* Return amount read, scream if there's a problem */
+/* Reading is done ZERO-OFFSET */
+{
+ int size, size2;
+ if (fread(&size,sizeof(int),1,f)!=1)
+ myerror("f77read(): Error reading begin delimiter.");
+ if (size>maxbytes)
+ myerror("f77read(): Block delimiter exceeds size of storage.");
+ if (size<maxbytes)
+ mywarn("f77read(): Block size is smaller than size of storage.");
+ if (fread(p,1,size,f)!=size) myerror("f77read(): Error reading data.");
+ if (fread(&size2,sizeof(int),1,f)!=1)
+ myerror("f77read(): Error reading end delimiter.");
+ if (size!=size2) myerror("f77read(): Delimiters do not match.");
+ return size;
+}
+
+/* =================================================================== */
+/* The following are public-domain routines from Numerical Repices in C,
+2nd edition, by Press, Teulkolsky, Vetterling, & Flannery, 1992, Cambridge
+Univ. Press */
+
+#define NR_END 1
+#define FREE_ARG char*
+
+float *vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+ float *v;
+
+ v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+ if (!v) myerror("allocation failure in vector()");
+ return v-nl+NR_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+ int *v;
+
+ v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+ if (!v) myerror("allocation failure in ivector()");
+ return v-nl+NR_END;
+}
+void free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+#ifdef NOT_USED
+/* =================================================================== */
+/* =================================================================== */
+
+/* Everything from here on down may be expendable */
+
+int read_header(FILE *f, Slice *s)
+/* Read the header from file f. Place information and calculated
+numbers into the Slice variable, unless *s is NULL, in which case just
+skip the header */
+/* Return 0 if successful, non-zero otherwise */
+{
+ int sizeheader[2];
+ float header[100];
+
+ f77read(f, header, 400);
+ f77read(f, sizeheader, 8);
+
+ if (s==NULL) return 0; /* We've been instructed just to skip the
+ header without recording the information */
+ /* Now read the interesting information out of these arrays */
+ s->numpart = sizeheader[0];
+ s->numblocks = sizeheader[1];
+ s->numperblock = sizeheader[0]/sizeheader[1];
+ if (s->numpart != s->numblocks*(s->numperblock))
+ myerror("Number of blocks not an even divisor of number of particles.");
+
+ s->z = header[0];
+ s->boxsize = header[1]*1000.0; /* We use kpc, not Mpc */
+ s->physsize = s->boxsize/(1.0+s->z); /* We use kpc, not Mpc */
+ s->velscale = 100.0*header[1]*sqrt(3.0/8.0/PI)/(1.0+s->z);
+ /* To go from data to pec vel */
+ s->omega = header[4];
+ if (header[6]!=0.0) myerror("HDM component listed in header.");
+ s->lambda = header[7];
+ s->h0 = header[8];
+ s->sigma8 = header[9]; /* At z=0 */
+
+ /* Now find some computed quantities. */
+ s->a = 1.0/(1.0+s->z);
+ s->curv = 1.0-s->omega-s->lambda;
+ s->gamma = s->omega*(s->h0);
+ s->specn = 1.0;
+ s->hubb = 0.1*sqrt(s->omega/CUBE(s->a)+s->curv/SQR(s->a)+s->lambda)*(s->a);
+
+ /* The following assume Omega = 1 */
+ s->masspart = RHOCRIT/s->numpart*CUBE(s->boxsize);
+ s->growth = s->a;
+ s->t = HTIME*(s->h0)*pow(s->a, 1.5);
+ return 0;
+}
+
+void normalizedata(Slice *s, int conp, int conv)
+/* Put raw data into comoving h^-1 kpc and km/s units */
+{
+ int j;
+ float velnorm;
+ if (conp) {
+ for (j=1;j<=s->numlist;j++) s->px[j] *= s->boxsize;
+ for (j=1;j<=s->numlist;j++) s->py[j] *= s->boxsize;
+ for (j=1;j<=s->numlist;j++) s->pz[j] *= s->boxsize;
+ }
+
+ if (conv) {
+ for (j=1;j<=s->numlist;j++) s->vx[j] *= s->velscale;
+ for (j=1;j<=s->numlist;j++) s->vy[j] *= s->velscale;
+ for (j=1;j<=s->numlist;j++) s->vz[j] *= s->velscale;
+ }
+ return;
+}
+
+/* ================================================================ */
+
+int read_alldata(FILE *f, FILE *ftag, Slice *s, int conp, int conv)
+/* Read all the data, including the tags if ftag!=NULL. */
+/* Store positions and velocities unless conp or conv = 0 */
+/* Assume that the data file header has been skipped, but read the
+tag file header. */
+{
+ int block;
+ float *dummylist;
+
+ dummylist = NULL;
+ if (!conp || !conv) dummylist = vector(1,s->numperblock);
+ if (s->pid!=NULL)
+ mywarn("Non-NULL s->pid[] passed to read_alldata(). Ignoring...");
+
+ s->numlist=s->numpart;
+ if (conp) {
+ s->px=vector(1,s->numlist);
+ s->py=vector(1,s->numlist);
+ s->pz=vector(1,s->numlist);
+ }
+ if (conv) {
+ s->vx=vector(1,s->numlist);
+ s->vy=vector(1,s->numlist);
+ s->vz=vector(1,s->numlist);
+ }
+ if (ftag!=NULL) {
+ s->ntag = ivector(1,s->numlist);
+ s->ID = ivector(1,s->numlist); /* S Skory */
+ }
+
+
+ printf("Reading data...");
+ for (block=0;block<s->numblocks;block++) {
+ /* Read the three position blocks */
+ if (conp) { /* Store the data */
+ f77read(f, s->px+s->numperblock*block+1, s->numperblock*sizeof(float));
+ f77read(f, s->py+s->numperblock*block+1, s->numperblock*sizeof(float));
+ f77read(f, s->pz+s->numperblock*block+1, s->numperblock*sizeof(float));
+ } else { /* Don't store the data */
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ }
+ /* Now read the three velocity blocks */
+ if (conv) { /* Store the data */
+ f77read(f, s->vx+s->numperblock*block+1, s->numperblock*sizeof(float));
+ f77read(f, s->vy+s->numperblock*block+1, s->numperblock*sizeof(float));
+ f77read(f, s->vz+s->numperblock*block+1, s->numperblock*sizeof(float));
+ } else { /* Don't store the data */
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ f77read(f, dummylist+1, s->numperblock*sizeof(float));
+ }
+ if (block%8==1) {printf("."); fflush(stdout);}
+ }
+ if (dummylist!=NULL) free_vector(dummylist, 1, s->numperblock);
+ normalizedata(s,conp,conv);
+
+ if (ftag!=NULL) {
+ printf("tags..."); fflush(stdout);
+ readalltags(ftag, s);
+ }
+
+ printf("done!"); fflush(stdout);
+ return 0;
+}
+
+int read_partdata(FILE *f, FILE *ftag, Slice *s)
+/* Read one block (128k particles) of data into Slice s. Allocate needed
+storage, erasing and freeing previous storage. */
+/* This cannot be done with s->pid!=NULL, so s->pid is ignored and s->numlist
+is reset to BLOCKSIZE */
+/* Unlike other routines, this stores both positions and velocities in
+all cases (since the storage requirements are already small */
+/* If ftag==NULL, don't read the tag file. Otherwise do read it. */
+{
+ if (s->pid!=NULL)
+ mywarn("Non-trivial pid[] not supported with incremental reads");
+ /* If we need to reallocate memory, do it. Otherwise, just write over */
+ if (s->px==NULL || s->vx==NULL || s->numlist!=s->numperblock) {
+ if (s->px!=NULL) free_vector(s->px,1,s->numlist);
+ if (s->py!=NULL) free_vector(s->py,1,s->numlist);
+ if (s->pz!=NULL) free_vector(s->pz,1,s->numlist);
+ if (s->vx!=NULL) free_vector(s->vx,1,s->numlist);
+ if (s->vy!=NULL) free_vector(s->vy,1,s->numlist);
+ if (s->vz!=NULL) free_vector(s->vz,1,s->numlist);
+ if (ftag!=NULL && s->ntag!=NULL) {
+ free_ivector(s->ntag, 1, s->numlist);
+ free_ivector(s->ID, 1, s->numlist); /* S Skory */
+ }
+ s->numlist = s->numperblock;
+ s->px = vector(1,s->numlist);
+ s->py = vector(1,s->numlist);
+ s->pz = vector(1,s->numlist);
+ s->vx = vector(1,s->numlist);
+ s->vy = vector(1,s->numlist);
+ s->vz = vector(1,s->numlist);
+ if (ftag!=NULL) {
+ s->ntag = ivector(1, s->numlist);
+ s->ID = ivector(1, s->numlist); /* S Skory */
+ }
+ s->offset=0;
+ /* fprintf(stderr, "Reallocating data arrays.\n"); */
+ }
+ else s->offset+=s->numlist;
+
+ f77read(f,s->px+1,sizeof(float)*s->numlist);
+ f77read(f,s->py+1,sizeof(float)*s->numlist);
+ f77read(f,s->pz+1,sizeof(float)*s->numlist);
+ f77read(f,s->vx+1,sizeof(float)*s->numlist);
+ f77read(f,s->vy+1,sizeof(float)*s->numlist);
+ f77read(f,s->vz+1,sizeof(float)*s->numlist);
+
+ if (ftag!=NULL) readtag(ftag, s->numlist, s->ntag);
+ normalizedata(s,1,1);
+ return 0;
+}
+
+/* =============================================================== */
+
+int readtag(FILE *f, int numread, int *ntag)
+/* Read numread values from FILE f and put the values in ntag */
+/* Return 0 if successful, 1 if not */
+/* The storage ntag[1..numread] must exist */
+/* Note: the first 8 bytes of the tag file contain the number of particles
+and the number of groups. These must be skipped before calling this routine. */
+{
+ if (fread(ntag+1, sizeof(int), numread, f)!=numread)
+ myerror("Error in reading tag file.");
+ return 0;
+}
+
+int skiptagheader(FILE *f, Slice *s)
+/* Read the first 8 bytes from the tag file. Check that the first int equals
+the number of particles. Return the second, which is the number of groups */
+{
+ int dummy[2];
+ if (fread(&dummy, sizeof(int), 2, f)!=2) myerror("Error in reading tag file.");
+ if (s->numpart!=0 && dummy[0]!=s->numpart)
+ myerror("First number in tag file doesn't match expected number of particles.");
+ s->numgroups = dummy[1];
+ return dummy[1];
+}
+
+int readalltags(FILE *f, Slice *s)
+/* Read the whole tag file. Allocate memory as needed */
+/* Return the number of groups */
+{
+ int dummy[2];
+ if (s->ntag==NULL || s->numlist!=s->numpart) {
+ if (s->ntag!=NULL) {
+ free_ivector(s->ntag, 1, s->numlist);
+ free_ivector(s->ID, 1, s->numlist); /* S Skory */
+ }
+ s->numlist = s->numpart;
+ s->ntag = ivector(1, s->numlist);
+ s->ID = ivector(1, s->numlist);
+ }
+ if (fread(&dummy, sizeof(int), 2, f)!=2) myerror("Error 1 in reading tag file.");
+ if (dummy[0]!=s->numpart)
+ myerror("First int of tag file doesn't match numpart.");
+ s->numgroups = dummy[1];
+
+ if (fread(s->ntag+1, sizeof(int), s->numlist, f)!=s->numlist)
+ myerror("Couldn't read entire tag file.");
+ return dummy[1];
+}
+
+#endif
+
+/* ===================== Warnings and Errors =========================== */
+
+/* Print a message and die */
+void myerror(char *message)
+{
+ fprintf(stderr, "%s\n", message);
+ exit(1); return;
+}
+
+/* Just print a message */
+void mywarn(char *message)
+{
+ fprintf(stderr, "%s\n", message);
+ fflush(NULL); /* Flush everything, so we know where we are */
+ return;
+}
Added: trunk/yt/lagos/hop/hop_smooth.c
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/hop_smooth.c Tue Jun 10 09:05:10 2008
@@ -0,0 +1,465 @@
+/* SMOOTH.C */
+/* This was written by Joachim Stadel and the NASA HPCC ESS at
+the University of Washington Department of Astronomy as part of
+the SMOOTH program, v2.0.1.
+URL: http://www-hpcc.astro.washington.edu/tools/SMOOTH */
+
+/* DJE--I have removed unneeded subroutines, notably those having
+to do with velocity field reconstructions (because they refer to
+particle data that I chose not to store) and output routines
+(because I wanted binary output). Also, the density subroutine
+was slightly customized to reduce memory consumption in
+the case of equal mass particles. */
+
+/* HOP Version 1.0 (12/15/97) -- Original Release */
+
+#include <stdio.h>
+//#include <malloc.h>
+#include <math.h>
+#include <assert.h>
+#include "smooth.h"
+#include "kd.h"
+
+#define ISYM "d"
+#define GSYM "g"
+
+//#include "macros_and_parameters.h"
+
+#define IMARK 1 /* All particles are marked to be included */
+
+int smInit(SMX *psmx,KD kd,int nSmooth,float *fPeriod)
+{
+ SMX smx;
+ PQ_STATIC;
+ int pi,j;
+ fprintf(stderr,"nSmooth = %d kd->nActive = %d\n", nSmooth, kd->nActive);
+ assert(nSmooth <= kd->nActive);
+ smx = (SMX)malloc(sizeof(struct smContext));
+ assert(smx != NULL);
+ smx->kd = kd;
+ smx->nSmooth = nSmooth;
+ smx->pq = (PQ *)malloc(nSmooth*sizeof(PQ));
+ assert(smx->pq != NULL);
+ PQ_INIT(smx->pq,nSmooth);
+ smx->pfBall2 = (float *)malloc((kd->nActive+1)*sizeof(int));
+ assert(smx->pfBall2 != NULL);
+ smx->iMark = (char *)malloc(kd->nActive*sizeof(char));
+ assert(smx->iMark);
+ smx->nListSize = smx->nSmooth+RESMOOTH_SAFE;
+ smx->fList = (float *)malloc(smx->nListSize*sizeof(float));
+ assert(smx->fList != NULL);
+ smx->pList = (int *)malloc(smx->nListSize*sizeof(int));
+ assert(smx->pList != NULL);
+ /*
+ ** Set for Periodic Boundary Conditions.
+ */
+ for (j=0;j<3;++j) smx->fPeriod[j] = fPeriod[j];
+ /*
+ ** Initialize arrays for calculated quantities.--DJE
+ */
+ for (pi=0;pi<smx->kd->nActive;++pi) {
+ smx->kd->p[pi].fDensity = 0.0;
+ smx->kd->p[pi].iHop = 0;
+ }
+ *psmx = smx;
+ return(1);
+ }
+
+
+void smFinish(SMX smx)
+{
+ free(smx->pfBall2);
+ free(smx->iMark);
+ free(smx->pq);
+ free(smx);
+ }
+
+
+void smBallSearch(SMX smx,float fBall2,float *ri)
+{
+ KDN *c;
+ PARTICLE *p;
+ int cell,cp,ct,pj;
+ float fDist2,dx,dy,dz,lx,ly,lz,sx,sy,sz,x,y,z;
+ PQ *pq;
+ PQ_STATIC;
+
+ c = smx->kd->kdNodes;
+ p = smx->kd->p;
+ pq = smx->pqHead;
+ x = ri[0];
+ y = ri[1];
+ z = ri[2];
+ lx = smx->fPeriod[0];
+ ly = smx->fPeriod[1];
+ lz = smx->fPeriod[2];
+ cell = ROOT;
+ /*
+ ** First find the "local" Bucket.
+ ** This could mearly be the closest bucket to ri[3].
+ */
+ while (cell < smx->kd->nSplit) {
+ if (ri[c[cell].iDim] < c[cell].fSplit) cell = LOWER(cell);
+ else cell = UPPER(cell);
+ }
+ /*
+ ** Now start the search from the bucket given by cell!
+ */
+ for (pj=c[cell].pLower;pj<=c[cell].pUpper;++pj) {
+ dx = x - p[pj].r[0];
+ dy = y - p[pj].r[1];
+ dz = z - p[pj].r[2];
+ fDist2 = dx*dx + dy*dy + dz*dz;
+ if (fDist2 < fBall2) {
+ if (smx->iMark[pj]) continue;
+ smx->iMark[pq->p] = 0;
+ smx->iMark[pj] = 1;
+ pq->fKey = fDist2;
+ pq->p = pj;
+ pq->ax = 0.0;
+ pq->ay = 0.0;
+ pq->az = 0.0;
+ PQ_REPLACE(pq);
+ fBall2 = pq->fKey;
+ }
+ }
+ while (cell != ROOT) {
+ cp = SIBLING(cell);
+ ct = cp;
+ SETNEXT(ct);
+ while (1) {
+ INTERSECT(c,cp,fBall2,lx,ly,lz,x,y,z,sx,sy,sz);
+ /*
+ ** We have an intersection to test.
+ */
+ if (cp < smx->kd->nSplit) {
+ cp = LOWER(cp);
+ continue;
+ }
+ else {
+ for (pj=c[cp].pLower;pj<=c[cp].pUpper;++pj) {
+ dx = sx - p[pj].r[0];
+ dy = sy - p[pj].r[1];
+ dz = sz - p[pj].r[2];
+ fDist2 = dx*dx + dy*dy + dz*dz;
+ if (fDist2 < fBall2) {
+ if (smx->iMark[pj]) continue;
+ smx->iMark[pq->p] = 0;
+ smx->iMark[pj] = 1;
+ pq->fKey = fDist2;
+ pq->p = pj;
+ pq->ax = sx - x;
+ pq->ay = sy - y;
+ pq->az = sz - z;
+ PQ_REPLACE(pq);
+ fBall2 = pq->fKey;
+ }
+ }
+ }
+ GetNextCell:
+ SETNEXT(cp);
+ if (cp == ct) break;
+ }
+ cell = PARENT(cell);
+ }
+ smx->pqHead = pq;
+ }
+
+
+int smBallGather(SMX smx,float fBall2,float *ri)
+{
+ KDN *c;
+ PARTICLE *p;
+ int pj,nCnt,cp,nSplit;
+ float dx,dy,dz,x,y,z,lx,ly,lz,sx,sy,sz,fDist2;
+
+ c = smx->kd->kdNodes;
+ p = smx->kd->p;
+ nSplit = smx->kd->nSplit;
+ lx = smx->fPeriod[0];
+ ly = smx->fPeriod[1];
+ lz = smx->fPeriod[2];
+ x = ri[0];
+ y = ri[1];
+ z = ri[2];
+ nCnt = 0;
+ cp = ROOT;
+ while (1) {
+ INTERSECT(c,cp,fBall2,lx,ly,lz,x,y,z,sx,sy,sz);
+ /*
+ ** We have an intersection to test.
+ */
+ if (cp < nSplit) {
+ cp = LOWER(cp);
+ continue;
+ }
+ else {
+ for (pj=c[cp].pLower;pj<=c[cp].pUpper;++pj) {
+ dx = sx - p[pj].r[0];
+ dy = sy - p[pj].r[1];
+ dz = sz - p[pj].r[2];
+ fDist2 = dx*dx + dy*dy + dz*dz;
+ if (fDist2 < fBall2) {
+ smx->fList[nCnt] = fDist2;
+ smx->pList[nCnt++] = pj;
+ /* Insert debugging flag here */
+ if (nCnt > smx->nListSize) {
+ fprintf(stderr,"nCnt too big.\n");
+ }
+ }
+ }
+ }
+ GetNextCell:
+ SETNEXT(cp);
+ if (cp == ROOT) break;
+ }
+ assert(nCnt <= smx->nListSize);
+ return(nCnt);
+ }
+
+
+void smSmooth(SMX smx,void (*fncSmooth)(SMX,int,int,int *,float *))
+{
+ KDN *c;
+ PARTICLE *p;
+ PQ *pq,*pqLast;
+ PQ_STATIC;
+ int cell;
+ int pi,pin,pj,pNext,nCnt,nSmooth;
+ float dx,dy,dz,x,y,z,h2,ax,ay,az;
+
+
+ for (pi=0;pi<smx->kd->nActive;++pi) {
+ if (IMARK) smx->pfBall2[pi] = -1.0;
+ else smx->pfBall2[pi] = 1.0; /* pretend it is already done! */
+ }
+ smx->pfBall2[smx->kd->nActive] = -1.0; /* stop condition */
+ for (pi=0;pi<smx->kd->nActive;++pi) {
+ smx->iMark[pi] = 0;
+ }
+ pqLast = &smx->pq[smx->nSmooth-1];
+ c = smx->kd->kdNodes;
+ p = smx->kd->p;
+ nSmooth = smx->nSmooth;
+ /*
+ ** Initialize Priority Queue.
+ */
+ pin = 0;
+ pNext = 1;
+ ax = 0.0;
+ ay = 0.0;
+ az = 0.0;
+ for (pq=smx->pq,pj=0;pq<=pqLast;++pq,++pj) {
+ smx->iMark[pj] = 1;
+ pq->p = pj;
+ pq->ax = ax;
+ pq->ay = ay;
+ pq->az = az;
+ }
+ while (1) {
+ if (smx->pfBall2[pin] >= 0) {
+ /*
+ ** Find next particle which is not done, and load the
+ ** priority queue with nSmooth number of particles.
+ */
+ while (smx->pfBall2[pNext] >= 0) ++pNext;
+ /*
+ ** Check if we are really finished.
+ */
+ if (pNext == smx->kd->nActive) break;
+ pi = pNext;
+ ++pNext;
+ x = p[pi].r[0];
+ y = p[pi].r[1];
+ z = p[pi].r[2];
+ /* printf("%"ISYM": %"GSYM" %"GSYM" %"GSYM"\n", pi, x, y, z); */
+ /*
+ ** First find the "local" Bucket.
+ ** This could mearly be the closest bucket to ri[3].
+ */
+ cell = ROOT;
+ while (cell < smx->kd->nSplit) {
+ if (p[pi].r[c[cell].iDim] < c[cell].fSplit)
+ cell = LOWER(cell);
+ else
+ cell = UPPER(cell);
+ }
+ /*
+ ** Remove everything from the queue.
+ */
+ smx->pqHead = NULL;
+ for (pq=smx->pq;pq<=pqLast;++pq) smx->iMark[pq->p] = 0;
+ /*
+ ** Add everything from pj up to and including pj+nSmooth-1.
+ */
+ pj = c[cell].pLower;
+ if (pj > smx->kd->nActive - nSmooth)
+ pj = smx->kd->nActive - nSmooth;
+ for (pq=smx->pq;pq<=pqLast;++pq) {
+ smx->iMark[pj] = 1;
+ dx = x - p[pj].r[0];
+ dy = y - p[pj].r[1];
+ dz = z - p[pj].r[2];
+ pq->fKey = dx*dx + dy*dy + dz*dz;
+ pq->p = pj++;
+ pq->ax = 0.0;
+ pq->ay = 0.0;
+ pq->az = 0.0;
+ }
+ PQ_BUILD(smx->pq,nSmooth,smx->pqHead);
+ }
+ else {
+ /*
+ ** Calculate the priority queue using the previous particles!
+ */
+ pi = pin;
+ x = p[pi].r[0];
+ y = p[pi].r[1];
+ z = p[pi].r[2];
+ smx->pqHead = NULL;
+ for (pq=smx->pq;pq<=pqLast;++pq) {
+ pq->ax -= ax;
+ pq->ay -= ay;
+ pq->az -= az;
+ dx = x + pq->ax - p[pq->p].r[0];
+ dy = y + pq->ay - p[pq->p].r[1];
+ dz = z + pq->az - p[pq->p].r[2];
+ pq->fKey = dx*dx + dy*dy + dz*dz;
+ }
+ PQ_BUILD(smx->pq,nSmooth,smx->pqHead);
+ ax = 0.0;
+ ay = 0.0;
+ az = 0.0;
+ }
+ smBallSearch(smx,smx->pqHead->fKey,p[pi].r);
+ smx->pfBall2[pi] = smx->pqHead->fKey;
+ /*
+ ** Pick next particle, 'pin'.
+ ** Create fList and pList for function 'fncSmooth'.
+ */
+ pin = pi;
+ nCnt = 0;
+ h2 = smx->pqHead->fKey;
+ for (pq=smx->pq;pq<=pqLast;++pq) {
+ if (pq == smx->pqHead) continue;
+ smx->pList[nCnt] = pq->p;
+ smx->fList[nCnt++] = pq->fKey;
+ if (smx->pfBall2[pq->p] >= 0) continue;
+ if (pq->fKey < h2) {
+ pin = pq->p;
+ h2 = pq->fKey;
+ ax = pq->ax;
+ ay = pq->ay;
+ az = pq->az;
+ }
+ }
+ (*fncSmooth)(smx,pi,nCnt,smx->pList,smx->fList);
+ }
+ }
+
+
+void smReSmooth(SMX smx,void (*fncSmooth)(SMX,int,int,int *,float *))
+{
+ PARTICLE *p;
+ int pi,nSmooth;
+
+ p = smx->kd->p;
+ for (pi=0;pi<smx->kd->nActive;++pi) {
+ if (IMARK == 0) continue;
+ /*
+ ** Do a Ball Gather at the radius of the most distant particle
+ ** which is smDensity sets in smx->pBall[pi].
+ */
+ nSmooth = smBallGather(smx,smx->pfBall2[pi],p[pi].r);
+ (*fncSmooth)(smx,pi,nSmooth,smx->pList,smx->fList);
+ }
+ }
+
+
+void smDensity(SMX smx,int pi,int nSmooth,int *pList,float *fList)
+{
+ float ih2,r2,rs,fDensity;
+ int i,pj;
+
+ ih2 = 4.0/smx->pfBall2[pi];
+ fDensity = 0.0;
+ for (i=0;i<nSmooth;++i) {
+ pj = pList[i];
+ r2 = fList[i]*ih2;
+ rs = 2.0 - sqrt(r2);
+ if (r2 < 1.0) rs = (1.0 - 0.75*rs*r2);
+ else rs = 0.25*rs*rs*rs;
+#ifdef DIFFERENT_MASSES
+ fDensity += rs*smx->kd->p[pj].fMass;
+#else
+ fDensity += rs*smx->kd->fMass;
+#endif
+ }
+ smx->kd->p[pi].fDensity = M_1_PI*sqrt(ih2)*ih2*fDensity;
+ }
+
+
+void smDensitySym(SMX smx,int pi,int nSmooth,int *pList,float *fList)
+{
+ float fNorm,ih2,r2,rs;
+ int i,pj;
+
+ ih2 = 4.0/smx->pfBall2[pi];
+ fNorm = 0.5*M_1_PI*sqrt(ih2)*ih2;
+ for (i=0;i<nSmooth;++i) {
+ pj = pList[i];
+ r2 = fList[i]*ih2;
+ rs = 2.0 - sqrt(r2);
+ if (r2 < 1.0) rs = (1.0 - 0.75*rs*r2);
+ else rs = 0.25*rs*rs*rs;
+ rs *= fNorm;
+#ifdef DIFFERENT_MASSES
+ smx->kd->p[pi].fDensity += rs*smx->kd->p[pj].fMass;
+ smx->kd->p[pj].fDensity += rs*smx->kd->p[pi].fMass;
+#else
+ smx->kd->p[pi].fDensity += rs*smx->kd->fMass;
+ smx->kd->p[pj].fDensity += rs*smx->kd->fMass;
+#endif
+ }
+ }
+
+/* I'm not using the following function, but I left it here in case someone
+wants the densities outputted in Tipsy format. But you're probably better
+off just fetching the smooth() program from the HPCC web site... */
+
+void smOutDensity(SMX smx,FILE *fp)
+{
+ int i,iCnt;
+
+ fprintf(fp,"%"ISYM"\n",smx->kd->nParticles);
+ iCnt = 0;
+ for (i=0;i<smx->kd->nGas;++i) {
+ if (smx->kd->bGas) {
+ if (IMARK)
+ fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ else fprintf(fp,"0\n");
+ ++iCnt;
+ }
+ else fprintf(fp,"0\n");
+ }
+ for (i=0;i<smx->kd->nDark;++i) {
+ if (smx->kd->bDark) {
+ if (IMARK)
+ fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ else fprintf(fp,"0\n");
+ ++iCnt;
+ }
+ else fprintf(fp,"0\n");
+ }
+ for (i=0;i<smx->kd->nStar;++i) {
+ if (smx->kd->bStar) {
+ if (IMARK)
+ fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ else fprintf(fp,"0\n");
+ ++iCnt;
+ }
+ else fprintf(fp,"0\n");
+ }
+}
+
+
Added: trunk/yt/lagos/hop/kd.h
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/kd.h Tue Jun 10 09:05:10 2008
@@ -0,0 +1,192 @@
+/* KD.H */
+/* This was written by Joachim Stadel and the NASA HPCC ESS at
+the University of Washington Department of Astronomy as part of
+the SMOOTH program, v2.0.1.
+URL: http://www-hpcc.astro.washington.edu/tools/SMOOTH */
+
+/* DJE--I have made a few alterations to the PARTICLE structure
+in order to reduce memory consumption. */
+
+/* HOP Version 1.0 (12/15/97) -- Original Release */
+
+/* GLB--set different masses on */
+
+#define DIFFERENT_MASSES
+
+//#include "macros_and_parameters.h"
+
+#ifndef KD_HINCLUDED
+#define KD_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)\
+{\
+ while (i&1) i=i>>1;\
+ ++i;\
+ }
+
+#define DARK 1
+#define GAS 2
+#define STAR 4
+
+typedef struct Particle {
+ float r[3];
+ int iOrder;
+ float fDensity;
+ int iID; /* the real ID of the particle S. Skory */
+ int iHop; /* DJE: The number of the highest-density neighbor;
+ Later, the group number. */
+#ifdef DIFFERENT_MASSES
+ float fMass;
+#endif
+ /* DJE: The following are unused and cost too much memory to keep */
+ /* float v[3]; */
+ /* float fMass; */
+ /* int iMark; */
+ /* float vMean[3]; */
+ /* float fVelDisp2; */
+ } PARTICLE;
+
+typedef struct bndBound {
+ float fMin[3];
+ float fMax[3];
+ } BND;
+
+typedef struct kdNode {
+ float fSplit;
+ BND bnd;
+ int iDim;
+ int pLower;
+ int pUpper;
+ } KDN;
+
+typedef struct kdContext {
+ int nBucket;
+ int nParticles;
+ int nDark;
+ int nGas;
+ int nStar;
+ int bDark;
+ int bGas;
+ int bStar;
+ int nActive;
+ float fTime;
+ BND bnd;
+ int nLevels;
+ int nNodes;
+ int nSplit;
+ float fMass; /* DJE: If all particles have the same mass */
+ PARTICLE *p;
+ KDN *kdNodes;
+ int uSecond;
+ int uMicro;
+ } * KD;
+
+
+#define INTERSECT(c,cp,fBall2,lx,ly,lz,x,y,z,sx,sy,sz)\
+{\
+ float dx,dy,dz,dx1,dy1,dz1,fDist2;\
+ dx = c[cp].bnd.fMin[0]-x;\
+ dx1 = x-c[cp].bnd.fMax[0];\
+ dy = c[cp].bnd.fMin[1]-y;\
+ dy1 = y-c[cp].bnd.fMax[1];\
+ dz = c[cp].bnd.fMin[2]-z;\
+ dz1 = z-c[cp].bnd.fMax[2];\
+ if (dx > 0.0) {\
+ dx1 += lx;\
+ if (dx1 < dx) {\
+ fDist2 = dx1*dx1;\
+ sx = x+lx;\
+ }\
+ else {\
+ fDist2 = dx*dx;\
+ sx = x;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else if (dx1 > 0.0) {\
+ dx += lx;\
+ if (dx < dx1) {\
+ fDist2 = dx*dx;\
+ sx = x-lx;\
+ }\
+ else {\
+ fDist2 = dx1*dx1;\
+ sx = x;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else {\
+ fDist2 = 0.0;\
+ sx = x;\
+ }\
+ if (dy > 0.0) {\
+ dy1 += ly;\
+ if (dy1 < dy) {\
+ fDist2 += dy1*dy1;\
+ sy = y+ly;\
+ }\
+ else {\
+ fDist2 += dy*dy;\
+ sy = y;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else if (dy1 > 0.0) {\
+ dy += ly;\
+ if (dy < dy1) {\
+ fDist2 += dy*dy;\
+ sy = y-ly;\
+ }\
+ else {\
+ fDist2 += dy1*dy1;\
+ sy = y;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else {\
+ sy = y;\
+ }\
+ if (dz > 0.0) {\
+ dz1 += lz;\
+ if (dz1 < dz) {\
+ fDist2 += dz1*dz1;\
+ sz = z+lz;\
+ }\
+ else {\
+ fDist2 += dz*dz;\
+ sz = z;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else if (dz1 > 0.0) {\
+ dz += lz;\
+ if (dz < dz1) {\
+ fDist2 += dz*dz;\
+ sz = z-lz;\
+ }\
+ else {\
+ fDist2 += dz1*dz1;\
+ sz = z;\
+ }\
+ if (fDist2 > fBall2) goto GetNextCell;\
+ }\
+ else {\
+ sz = z;\
+ }\
+ }
+
+
+void kdTime(KD,int *,int *);
+int kdInit(KD *,int);
+int kdReadTipsy(KD,FILE *,int,int,int);
+void kdInMark(KD,char *);
+int kdBuildTree(KD);
+void kdOrder(KD);
+void kdFinish(KD);
+
+#endif
Added: trunk/yt/lagos/hop/setup.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/setup.py Tue Jun 10 09:05:10 2008
@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+import setuptools
+import os, sys, os.path
+
+import os.path
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('hop',parent_package,top_path)
+ config.make_config_py() # installs __config__.py
+ config.add_extension("EnzoHop", sources=
+ ["EnzoHop.c",
+ "hop_hop.c",
+ "hop_kd.c",
+ "hop_regroup.c",
+ "hop_slice.c",
+ "hop_smooth.c",])
+ return config
Added: trunk/yt/lagos/hop/slice.h
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/slice.h Tue Jun 10 09:05:10 2008
@@ -0,0 +1,97 @@
+/* SLICE.H, Daniel Eisenstein, 1997 */
+/* Based on a paper by Daniel Eisenstein & Piet Hut,
+"HOP: A New Group-Finding Algorithm for N-body Simulations."
+See the included documentation or view it at
+http://www.sns.ias.edu/~eisenste/hop/hop_doc.html */
+
+/* Version 1.0 (12/15/97) -- Original Release */
+
+#ifndef NBODYUTIL_H
+#define NBODYUTIL_H
+
+#define RHOCRIT 277.5 /* in h^2 Msun/kpc^3 */
+#define HTIME 9.7776 /* in h^-1 Gyr */
+#define GNEWT 4.51e-6 /* in kpc^3/Msun/Gyr^2, h cancels */
+#define KMS 0.975 /* to convert from kpc/Gyr to km/s */
+
+#define PI 3.141592654
+#define ROOT2 1.414213562
+#define ROOTPI2 1.253314137 /* sqrt(Pi/2) */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+//#include "macros_and_parameters.h"
+
+/* Structure to hold data from one time slice and info about the slice */
+/* This needn't hold all of the points; it might just hold some subset */
+/* Hence, data vectors are just given as arrays here */
+/* Usage: pid[] Array goes from 0 to pid[0]. pid is NULL if the particle
+numbering is just the trivial mapping. Kinematic arrays go from 1 to numlist
+and can be NULL if the data hasn't been read. Numlist==pid[0] often */
+
+typedef struct slicestruct {
+#ifdef NOT_USED
+ /* First, some generic stuff about the simulation */
+ float omega, lambda, curv;
+ float h0; /* H_0 = 100 km/s/Mpc * h0 */
+ float specn, gamma; /* The spectral index and BBKS PS Gamma */
+ float sigma8; /* At z=0 */
+
+ /* Next, some information about this slice in particular */
+ float z, a, t, growth;
+ /* a=1 at z=0, g=a at early times */
+ float hubb; /* This is a*H(z), but without h0. So it's 0.1 km/s/kpc
+ redshifted appropriately. This is used to relate
+ comoving positions and peculiar velocities */
+
+ /* Now some information about the data */
+ int numblocks; /* Number of blocks in the data file */
+ int numperblock; /* Number of particles per block */
+ float masspart; /* Mass per particle in h^-1 Msun */
+ float boxsize; /* Comoving size of box in h^-1 kpc */
+ float physsize; /* Physical Size in h^-1 kpc */
+ float velscale; /* To turn raw velocity data into peculiar vel in km/s*/
+#endif
+
+ int numpart; /* Total number of particles in the simulation */
+ /* Now the data itself */
+ int *pid; /* The id number of the particle */
+ /* pid[0] holds the number of particles in the list */
+ int offset; /* If pid==NULL, then the arrays are consecutively
+ numbered, starting from offset+1 */
+ int numlist; /* Length of arrays below, set when allocated */
+ float *px, *py, *pz, *vx, *vy, *vz; /* The kinematic information */
+
+ /* And here's the group tag information */
+ int *ntag; /* Only stored for the numlist above */
+ int *ID; /* The real, true ID of the particle. S Skory */
+ int numgroups; /* The number of groups read out of the tag file */
+} Slice; /* Type Slice is defined */
+
+/* Prototypes */
+Slice *newslice();
+void free_tags(Slice *s);
+void free_data(Slice *s);
+void free_slice(Slice *s);
+int f77write(FILE *f, void *p, int len);
+int f77read(FILE *f, void *p, int len);
+float *vector(long nl, long nh);
+int *ivector(long nl, long nh);
+void free_vector(float *v, long nl, long nh);
+void free_ivector(int *v, long nl, long nh);
+
+void myerror(char *message);
+void mywarn(char *message);
+
+int read_header(FILE *f, Slice *s);
+void normalizedata(Slice *s, int conp, int conv);
+int read_alldata(FILE *f, FILE *ftag, Slice *s, int conp, int conv);
+int read_partdata(FILE *f, FILE *ftag, Slice *s);
+
+int readtag(FILE *f, int numread, int *ntag);
+int skiptagheader(FILE *f, Slice *s);
+int readalltags(FILE *f, Slice *s);
+
+
+#endif /* NBODYUTIL_H */
Added: trunk/yt/lagos/hop/smooth.h
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/smooth.h Tue Jun 10 09:05:10 2008
@@ -0,0 +1,140 @@
+/* SMOOTH.H */
+/* This was written by Joachim Stadel and the NASA HPCC ESS at
+the University of Washington Department of Astronomy as part of
+the SMOOTH program, v2.0.1.
+URL: http://www-hpcc.astro.washington.edu/tools/SMOOTH */
+
+/* DJE--I have made a few additions to the SMX structure
+in order to store information necessary for HOP. I have also
+added the Boundary structure. */
+
+/* HOP Version 1.0 (12/15/97) -- Original Release */
+
+#ifndef SMOOTH_HINCLUDED
+#define SMOOTH_HINCLUDED
+
+#include "kd.h"
+//#include "macros_and_parameters.h"
+
+#define RESMOOTH_SAFE 30
+
+/* DJE: Define this structure to hold the boundary data. */
+typedef struct boundarystruct {
+ int nGroup1, nGroup2; /* The two groups involved, ordered such
+ that nGroup1<nGroup2 */
+ float fDensity; /* The highest average density of two
+ boundary particles */
+ } Boundary; /* Type Boundary is defined */
+
+
+typedef struct pqNode {
+ float fKey;
+ struct pqNode *pqLoser;
+ struct pqNode *pqFromInt;
+ struct pqNode *pqFromExt;
+ struct pqNode *pqWinner; /* Only used when building initial tree */
+ int p;
+ float ax;
+ float ay;
+ float az;
+ } PQ;
+
+
+typedef struct smContext {
+ KD kd;
+ int nSmooth;
+ float fPeriod[3];
+ PQ *pq;
+ PQ *pqHead;
+ float *pfBall2;
+ char *iMark;
+ int nListSize;
+ float *fList;
+ int *pList;
+ /* DJE -- Added the following fields to SMX */
+ int nDens; /* The number of neighbors for calculating density */
+ int nHop; /* How many neighbors to consider when hopping*/
+ int nMerge; /* The # of neighbors to look at when merging */
+ int nGroups; /* Groups number 0 to nGroups - 1 */
+ int *nmembers; /* The number of particles in the groups */
+ /* nmembers[nGroups] is the # not in groups */
+ int *densestingroup; /* The ID # of the densest particle in the gr */
+ int nHashLength; /* The length of the hash table */
+ Boundary *hash; /* The hash table for boundaries */
+ float fDensThresh; /* Density Threshold for group finding */
+ } * SMX;
+
+#define PQ_STATIC PQ *PQ_t,*PQ_lt;int PQ_j,PQ_i
+
+
+#define PQ_INIT(pq,n)\
+{\
+ for (PQ_j=0;PQ_j<(n);++PQ_j) {\
+ if (PQ_j < 2) (pq)[PQ_j].pqFromInt = NULL;\
+ else (pq)[PQ_j].pqFromInt = &(pq)[PQ_j>>1];\
+ (pq)[PQ_j].pqFromExt = &(pq)[(PQ_j+(n))>>1];\
+ }\
+ }
+
+
+#define PQ_BUILD(pq,n,q)\
+{\
+ for (PQ_j=(n)-1;PQ_j>0;--PQ_j) {\
+ PQ_i = (PQ_j<<1);\
+ if (PQ_i < (n)) PQ_t = (pq)[PQ_i].pqWinner;\
+ else PQ_t = &(pq)[PQ_i-(n)];\
+ ++PQ_i;\
+ if (PQ_i < (n)) PQ_lt = (pq)[PQ_i].pqWinner;\
+ else PQ_lt = &(pq)[PQ_i-(n)];\
+ if (PQ_t->fKey < PQ_lt->fKey) {\
+ (pq)[PQ_j].pqLoser = PQ_t;\
+ (pq)[PQ_j].pqWinner = PQ_lt;\
+ }\
+ else {\
+ (pq)[PQ_j].pqLoser = PQ_lt;\
+ (pq)[PQ_j].pqWinner = PQ_t;\
+ }\
+ }\
+ (q) = (pq)[1].pqWinner;\
+ }
+
+
+#define PQ_REPLACE(q)\
+{\
+ PQ_t = (q)->pqFromExt;\
+ while (PQ_t) {\
+ if (PQ_t->pqLoser->fKey > (q)->fKey) {\
+ PQ_lt = PQ_t->pqLoser;\
+ PQ_t->pqLoser = (q);\
+ (q) = PQ_lt;\
+ }\
+ PQ_t = PQ_t->pqFromInt;\
+ }\
+ }
+
+
+
+int smInit(SMX *,KD,int,float *);
+void smFinish(SMX);
+void smBallSearch(SMX,float,float *);
+int smBallGather(SMX,float,float *);
+void smSmooth(SMX,void (*)(SMX,int,int,int *,float *));
+void smReSmooth(SMX,void (*)(SMX,int,int,int *,float *));
+void smDensity(SMX,int,int,int *,float *);
+void smDensitySym(SMX,int,int,int *,float *);
+void smMeanVel(SMX,int,int,int *,float *);
+void smMeanVelSym(SMX,int,int,int *,float *);
+void smVelDisp(SMX,int,int,int *,float *);
+void smVelDispSym(SMX,int,int,int *,float *);
+void smNull(SMX,int,int,int *,float *);
+void smOutDensity(SMX,FILE *);
+void smOutMeanVel(SMX,FILE *);
+void smOutVelDisp(SMX,FILE *);
+void smOutPhase(SMX,FILE *);
+void smOutMach(SMX,FILE *);
+void smOutSpeed(SMX,FILE *);
+
+#endif
+
+
+
Modified: trunk/yt/lagos/setup.py
==============================================================================
--- trunk/yt/lagos/setup.py (original)
+++ trunk/yt/lagos/setup.py Tue Jun 10 09:05:10 2008
@@ -6,14 +6,17 @@
try:
H5dir = open("hdf5.cfg").read().strip().rstrip()
except:
- print "Reading HDF5 location from hdf5.cfg failed: Not using HDF5 wrapper!"
- H5dir = None
+ print "Reading HDF5 location from hdf5.cfg failed."
+ print "Please place the base directory of your HDF5 install in hdf5.cfg and restart."
+ print "(ex: \"echo '/usr/local/' > hdf5.cfg\" )"
+ sys.exit(1)
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('lagos',parent_package,top_path)
config.make_config_py() # installs __config__.py
config.add_extension("PointCombine", "yt/lagos/PointCombine.c", libraries=["m"])
+ config.add_subpackage("hop")
if H5dir is not None:
include_dirs=[os.path.join(H5dir,"include")]
library_dirs=[os.path.join(H5dir,"lib")]
More information about the yt-svn
mailing list