[Yt-svn] yt-commit r1760 - in trunk/yt: . extensions extensions/kdtree lagos lagos/parallelHOP
sskory at wrangler.dreamhost.com
sskory at wrangler.dreamhost.com
Thu Jun 24 17:13:01 PDT 2010
Author: sskory
Date: Thu Jun 24 17:12:57 2010
New Revision: 1760
URL: http://yt.enzotools.org/changeset/1760
Log:
Adding patch from hg, renaming StructureFunctionGenerator->TwoPointFunctions.
Added:
trunk/yt/lagos/TwoPointFunctions.py
- copied, changed from r1759, /trunk/yt/lagos/StructureFunctionGenerator.py
Removed:
trunk/yt/lagos/StructureFunctionGenerator.py
Modified:
trunk/yt/extensions/StarAnalysis.py
trunk/yt/extensions/halo_mass_function.py
trunk/yt/extensions/kdtree/fKD.f90
trunk/yt/extensions/kdtree/fKD.v
trunk/yt/extensions/kdtree/fKD_source.f90
trunk/yt/extensions/merger_tree.py
trunk/yt/lagos/HDF5LightReader.c
trunk/yt/lagos/HaloFinding.py
trunk/yt/lagos/ParticleIO.py
trunk/yt/lagos/__init__.py
trunk/yt/lagos/parallelHOP/parallelHOP.py
trunk/yt/mods.py
Modified: trunk/yt/extensions/StarAnalysis.py
==============================================================================
--- trunk/yt/extensions/StarAnalysis.py (original)
+++ trunk/yt/extensions/StarAnalysis.py Thu Jun 24 17:12:57 2010
@@ -127,7 +127,7 @@
# Use the center of the time_bin, not the left edge.
fp.write("#time\tlookback\tredshift\tMsol/yr\tMsol/yr/Mpc3\tMsol\tcumMsol\t\n")
for i, time in enumerate((self.time_bins[1:] + self.time_bins[:-1])/2.):
- line = "%1.5e\t%1.5e\t%1.5e\t%1.5e\t%1.5e\t%1.5e\t%1.5e\n" % \
+ line = "%1.5e %1.5e %1.5e %1.5e %1.5e %1.5e %1.5e\n" % \
(time * tc / YEAR, # Time
(self.time_now - time * tc)/YEAR, # Lookback time
self.cosm.ComputeRedshiftFromTime(time * tc), # Redshift
@@ -320,7 +320,7 @@
# Add this flux to the total, weighted by mass.
self.final_spec += na.power(10., int_flux) * star[4]
# Normalize.
- self.total_mass = sum(self.star_mass)
+ self.total_mass = na.sum(self.star_mass)
self.avg_mass = na.mean(self.star_mass)
tot_metal = sum(self.star_metal * self.star_mass)
self.avg_metal = math.log10(tot_metal / self.total_mass / Zsun)
Modified: trunk/yt/extensions/halo_mass_function.py
==============================================================================
--- trunk/yt/extensions/halo_mass_function.py (original)
+++ trunk/yt/extensions/halo_mass_function.py Thu Jun 24 17:12:57 2010
@@ -70,6 +70,7 @@
the calculations and generated fit. Default=360.
:param fitting_function (int): Which fitting function to use.
1 = Press-schechter, 2 = Jenkins, 3 = Sheth-Tormen, 4 = Warren fit
+ 5 = Tinker
Default=4.
:param mass_column (int): The column of halo_file that contains the
masses of the haloes. Default=4.
@@ -419,6 +420,15 @@
b=0.2538;
c=1.1982;
thismult = A*( math.pow(sigma, -1.0*a) + b)*math.exp(-1.0*c / sigma / sigma );
+
+ elif self.fitting_function==5:
+ # Tinker et al. 2008, eqn 3, \Delta=300 # \Delta=200
+ A = 0.2 #0.186
+ a = 1.52 #1.47
+ b = 2.25 #2.57
+ c = 1.27 #1.19
+ thismult = A * ( math.pow((sigma / b), -a) + 1) * \
+ math.exp(-1 * c / sigma / sigma)
else:
mylog.error("Don't understand this. Fitting function requested is %d\n",
Modified: trunk/yt/extensions/kdtree/fKD.f90
==============================================================================
--- trunk/yt/extensions/kdtree/fKD.f90 (original)
+++ trunk/yt/extensions/kdtree/fKD.f90 Thu Jun 24 17:12:57 2010
@@ -1,21 +1,59 @@
-subroutine create_tree()
+subroutine create_tree(treeID)
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
+ integer, intent (in), optional :: treeID
+
+ integer :: ID
+
+ if (present(treeID)) then
+ ID = treeID
+ else
+ ID = -1
+ end if
+
! create a kd tree object
+
+ if (ID == 1) then
+ t1%tree2 => kdtree2_create(t1%pos,sort=t1%sort,rearrange=t1%rearrange)
+ elseif (ID == 2) then
+ t2%tree2 => kdtree2_create(t2%pos,sort=t2%sort,rearrange=t2%rearrange)
+ else
+ tree2 => kdtree2_create(pos,sort=sort,rearrange=rearrange)
+ end if
- tree2 => kdtree2_create(pos,sort=sort,rearrange=rearrange) ! this is how you create a tree.
- return
+ return
end subroutine create_tree
+subroutine add_tree(treeID)
+ use kdtree2_module
+ use fKD_module
+ use tree_setmodule
+ use kdtree2module
+ use tree_nodemodule
+ use intervalmodule
+
+ integer :: treeID
+
+ if (treeID == 1) then
+ t1 => Newtree_set()
+ elseif (treeID == 2) then
+ t2 => Newtree_set()
+ end if
+
+ return
+
+end subroutine add_tree
subroutine find_nn_nearest_neighbors()
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
@@ -48,13 +86,14 @@
! nearest neighbors.
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
integer :: k, number
type(kdtree2_result),allocatable :: results(:)
-
+
allocate(results(nn))
number = size(qv_many,2)
@@ -70,15 +109,114 @@
end subroutine find_many_nn_nearest_neighbors
+subroutine find_r_nearest()
+ ! Given an input array of a point (qv), find the number, fortran indices and
+ ! squared distances to all neighbors within (radius) of (qv).
+ ! The number of neighbors with indices and radius returned is set by
+ ! (radius_n), and if this is smaller than (nfound) when returned,
+ ! there are neighbors missing from the results.
+ ! Store the results in (dist) and (tags).
+ use kdtree2_module
+ use fKD_module
+ use tree_setmodule
+ use kdtree2module
+ use tree_nodemodule
+ use intervalmodule
+
+ integer :: k
+ type(kdtree2_result),allocatable :: results(:) ! nearest neighbors
+
+ allocate(results(radius_n))
+ nfound = 0
+
+ ! do it.
+ call kdtree2_r_nearest(tp=tree2,qv=qv,r2=radius,nfound=nfound,nalloc=radius_n,results=results)
+
+ tags(:) = results%idx
+ dist(:) = results%dis
+
+ deallocate(results)
+ return
+
+end subroutine find_r_nearest
+
+subroutine find_many_r_nearest(treeID)
+ ! Given an input array of a points (qv_many), find the number, fortran indices and
+ ! squared distances to all neighbors within (radius) of (qv).
+ ! The number of neighbors with indices and radius returned is set by
+ ! (radius_n), and if this is smaller than (nfound) when returned,
+ ! there are neighbors missing from the results.
+ ! Store the results in (dist) and (tags).
+ use kdtree2_module
+ use fKD_module
+ use tree_setmodule
+ use kdtree2module
+ use tree_nodemodule
+ use intervalmodule
+
+ integer, intent (in), optional :: treeID
+
+ integer :: ID
+ integer :: k, number
+ type(kdtree2_result),allocatable :: results(:) ! nearest neighbors
+
+
+ if (present(treeID)) then
+ ID = treeID
+ else
+ ID = -1
+ end if
+
+ if (ID==1) then
+ allocate(results(t1%radius_n))
+ number = size(t1%qv_many,2)
+ do k=1, number
+ t1%qv(:) = t1%qv_many(:,k)
+ t1%nfound = t1%nfound_many(k)
+ call kdtree2_r_nearest(tp=t1%tree2,qv=t1%qv,r2=t1%radius,nfound=t1%nfound,nalloc=t1%radius_n,results=results)
+ t1%nfound_many(k) = t1%nfound
+ t1%nn_tags(:, k) = results%idx
+ t1%nn_dist(:, k) = results%dis
+ end do
+ elseif (ID==2) then
+ allocate(results(t2%radius_n))
+ number = size(t2%qv_many,2)
+ do k=1, number
+ t2%qv(:) = t2%qv_many(:,k)
+ t2%nfound = t2%nfound_many(k)
+ call kdtree2_r_nearest(tp=t2%tree2,qv=t2%qv,r2=t2%radius,nfound=t2%nfound,nalloc=t2%radius_n,results=results)
+ t2%nfound_many(k) = t2%nfound
+ t2%nn_tags(:, k) = results%idx
+ t2%nn_dist(:, k) = results%dis
+ end do
+ else
+ allocate(results(radius_n))
+ number = size(qv_many,2)
+ do k=1, number
+ qv(:) = qv_many(:,k)
+ nfound = nfound_many(k)
+ call kdtree2_r_nearest(tp=tree2,qv=qv,r2=radius,nfound=nfound,nalloc=radius_n,results=results)
+ nfound_many(k) = nfound
+ nn_tags(:, k) = results%idx
+ nn_dist(:, k) = results%dis
+ end do
+ end if
+
+ deallocate(results)
+ return
+
+end subroutine find_many_r_nearest
+
subroutine find_all_nn_nearest_neighbors()
! for all particles in pos, find their nearest neighbors and return the
! indexes and distances as big arrays
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
-
+
integer :: k
type(kdtree2_result),allocatable :: results(:) ! nearest neighbors
allocate(results(nn))
@@ -99,6 +237,7 @@
! for a chunk of the full number of particles, find their nearest neighbors
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
@@ -123,6 +262,7 @@
! their density. Return only nMerge nearest neighbors.
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
@@ -172,16 +312,33 @@
end subroutine chainHOP_tags_dens
-subroutine free_tree()
+subroutine free_tree(treeID)
use kdtree2_module
use fKD_module
+ use tree_setmodule
use kdtree2module
use tree_nodemodule
use intervalmodule
+
+ integer, intent (in), optional :: treeID
+
+ integer :: ID
+
+ if (present(treeID)) then
+ ID = treeID
+ else
+ ID = -1
+ end if
! this releases memory for the tree BUT NOT THE ARRAY OF DATA YOU PASSED
! TO MAKE THE TREE.
- call kdtree2_destroy(tree2)
+ if (ID == 1) then
+ call kdtree2_destroy(t1%tree2)
+ elseif (ID == 2) then
+ call kdtree2_destroy(t2%tree2)
+ else
+ call kdtree2_destroy(tree2)
+ end if
! The data to make the tree has to be deleted in python BEFORE calling
! this!
Modified: trunk/yt/extensions/kdtree/fKD.v
==============================================================================
--- trunk/yt/extensions/kdtree/fKD.v (original)
+++ trunk/yt/extensions/kdtree/fKD.v Thu Jun 24 17:12:57 2010
@@ -1,7 +1,10 @@
fKD
****** fKD_module vars:
-# Not all of these are being used, but they only take memory
+t1 _tree_set
+t2 _tree_set
+# This contains the objects that any given kD tree might wish to use.
+# Not all of these are being used all the time, but they only take memory
# if they're initialized in python.
tags(:) _integer # particle ID tags
dist(:) _real # interparticle spacings
@@ -21,6 +24,10 @@
tree2 _kdtree2
sort logical /.false./
rearrange logical /.true./
+radius real # the unsquared radius for radius searches
+radius_n integer # the number of results to return
+nfound integer # number of neighbors within radius
+nfound_many(:) _integer # an array of number of neighbors within radius
%%%% interval:
lower real
@@ -92,12 +99,41 @@
# root pointer of the tree
+%%%% tree_set:
+# This contains the objects that any given kD tree might wish to use.
+# Not all of these are being used all the time, but they only take memory
+# if they're initialized in python.
+tags(:) _integer # particle ID tags
+dist(:) _real # interparticle spacings
+nn_tags(:,:) _integer # for all particles at once, [nth neighbor, index]
+chunk_tags(:,:) _integer # for finding only a chunk of the nearest neighbors
+nn_dist(:,:) _real
+pos(3,:) _real
+dens(:) _real
+mass(:) _real
+qv(3) real
+qv_many(3,:) _real
+nparts integer
+nn integer
+nMerge integer # number of nearest neighbors used in chain merging
+start integer
+finish integer
+tree2 _kdtree2
+sort logical /.false./
+rearrange logical /.true./
+radius real # the unsquared radius for radius searches
+radius_n integer # the number of results to return
+nfound integer # number of neighbors within radius
+nfound_many(:) _integer # an array of number of neighbors within radius
***** Subroutines:
find_nn_nearest_neighbors subroutine
-create_tree() subroutine
-free_tree() subroutine
+create_tree(treeID:integer) subroutine
+add_tree(treeID:integer) subroutine
+free_tree(treeID:integer) subroutine
find_all_nn_nearest_neighbors subroutine
+find_r_nearest subroutine
+find_many_r_nearest(treeID:integer) subroutine
find_many_nn_nearest_neighbors subroutine
find_chunk_nearest_neighbors subroutine
chainHOP_tags_dens subroutine
Modified: trunk/yt/extensions/kdtree/fKD_source.f90
==============================================================================
--- trunk/yt/extensions/kdtree/fKD_source.f90 (original)
+++ trunk/yt/extensions/kdtree/fKD_source.f90 Thu Jun 24 17:12:57 2010
@@ -1166,11 +1166,11 @@
call kdtree2_sort_results(nfound, results)
endif
- if (sr%overflow) then
- write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors'
- write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball'
- write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.'
- endif
+ !if (sr%overflow) then
+ ! write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors'
+ ! write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball'
+ ! write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.'
+ !endif
return
end subroutine kdtree2_r_nearest
Modified: trunk/yt/extensions/merger_tree.py
==============================================================================
--- trunk/yt/extensions/merger_tree.py (original)
+++ trunk/yt/extensions/merger_tree.py Thu Jun 24 17:12:57 2010
@@ -24,10 +24,13 @@
"""
import yt.lagos as lagos
-from yt.extensions.kdtree import *
from yt.lagos.HaloFinding import HaloFinder
from yt.logger import lagosLogger as mylog
import yt.extensions.HaloProfiler as HP
+try:
+ from yt.extensions.kdtree import *
+except ImportError:
+ mylog.debug("The Fortran kD-Tree did not import correctly.")
import numpy as na
import os, glob, md5, time, gc, sys
@@ -303,7 +306,7 @@
fKD.nn = 5
fKD.sort = True
fKD.rearrange = True
- create_tree()
+ create_tree(0)
# Find the parent points from the database.
parent_pf = lagos.EnzoStaticOutput(parentfile)
@@ -330,7 +333,7 @@
candidates[row[0]] = nIDs
del fKD.pos, fKD.tags, fKD.dist
- free_tree() # Frees the kdtree object.
+ free_tree(0) # Frees the kdtree object.
self.candidates = candidates
Modified: trunk/yt/lagos/HDF5LightReader.c
==============================================================================
--- trunk/yt/lagos/HDF5LightReader.c (original)
+++ trunk/yt/lagos/HDF5LightReader.c Thu Jun 24 17:12:57 2010
@@ -80,9 +80,18 @@
npy_float64 radius;
} sphere_validation;
+typedef struct cylinder_validation_ {
+ /* These cannot contain any pointers */
+ npy_float64 center[3];
+ npy_float64 normal[3];
+ npy_float64 radius;
+ npy_float64 height;
+} cylinder_validation;
+
/* Forward declarations */
int setup_validator_region(particle_validation *data, PyObject *InputData);
int setup_validator_sphere(particle_validation *data, PyObject *InputData);
+int setup_validator_cylinder(particle_validation *data, PyObject *InputData);
int run_validators(particle_validation *pv, char *filename,
int grid_id, const int read, const int packed,
int grid_index);
@@ -97,6 +106,10 @@
int count_particles_sphere_DOUBLE(particle_validation *data);
int count_particles_sphere_LONGDOUBLE(particle_validation *data);
+int count_particles_cylinder_FLOAT(particle_validation *data);
+int count_particles_cylinder_DOUBLE(particle_validation *data);
+int count_particles_cylinder_LONGDOUBLE(particle_validation *data);
+
int get_my_desc_type(hid_t native_type_id){
/*
@@ -820,6 +833,10 @@
/* Sphere type */
setup_validator_sphere(&pv, vargs);
break;
+ case 2:
+ /* Cylinder type */
+ setup_validator_cylinder(&pv, vargs);
+ break;
default:
PyErr_Format(_hdf5ReadError,
"Unrecognized data source.\n");
@@ -995,6 +1012,40 @@
return 1;
}
+int setup_validator_cylinder(particle_validation *data, PyObject *InputData)
+{
+ int i;
+ /* These are borrowed references */
+ PyArrayObject *center = (PyArrayObject *) PyTuple_GetItem(InputData, 0);
+ PyArrayObject *normal = (PyArrayObject *) PyTuple_GetItem(InputData, 1);
+ PyObject *radius = (PyObject *) PyTuple_GetItem(InputData, 2);
+ PyObject *height = (PyObject *) PyTuple_GetItem(InputData, 3);
+
+ /* This will get freed in the finalization of particle validation */
+ cylinder_validation *cv = (cylinder_validation *)
+ malloc(sizeof(cylinder_validation));
+ data->validation_reqs = (void *) cv;
+
+ for (i = 0; i < 3; i++){
+ cv->center[i] = *(npy_float64*) PyArray_GETPTR1(center, i);
+ }
+
+ for (i = 0; i < 3; i++){
+ cv->normal[i] = *(npy_float64*) PyArray_GETPTR1(normal, i);
+ }
+
+ cv->radius = (npy_float64) PyFloat_AsDouble(radius);
+ cv->height = (npy_float64) PyFloat_AsDouble(height);
+
+ data->count_func = NULL;
+ data->count_func_float = count_particles_cylinder_FLOAT;
+ data->count_func_double = count_particles_cylinder_DOUBLE;
+ data->count_func_longdouble = count_particles_cylinder_LONGDOUBLE;
+
+ return 1;
+}
+
+
int run_validators(particle_validation *pv, char *filename,
int grid_id, const int read, const int packed,
int grid_index)
@@ -1563,6 +1614,149 @@
return n;
}
+int count_particles_cylinder_FLOAT(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ cylinder_validation *vdata;
+
+ vdata = (cylinder_validation*) data->validation_reqs;
+
+ float **particle_data = (float **) data->particle_position;
+
+ float *particle_position_x = particle_data[0];
+ float *particle_position_y = particle_data[1];
+ float *particle_position_z = particle_data[2];
+ float temph, tempd, d;
+
+ d = -1. * (vdata->normal[0] * vdata->center[0] +
+ vdata->normal[1] * vdata->center[1] +
+ vdata->normal[2] * vdata->center[2]);
+
+ double pradius, ph, pd;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0; ph = 0.0; pd = 0.0;
+
+ temph = (particle_position_x[ind] * vdata->normal[0]); ph += temph;
+ temph = (particle_position_y[ind] * vdata->normal[1]); ph += temph;
+ temph = (particle_position_z[ind] * vdata->normal[2]); ph += temph;
+ ph += d;
+
+ tempd = (particle_position_x[ind] - vdata->center[0]); pd += tempd*tempd;
+ tempd = (particle_position_y[ind] - vdata->center[1]); pd += tempd*tempd;
+ tempd = (particle_position_z[ind] - vdata->center[2]); pd += tempd*tempd;
+
+ pradius = pow(pd - ph*ph, 0.5);
+ if ((pradius <= vdata->radius) && (fabs(ph) <= vdata->height)) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
+
+int count_particles_cylinder_DOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ cylinder_validation *vdata;
+
+ vdata = (cylinder_validation*) data->validation_reqs;
+
+ double **particle_data = (double **) data->particle_position;
+
+ double *particle_position_x = particle_data[0];
+ double *particle_position_y = particle_data[1];
+ double *particle_position_z = particle_data[2];
+ double temph, tempd, d;
+
+ d = -1. * (vdata->normal[0] * vdata->center[0] +
+ vdata->normal[1] * vdata->center[1] +
+ vdata->normal[2] * vdata->center[2]);
+
+ double pradius, ph, pd;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0; ph = 0.0; pd = 0.0;
+
+ temph = (particle_position_x[ind] * vdata->normal[0]); ph += temph;
+ temph = (particle_position_y[ind] * vdata->normal[1]); ph += temph;
+ temph = (particle_position_z[ind] * vdata->normal[2]); ph += temph;
+ ph += d;
+
+ tempd = (particle_position_x[ind] - vdata->center[0]); pd += tempd*tempd;
+ tempd = (particle_position_y[ind] - vdata->center[1]); pd += tempd*tempd;
+ tempd = (particle_position_z[ind] - vdata->center[2]); pd += tempd*tempd;
+
+ pradius = pow(pd - ph*ph, 0.5);
+ if ((pradius <= vdata->radius) && (fabs(ph) <= vdata->height)) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
+
+int count_particles_cylinder_LONGDOUBLE(particle_validation *data)
+{
+ /* Our data comes packed in a struct, off which our pointers all hang */
+
+ /* First is our validation requirements, which are a set of three items: */
+
+ int ind, n=0;
+ cylinder_validation *vdata;
+
+ vdata = (cylinder_validation*) data->validation_reqs;
+
+ long double **particle_data = (long double **) data->particle_position;
+
+ long double *particle_position_x = particle_data[0];
+ long double *particle_position_y = particle_data[1];
+ long double *particle_position_z = particle_data[2];
+ long double temph, tempd, d;
+
+ d = -1. * (vdata->normal[0] * vdata->center[0] +
+ vdata->normal[1] * vdata->center[1] +
+ vdata->normal[2] * vdata->center[2]);
+
+ long double pradius, ph, pd;
+
+ for (ind = 0; ind < data->particles_to_check; ind++) {
+ pradius = 0.0; ph = 0.0; pd = 0.0;
+
+ temph = (particle_position_x[ind] * vdata->normal[0]); ph += temph;
+ temph = (particle_position_y[ind] * vdata->normal[1]); ph += temph;
+ temph = (particle_position_z[ind] * vdata->normal[2]); ph += temph;
+ ph += d;
+
+ tempd = (particle_position_x[ind] - vdata->center[0]); pd += tempd*tempd;
+ tempd = (particle_position_y[ind] - vdata->center[1]); pd += tempd*tempd;
+ tempd = (particle_position_z[ind] - vdata->center[2]); pd += tempd*tempd;
+
+ pradius = pow(pd - ph*ph, 0.5);
+ if ((pradius <= vdata->radius) && (fabsl(ph) <= vdata->height)) {
+ if(data->update_count == 1) data->total_valid_particles++;
+ data->mask[ind] = 1;
+ n++;
+ } else {
+ data->mask[ind] = 0;
+ }
+ }
+ return n;
+}
static PyMethodDef _hdf5LightReaderMethods[] = {
{"ReadData", Py_ReadHDF5DataSet, METH_VARARGS},
Modified: trunk/yt/lagos/HaloFinding.py
==============================================================================
--- trunk/yt/lagos/HaloFinding.py (original)
+++ trunk/yt/lagos/HaloFinding.py Thu Jun 24 17:12:57 2010
@@ -1141,12 +1141,14 @@
class parallelHF(GenericHaloFinder, parallelHOPHaloList):
def __init__(self, pf, threshold=160, dm_only=True, resize=True, rearrange=True,\
- fancy_padding=True, safety=1.5, premerge=True):
+ fancy_padding=True, safety=1.5, premerge=True, sample=0.03):
GenericHaloFinder.__init__(self, pf, dm_only, padding=0.0)
- self.padding = 0.0
+ self.padding = 0.0
self.num_neighbors = 65
self.safety = safety
+ self.sample = sample
period = pf["DomainRightEdge"] - pf["DomainLeftEdge"]
+ topbounds = na.array([[0., 0., 0.], period])
# Cut up the volume evenly initially, with no padding.
padded, LE, RE, self._data_source = self._partition_hierarchy_3d(padding=self.padding)
# also get the total mass of particles
@@ -1154,40 +1156,17 @@
# Adaptive subregions by bisection.
ds_names = ["particle_position_x","particle_position_y","particle_position_z"]
if ytcfg.getboolean("yt","inline") == False and \
- resize and self._mpi_get_size() != 1:
+ resize and self._mpi_get_size() != 1:
+ random.seed(self._mpi_get_rank())
cut_list = self._partition_hierarchy_3d_bisection_list()
- for i,cut in enumerate(cut_list):
- dim = cut[0]
- if i == 0:
- width = pf["DomainRightEdge"][dim] - pf["DomainLeftEdge"][dim]
- new_LE = pf["DomainLeftEdge"]
- else:
- new_LE, new_RE = new_top_bounds
- width = new_RE[dim] - new_LE[dim]
- if ytcfg.getboolean("yt","inline") == False:
- data = self._data_source[ds_names[dim]]
- else:
- data = self._data_source[ds_names[dim]]
- if i == 0:
- local_parts = data.size
- n_parts = self._mpi_allsum(local_parts)
- min = self._mpi_allmin(local_parts)
- max = self._mpi_allmax(local_parts)
- mylog.info("Initial distribution: (min,max,tot) (%d,%d,%d)" \
- % (min, max, n_parts))
- num_bins = 1000
- bin_width = float(width)/float(num_bins)
- bins = na.arange(num_bins+1, dtype='float64') * bin_width + new_LE[dim]
- counts, bins = na.histogram(data, bins, new=True)
- if i == 0:
- new_group, new_comm, LE, RE, new_top_bounds, new_cc, self._data_source= \
- self._partition_hierarchy_3d_bisection(dim, bins, counts, top_bounds = None,\
- old_group = None, old_comm = None, cut=cut)
- else:
- new_group, new_comm, LE, RE, new_top_bounds, new_cc, self._data_source = \
- self._partition_hierarchy_3d_bisection(dim, bins, counts, top_bounds = new_top_bounds,\
- old_group = new_group, old_comm = new_comm, cut=cut, old_cc=new_cc)
- del bins, counts
+ root_points = self._subsample_points()
+ self.bucket_bounds = []
+ if self._mpi_get_rank() == 0:
+ self._recursive_divide(root_points, topbounds, 0, cut_list)
+ self.bucket_bounds = self._mpi_bcast_pickled(self.bucket_bounds)
+ my_bounds = self.bucket_bounds[self._mpi_get_rank()]
+ LE, RE = my_bounds[0], my_bounds[1]
+ self._data_source = self.hierarchy.region_strict([0.]*3, LE, RE)
# If this isn't parallel, define the region as an AMRRegionStrict so
# particle IO works.
if self._mpi_get_size() == 1:
@@ -1215,6 +1194,8 @@
# Another approach to padding, perhaps more accurate.
elif fancy_padding and self._distributed:
LE_padding, RE_padding = na.empty(3,dtype='float64'), na.empty(3,dtype='float64')
+ avg_spacing = (float(vol) / data.size)**(1./3.)
+ base_padding = (self.num_neighbors)**(1./3.) * self.safety * avg_spacing
for dim in xrange(3):
if ytcfg.getboolean("yt","inline") == False:
data = self._data_source[ds_names[dim]]
@@ -1224,7 +1205,8 @@
width = self._data_source.right_edge[dim] - self._data_source.left_edge[dim]
area = (self._data_source.right_edge[(dim+1)%3] - self._data_source.left_edge[(dim+1)%3]) * \
(self._data_source.right_edge[(dim+2)%3] - self._data_source.left_edge[(dim+2)%3])
- bin_width = float(width)/float(num_bins)
+ bin_width = base_padding
+ num_bins = int(math.ceil(width / bin_width))
bins = na.arange(num_bins+1, dtype='float64') * bin_width + self._data_source.left_edge[dim]
counts, bins = na.histogram(data, bins, new=True)
# left side.
@@ -1266,6 +1248,89 @@
self._join_halolists()
yt_counters("Final Grouping")
+ def _subsample_points(self):
+ # Read in a random subset of the points in each domain, and then
+ # collect them on the root task.
+ xp = self._data_source["particle_position_x"]
+ n_parts = self._mpi_allsum(xp.size)
+ local_parts = xp.size
+ random_points = int(self.sample * n_parts)
+ # We want to get a representative selection of random particles in
+ # each subvolume.
+ adjust = float(local_parts) / ( float(n_parts) / self._mpi_get_size())
+ n_random = int(adjust * float(random_points) / self._mpi_get_size())
+ mylog.info("Reading in %d random particles." % n_random)
+ # Get unique random particles.
+ my_points = na.empty((n_random, 3), dtype='float64')
+ uni = na.array(random.sample(xrange(xp.size), n_random))
+ uni = uni[uni.argsort()]
+ my_points[:,0] = xp[uni]
+ del xp
+ self._data_source.clear_data()
+ my_points[:,1] = self._data_source["particle_position_y"][uni]
+ self._data_source.clear_data()
+ my_points[:,2] = self._data_source["particle_position_z"][uni]
+ self._data_source.clear_data()
+ del uni
+ # Collect them on the root task.
+ mine, sizes = self._mpi_info_dict(n_random)
+ if mine == 0:
+ tot_random = sum(sizes.values())
+ root_points = na.empty((tot_random, 3), dtype='float64')
+ root_points.shape = (1, 3*tot_random)
+ else:
+ root_points = na.empty([])
+ my_points.shape = (1, n_random*3)
+ root_points = self._mpi_concatenate_array_on_root_double(my_points[0])
+ del my_points
+ if mine == 0:
+ root_points.shape = (tot_random, 3)
+ return root_points
+
+ def _recursive_divide(self, points, bounds, level, cut_list):
+ dim = cut_list[level][0]
+ parts = points.shape[0]
+ num_bins = 1000
+ width = bounds[1][dim] - bounds[0][dim]
+ bin_width = width / num_bins
+ bins = na.arange(num_bins+1, dtype='float64') * bin_width + bounds[0][dim]
+ counts, bins = na.histogram(points[:,dim], bins)
+ # Find the bin that passes the cut points.
+ midpoints = [bounds[0][dim]]
+ sum = 0
+ bin = 0
+ for step in xrange(1,cut_list[level][1]):
+ while sum < ((parts*step)/cut_list[level][1]):
+ lastsum = sum
+ sum += counts[bin]
+ bin += 1
+ # Bin edges
+ left_edge = bins[bin-1]
+ right_edge = bins[bin]
+ # Find a better approx of the midpoint cut line using a linear approx.
+ a = float(sum - lastsum) / (right_edge - left_edge)
+ midpoints.append(left_edge + (0.5 - (float(lastsum) / parts / 2)) / a)
+ midpoints.append(bounds[1][dim])
+
+ # Split the points & update the bounds.
+ subpoints = []
+ subbounds = []
+ for pair in zip(midpoints[:-1],midpoints[1:]):
+ select = na.bitwise_and(points[:,dim] >= pair[0],
+ points[:,dim] < pair[1])
+ subpoints.append(points[select])
+ nb = bounds.copy()
+ nb[0][dim] = pair[0]
+ nb[1][dim] = pair[1]
+ subbounds.append(nb)
+ # If we're at the maxlevel, make a bucket. Otherwise, recurse down.
+ maxlevel = len(cut_list) - 1
+ for pair in zip(subpoints, subbounds):
+ if level == maxlevel:
+ self.bucket_bounds.append(pair[1])
+ else:
+ self._recursive_divide(pair[0], pair[1], level+1, cut_list)
+
def _join_halolists(self):
if self.group_count == 0:
mylog.info("There are no halos found.")
Modified: trunk/yt/lagos/ParticleIO.py
==============================================================================
--- trunk/yt/lagos/ParticleIO.py (original)
+++ trunk/yt/lagos/ParticleIO.py Thu Jun 24 17:12:57 2010
@@ -128,3 +128,20 @@
def _get_args(self):
return (1, (na.array(self.center, dtype='float64'), self.radius))
+
+class ParticleIOHandlerDisk(ParticleIOHandlerImplemented):
+ _source_type = "disk"
+
+ def __init__(self, pf, source):
+ self.center = source.center
+ self.normal = source._norm_vec
+ self.radius = source._radius
+ self.height = source._height
+ ParticleIOHandler.__init__(self, pf, source)
+
+ def _get_args(self):
+ args = (na.array(self.center, dtype='float64'),
+ na.array(self.normal, dtype='float64'),
+ self.radius, self.height)
+ return (2, args)
+
\ No newline at end of file
Copied: trunk/yt/lagos/TwoPointFunctions.py (from r1759, /trunk/yt/lagos/StructureFunctionGenerator.py)
==============================================================================
--- /trunk/yt/lagos/StructureFunctionGenerator.py (original)
+++ trunk/yt/lagos/TwoPointFunctions.py Thu Jun 24 17:12:57 2010
@@ -1,5 +1,5 @@
"""
-Structure Function Generator
+Two Point Functions Framework.
Author: Stephen Skory <stephenskory at yahoo.com>
Affiliation: UCSD Physics/CASS
@@ -29,19 +29,21 @@
try:
from yt.extensions.kdtree import *
except ImportError:
- mylog.debug("The Fortran kD-Tree did not import correctly. The structure function generator will not work correctly.")
+ mylog.debug("The Fortran kD-Tree did not import correctly.")
import math, sys, itertools, inspect, types, time
from collections import defaultdict
-class StructFcnGen(ParallelAnalysisInterface):
+sep = 12
+
+class TwoPointFunctions(ParallelAnalysisInterface):
def __init__(self, pf, fields, left_edge=None, right_edge=None,
total_values=1000000, comm_size=10000, length_type="lin",
length_number=10, length_range=None, vol_ratio = 1,
salt=0):
"""
*total_values* (int) How many total (global) pair calculations to run
- for each of the structure functions specified.
+ for each of the functions specified.
A single integer. Default: 1,000,000.
*comm_size* (int) How entries are sent during communication.
Default: 10,000.
@@ -57,10 +59,21 @@
the simulational volume.
A single pair (a list or array). Default: [sqrt(3)dx, 1/2*shortest box edge],
where dx is the smallest grid cell size.
+ *vol_ratio* (int) How to multiply-assign subvolumes to the parallel
+ tasks. This number must be an integer factor of the total number of tasks or
+ very bad things will happen. The default value of 1 will assign one task
+ to each subvolume, and there will be an equal number of subvolumes as tasks.
+ A value of 2 will assign two tasks to each subvolume and there will be
+ one-half as many subvolumes as tasks.
+ A value equal to the number of parallel tasks will result in each task
+ owning a complete copy of all the fields data, meaning each task will be
+ operating on the identical full volume.
+ Setting it to -1 will automatically adjust it such that each task
+ owns the entire volume.
*salt* (int) a number that will be added to the random number generator
seed. Use this if a different random series of numbers is desired when
keeping everything else constant from this set: (MPI task count,
- number of ruler lengths, ruler min/max, number of structure functions,
+ number of ruler lengths, ruler min/max, number of functions,
number of point pairs per ruler length). Default: 0.
"""
try:
@@ -73,6 +86,8 @@
self.size = self._mpi_get_size()
self.mine = self._mpi_get_rank()
self.vol_ratio = vol_ratio
+ if self.vol_ratio == -1:
+ self.vol_ratio = self.size
self.total_values = int(total_values / self.size)
# For communication.
self.recv_hooks = []
@@ -132,14 +147,14 @@
self.width = self.ds.right_edge - self.ds.left_edge
self.mt = na.random.mtrand.RandomState(seed = 1234 * self.mine + salt)
- def add_function(self, function, out_labels, sqrt):
+ def add_function(self, function, out_labels, sqrt, corr_norm=None):
"""
- Add a structure function to the list that will be evaluated at the
+ Add a function to the list that will be evaluated at the
generated pairs of points.
"""
fargs = inspect.getargspec(function)
if len(fargs.args) != 5:
- raise SyntaxError("The structure function %s needs five arguments." %\
+ raise SyntaxError("The function %s needs five arguments." %\
function.__name__)
out_labels = list(out_labels)
if len(out_labels) < 1:
@@ -149,8 +164,8 @@
if len(sqrt) != len(out_labels):
raise SyntaxError("Please have the same number of elements in out_labels as in sqrt for function %s." %\
function.__name__)
- self._fsets.append(StructSet(self, function, self.min_edge,
- out_labels, sqrt))
+ self._fsets.append(FcnSet(self, function, self.min_edge,
+ out_labels, sqrt,corr_norm))
return self._fsets[-1]
def __getitem__(self, key):
@@ -158,15 +173,15 @@
def run_generator(self):
"""
- After all the structure functions have been added, run the generator.
+ After all the functions have been added, run the generator.
"""
yt_counters("run_generator")
# We need a function!
if len(self._fsets) == 0:
- mylog.error("You need to add at least one structure function!")
+ mylog.error("You need to add at least one function!")
return None
# Do all the startup tasks to get the grid points.
- if self.nlevels == 1:
+ if self.nlevels == 0:
yt_counters("build_sort")
self._build_sort_array()
self.sort_done = False
@@ -227,9 +242,9 @@
mylog.info("Length (%d of %d) %1.5e took %d communication cycles to complete." % \
(bigloop+1, len(self.lengths), length, self.comm_cycle_count))
yt_counters("big loop over lengths")
- if self.nlevels > 1:
+ if self.nlevels >= 1:
del fKD.pos, fKD.qv_many, fKD.nn_tags
- free_tree() # Frees the kdtree object.
+ free_tree(0) # Frees the kdtree object.
yt_counters("allsum")
self._allsum_bin_hits()
mylog.info("Spent %f seconds waiting for communication." % t_waiting)
@@ -252,7 +267,7 @@
fKD.nn = 1
fKD.sort = False
fKD.rearrange = True
- create_tree()
+ create_tree(0)
def _build_sort_array(self):
"""
@@ -422,7 +437,7 @@
"""
Finds the closest grid cell for each point in a vectorized manner.
"""
- if self.nlevels == 1:
+ if self.nlevels == 0:
pos = (points - self.ds.left_edge) / self.width
n = (self.sizes[2] * pos[:,2]).astype('int32')
n += self.sizes[2] * (self.sizes[1] * pos[:,1]).astype('int32')
@@ -501,7 +516,7 @@
self.fields_vals.shape = (self.comm_size, len(self.fields)*2)
self.points.shape = (self.comm_size, 6)
- # To run the structure functions, what is key is that the
+ # To run the functions, what is key is that the
# second point in the pair is ours.
second_points = self.points[:,3:]
select = na.bitwise_or((second_points < self.ds.left_edge).any(axis=1),
@@ -535,12 +550,11 @@
@parallel_blocking_call
def write_out_means(self):
"""
- Writes out the weighted-average value for each structure function for
+ Writes out the weighted-average value for each function for
each dimension for each ruler length to a text file. The data is written
to files of the name 'function_name.txt' in the current working
directory.
"""
- sep = 12
for fset in self._fsets:
fp = self._write_on_root("%s.txt" % fset.function.__name__)
fset._avg_bin_hits()
@@ -567,7 +581,7 @@
def write_out_arrays(self):
"""
Writes out the raw probability bins and the bin edges to an HDF5 file
- for each of the structure functions. The files are named
+ for each of the functions. The files are named
'function_name.txt' and saved in the current working directory.
"""
if self.mine == 0:
@@ -594,13 +608,35 @@
f.create_dataset("/counts", data=bin_counts)
f.close()
-class StructSet(StructFcnGen):
- def __init__(self,sfg, function, min_edge, out_labels, sqrt):
- self.sfg = sfg # The overarching SFG class
+ @parallel_root_only
+ def write_out_correlation(self):
+ """
+ A special output function for doing two point correlation functions.
+ Outputs the correlation function xi(r) in a text file
+ 'function_name_corr.txt' in the current working directory.
+ """
+ for fset in self._fsets:
+ # Only operate on correlation functions.
+ if fset.corr_norm == None: continue
+ fp = self._write_on_root("%s_correlation.txt" % fset.function.__name__)
+ line = "# length".ljust(sep)
+ line += "\\xi".ljust(sep)
+ fp.write(line + "\n")
+ xi = fset._corr_sum_norm()
+ for length in self.lengths:
+ line = ("%1.5e" % length).ljust(sep)
+ line += ("%1.5e" % xi[length]).ljust(sep)
+ fp.write(line + "\n")
+ fp.close()
+
+class FcnSet(TwoPointFunctions):
+ def __init__(self,tpf, function, min_edge, out_labels, sqrt, corr_norm):
+ self.tpf = tpf # The overarching TPF class
self.function = function # Function to eval between the two points.
self.min_edge = min_edge # The length of the minimum edge of the box.
self.out_labels = out_labels # For output.
self.sqrt = sqrt # which columns to sqrt on output.
+ self.corr_norm = corr_norm # A number used to normalize a correlation function.
# These below are used to track how many times the function returns
# unbinned results.
self.too_low = na.zeros(len(self.out_labels), dtype='int32')
@@ -622,7 +658,7 @@
Default: None.
"""
# This should be called after setSearchParams.
- if not hasattr(self.sfg, "lengths"):
+ if not hasattr(self.tpf, "lengths"):
mylog.error("Please call setSearchParams() before calling setPDFParams().")
return None
# Make sure they're either all lists or only one is.
@@ -657,7 +693,7 @@
# Create the dict that stores the arrays to store the bin hits, and
# the arrays themselves.
self.length_bin_hits = {}
- for length in self.sfg.lengths:
+ for length in self.tpf.lengths:
# It's easier to index flattened, but will be unflattened later.
self.length_bin_hits[length] = na.zeros(self.bin_number,
dtype='int64').flatten()
@@ -682,10 +718,10 @@
def _eval_st_fcn(self, results, points, vec):
"""
- Return the value of the structure function using the provided results.
+ Return the value of the function using the provided results.
"""
- return self.function(results[:,:len(self.sfg.fields)],
- results[:,len(self.sfg.fields):], points[:,:3], points[:,3:], vec)
+ return self.function(results[:,:len(self.tpf.fields)],
+ results[:,len(self.tpf.fields):], points[:,:3], points[:,3:], vec)
"""
NOTE - A function looks like:
def stuff(a,b,r1,r2, vec):
@@ -742,9 +778,19 @@
For each dimension and length of bin_hits return the weighted average.
"""
self.length_avgs = defaultdict(dict)
- for length in self.sfg.lengths:
+ for length in self.tpf.lengths:
for dim in self.dims:
self.length_avgs[length][dim] = \
(self._dim_sum(self.length_bin_hits[length], dim) * \
((self.bin_edges[dim][:-1] + self.bin_edges[dim][1:]) / 2.)).sum()
+ def _corr_sum_norm(self):
+ """
+ Return the correlations xi for this function. We are tacitly assuming
+ that all correlation functions are one dimensional.
+ """
+ xi = {}
+ for length in self.tpf.lengths:
+ xi[length] = -1 + na.sum(self.length_bin_hits[length] * \
+ self.bin_edges[0][:-1]) / self.corr_norm
+ return xi
\ No newline at end of file
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Thu Jun 24 17:12:57 2010
@@ -100,7 +100,7 @@
from HaloFinding import *
-from StructureFunctionGenerator import *
+from TwoPointFunctions import *
# We load plugins. Keep in mind, this can be fairly dangerous -
# the primary purpose is to allow people to have a set of functions
Modified: trunk/yt/lagos/parallelHOP/parallelHOP.py
==============================================================================
--- trunk/yt/lagos/parallelHOP/parallelHOP.py (original)
+++ trunk/yt/lagos/parallelHOP/parallelHOP.py Thu Jun 24 17:12:57 2010
@@ -28,8 +28,11 @@
from yt.lagos import *
from yt.funcs import *
-from yt.extensions.kdtree import *
from yt.performance_counters import yt_counters, time_function
+try:
+ from yt.extensions.kdtree import *
+except ImportError:
+ mylog.debug("The Fortran kD-Tree did not import correctly.")
class RunParallelHOP(ParallelAnalysisInterface):
def __init__(self,period, padding, num_neighbors, bounds,
@@ -348,7 +351,7 @@
fKD.sort = True # Slower, but needed in _connect_chains
fKD.rearrange = self.rearrange # True is faster, but uses more memory
# Now call the fortran.
- create_tree()
+ create_tree(0)
self.__max_memory()
yt_counters("init kd tree")
@@ -1449,7 +1452,7 @@
self._connect_chains()
del fKD.dens, fKD.mass, fKD.dens
del fKD.pos, fKD.chunk_tags
- free_tree() # Frees the kdtree object.
+ free_tree(0) # Frees the kdtree object.
del self.densestNN
mylog.info('Communicating group links globally...')
self._make_global_chain_densest_n()
Modified: trunk/yt/mods.py
==============================================================================
--- trunk/yt/mods.py (original)
+++ trunk/yt/mods.py Thu Jun 24 17:12:57 2010
@@ -47,7 +47,7 @@
Clump, write_clump_hierarchy, find_clumps, write_clumps, \
get_lowest_clumps, \
OrionStaticOutput, HaloFinder, HOPHaloFinder, FOFHaloFinder, parallelHF, \
- axis_names, x_dict, y_dict, StructFcnGen, StructSet
+ axis_names, x_dict, y_dict, TwoPointFunctions, FcnSet
# This is a temporary solution -- in the future, we will allow the user to
# select this via ytcfg.
More information about the yt-svn
mailing list