[Yt-svn] yt-commit r551 - trunk/yt/lagos/hop
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Jun 12 16:51:17 PDT 2008
Author: mturk
Date: Thu Jun 12 16:51:16 2008
New Revision: 551
URL: http://yt.spacepope.org/changeset/551
Log:
It works and reproduces existing HOP output. More tests are needed, and I need
to add the HopAnalysis.out output. Otherwise, I am pretty pleased with where
it's at.
To test it, for example:
>>> pf = lagos.EnzoStaticOutput("my_data.dir/my_data")
>>> sp = pf.h.sphere([0.5]*3,1.0)
>>> hop = yt.lagos.hop.HopList(sp, 80.0)
>>> for i,group in enumerate(hop):
... md = group.maximum_density_location()
... print "% 5i\t%0.9e\t% 8i\t%0.8f %0.8f %0.8f" % \
... (i, group.total_mass()*m_conv, len(group.indices), \
... md[0], md[1], md[2])
Added:
trunk/yt/lagos/hop/HopOutput.py
Modified:
trunk/yt/lagos/hop/EnzoHop.c
trunk/yt/lagos/hop/__init__.py
trunk/yt/lagos/hop/hop_hop.c
trunk/yt/lagos/hop/hop_regroup.c
Modified: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- trunk/yt/lagos/hop/EnzoHop.c (original)
+++ trunk/yt/lagos/hop/EnzoHop.c Thu Jun 12 16:51:16 2008
@@ -48,11 +48,12 @@
*mass, *ID;
xpos=ypos=zpos=mass=ID=NULL;
npy_float64 totalmass = 0;
+ float thresh = 160.0;
int i;
- if (!PyArg_ParseTuple(args, "OOOOO",
- &oxpos, &oypos, &ozpos, &omass, &oID))
+ if (!PyArg_ParseTuple(args, "OOOOO|f",
+ &oxpos, &oypos, &ozpos, &omass, &oID, &thresh))
return PyErr_Format(_HOPerror,
"EnzoHop: Invalid parameters.");
@@ -135,7 +136,6 @@
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");
Added: trunk/yt/lagos/hop/HopOutput.py
==============================================================================
--- (empty file)
+++ trunk/yt/lagos/hop/HopOutput.py Thu Jun 12 16:51:16 2008
@@ -0,0 +1,135 @@
+"""
+HOP-output data handling
+
+ at author: U{Matthew Turk<http://www.stanford.edu/~mturk/>}
+ at organization: U{KIPAC<http://www-group.slac.stanford.edu/KIPAC/>}
+ at contact: U{mturk at slac.stanford.edu<mailto:mturk at slac.stanford.edu>}
+ at license:
+ 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/>.
+"""
+
+from yt.lagos.hop import *
+
+class HopList(object):
+ def __init__(self, data_source, threshold=160.0,
+ dm_only = True):
+ self.data_source = data_source
+ self.dm_only = dm_only
+ self.threshold = threshold
+ self._groups = []
+ self._max_dens = {}
+ mylog.info("Initializing HOP")
+ self.__run_hop()
+ mylog.info("Parsing outputs")
+ self.__parse_output()
+ mylog.debug("Finished.")
+
+ def __run_hop(self):
+ if self.dm_only: ii = self.__get_dm_indices()
+ else: ii = slice(None)
+ self.densities, self.tags = \
+ RunHOP(self.data_source["particle_position_x"][ii],
+ self.data_source["particle_position_y"][ii],
+ self.data_source["particle_position_z"][ii],
+ self.data_source["ParticleMassMsun"][ii],
+ self.data_source["particle_index"][ii].astype('int64'),
+ self.threshold)
+
+ def __get_dm_indices(self):
+ if 'creation_time' in self.data_source.hierarchy.field_list:
+ mylog.debug("Differentiating based on creation time")
+ return (self.data_source["creation_time"] < 0)
+ elif 'particle_type' in self.data_source.hierarchy.field_list:
+ mylog.debug("Differentiating based on particle type")
+ return (self.data_source["particle_type"] == 1)
+ else:
+ raise KeyError
+
+ def __parse_output(self):
+ unique_ids = na.unique(self.tags)
+ bins = na.digitize(self.tags, unique_ids)
+ counts = na.bincount(self.tags+1)
+ sort_indices = na.argsort(self.tags)
+ grab_indices = na.indices(self.tags.shape).ravel()[sort_indices]
+ cp = 0
+ for i in unique_ids:
+ cp_c = cp + counts[i+1]
+ group_indices = grab_indices[cp:cp_c]
+ self._groups.append(HopGroup(self, i, group_indices))
+ md_i = na.argmax(self.densities[sort_indices][cp:cp_c])
+ px, py, pz = [self.data_source['particle_position_%s'%ax][group_indices]
+ for ax in 'xyz']
+ self._max_dens[i] = (self.densities[sort_indices][cp:cp_c][md_i],
+ px[md_i], py[md_i], pz[md_i])
+ cp += counts[i+1]
+
+ def __len__(self):
+ return len(self._groups)
+
+ def __iter__(self):
+ self.__index = 1 # Start with the first non-zero group
+ return self
+
+ def next(self):
+ if self.__index == len(self): raise StopIteration
+ tr = self[self.__index]
+ self.__index += 1
+ return tr
+
+ def __getitem__(self, key):
+ return self._groups[key]
+
+ def write_out(self, filename="HopAnalysis.out"):
+ pass
+
+class HopGroup(object):
+
+ def __init__(self, hop_output, id, indices):
+ self.hop_output = hop_output
+ self.id = id
+ self.data = hop_output.data_source
+ self.indices = indices
+
+ def center_of_mass(self):
+ pm = self["ParticleMassMsun"]
+ cx = (self["particle_position_x"] * pm).sum()
+ cy = (self["particle_position_x"] * pm).sum()
+ cz = (self["particle_position_x"] * pm).sum()
+ return na.array([cx,cy,cz])/pm.sum()
+
+ def maximum_density(self):
+ return self.hop_output._max_dens[self.id][0]
+
+ def maximum_density_location(self):
+ return (self.hop_output._max_dens[self.id][1],
+ self.hop_output._max_dens[self.id][2],
+ self.hop_output._max_dens[self.id][3])
+
+ def total_mass(self):
+ return self["ParticleMassMsun"].sum()
+
+ def bulk_velocity(self):
+ pm = self["ParticleMassMsun"]
+ vx = (self["particle_velocity_x"] * pm).sum()
+ vy = (self["particle_velocity_x"] * pm).sum()
+ vz = (self["particle_velocity_x"] * pm).sum()
+ return na.array([vx,vy,vz])/pm.sum()
+
+
+ def __getitem__(self, key):
+ return self.data[key][self.indices]
Modified: trunk/yt/lagos/hop/__init__.py
==============================================================================
--- trunk/yt/lagos/hop/__init__.py (original)
+++ trunk/yt/lagos/hop/__init__.py Thu Jun 12 16:51:16 2008
@@ -1 +1,4 @@
+from yt.lagos import *
+
from EnzoHop import *
+from HopOutput import *
Modified: trunk/yt/lagos/hop/hop_hop.c
==============================================================================
--- trunk/yt/lagos/hop/hop_hop.c (original)
+++ trunk/yt/lagos/hop/hop_hop.c Thu Jun 12 16:51:16 2008
@@ -47,8 +47,6 @@
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) */
@@ -67,7 +65,8 @@
nSmooth = 64;
nDens = 64;
nHop = -1;
- fDensThresh = -1.0;
+/* fDensThresh = 3.0; */
+ fDensThresh = -1.0;
bDensity = 3;
bGroup = 3;
bMerge = 3;
@@ -103,7 +102,6 @@
PrepareKD(kd);
smInit(&smx,kd,nSmooth,fPeriod);
-
smx->nHop = nHop;
smx->nDens = nDens;
smx->nMerge = nMerge;
@@ -536,38 +534,6 @@
/* ----------------------------------------------------------------- */
-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 */
Modified: trunk/yt/lagos/hop/hop_regroup.c
==============================================================================
--- trunk/yt/lagos/hop/hop_regroup.c (original)
+++ trunk/yt/lagos/hop/hop_regroup.c Thu Jun 12 16:51:16 2008
@@ -233,8 +233,9 @@
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);
+ //if (c.qsort) sort_groups(s, gl, c.mingroupsize, c.outsizename);
+ if (c.qsort) sort_groups(s, gl, c.mingroupsize, NULL);
+ //writegmerge(s, gl, c.outgmergename, c.peak_thresh, c.saddle_thresh);
translatetags(s,gl);
}
else if (c.qgmerge_given) {
@@ -672,8 +673,8 @@
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);
More information about the yt-svn
mailing list