[Yt-svn] yt-commit r987 - trunk/yt/extensions/lightcone

britton at wrangler.dreamhost.com britton at wrangler.dreamhost.com
Mon Dec 8 21:12:30 PST 2008


Author: britton
Date: Mon Dec  8 21:12:29 2008
New Revision: 987
URL: http://yt.spacepope.org/changeset/987

Log:
Added halo mask generator to cut out circles centered on halos with radius of 
the virial from light cones.  This can be called as a stand-alone or from 
within the light cone generator.  An addition to the sample script will follow.

Added various improvements to the main light cone projection routine, including 
removing redundant calculations of various boolean masks used in the shifting 
procedure.  The depth problem is now handle by performing field_cuts on the 
line of site coordinate field and corresponding dx, instead of extracting a 
periodic region and throwing that as a source to the raven projection routine.  
This will allow the light cones to be run in parallel and has also proven to 
be much faster than the original method.


Added:
   trunk/yt/extensions/lightcone/HaloMask.py
Modified:
   trunk/yt/extensions/lightcone/LightCone.py
   trunk/yt/extensions/lightcone/LightConeProjection.py
   trunk/yt/extensions/lightcone/__init__.py

Added: trunk/yt/extensions/lightcone/HaloMask.py
==============================================================================
--- (empty file)
+++ trunk/yt/extensions/lightcone/HaloMask.py	Mon Dec  8 21:12:29 2008
@@ -0,0 +1,225 @@
+"""
+Light cone halo mask functions.
+
+Author: Britton Smith <brittons at origins.colorado.edu>
+Affiliation: CASA/University of CO, Boulder
+Homepage: http://yt.enzotools.org/
+License:
+  Copyright (C) 2008 Britton Smith.  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.extensions.HaloProfiler import *
+import yt.lagos as lagos
+import copy
+import numpy as na
+import tables as h5
+
+#### Note: assumption of box width 1.  I'll fix it someday.
+
+def MakeLightConeHaloMask(lightCone,HaloMaskParameterFile,cube_file=None,mask_file=None):
+    "Make a boolean mask to cut clusters out of light cone projections."
+
+    pixels = int(lightCone.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
+            lightCone.lightConeParameters['ImageResolutionInArcSeconds'])
+
+    lightConeMask = []
+
+    # Loop through files in light cone solution and get virial quantities.
+    for slice in lightCone.lightConeSolution:
+        hp = HaloProfiler(slice['filename'],HaloMaskParameterFile)
+        hp._LoadVirialData()
+
+        lightConeMask.append(_MakeHaloMask(slice,hp.virialQuantities,pixels))
+
+    # Write out cube of masks from each slice.
+    if cube_file is not None:
+        output = h5.openFile(cube_file,'a')
+        output.createArray("/",'haloMaskCube',na.array(lightConeMask))
+        output.close()
+
+    # Write out final mask.
+    if mask_file is not None:
+        # Final mask is simply the product of the mask from each slice.
+        finalMask = na.ones(shape=(pixels,pixels))
+        for mask in lightConeMask:
+            finalMask *= mask
+
+        output = h5.openFile(mask_file,'a')
+        output.createArray("/",'haloMask',na.array(finalMask))
+        output.close()
+
+    return lightConeMask
+
+def _MakeHaloMask(slice,virialQuantities,pixels):
+    "Make halo mask for one slice in light cone solution."
+
+    # Make numpy arrays for halo centers and virial radii.
+    halo_x = []
+    halo_y = []
+    halo_depth = []
+    halo_radius = []
+
+    # Get units to convert virial radii to code units.
+    dataset_object = lagos.EnzoStaticOutput(slice['filename'])
+    Mpc_units = dataset_object.units['mpc']
+    del dataset_object
+
+    for halo in virialQuantities:
+        if halo is not None:
+            center = copy.deepcopy(halo['center'])
+            halo_depth.append(center.pop(slice['ProjectionAxis']))
+            halo_x.append(center[0])
+            halo_y.append(center[1])
+            halo_radius.append(halo['RadiusMpc']/Mpc_units)
+
+    halo_x = na.array(halo_x)
+    halo_y = na.array(halo_y)
+    halo_depth = na.array(halo_depth)
+    halo_radius = na.array(halo_radius)
+
+    # Adjust halo centers along line of sight.
+    depthCenter = slice['ProjectionCenter'][slice['ProjectionAxis']]
+    depthLeft = depthCenter - 0.5 * slice['DepthBoxFraction']
+    depthRight = depthCenter + 0.5 * slice['DepthBoxFraction']
+
+    # Make boolean mask to pick out centers in region along line of sight.
+    # Halos near edges may wrap around to other side.
+    add_left = (halo_depth + halo_radius) > 1 # should be box width
+    add_right = (halo_depth - halo_radius) < 0
+
+    halo_depth = na.concatenate([halo_depth,(halo_depth[add_left]-1),(halo_depth[add_right]+1)])
+    halo_x = na.concatenate([halo_x,halo_x[add_left],halo_x[add_right]])
+    halo_y = na.concatenate([halo_y,halo_y[add_left],halo_y[add_right]])
+    halo_radius = na.concatenate([halo_radius,halo_radius[add_left],halo_radius[add_right]])
+
+    del add_left,add_right
+
+    # Cut out the halos outside the region of interest.
+    if (slice['DepthBoxFraction'] < 1):
+        if (depthLeft < 0):
+            mask = ((halo_depth + halo_radius >= 0) & (halo_depth - halo_radius <= depthRight)) | \
+                ((halo_depth + halo_radius >= depthLeft + 1) & (halo_depth - halo_radius <= 1))
+        elif (depthRight > 1):
+            mask = ((halo_depth + halo_radius >= 0) & (halo_depth - halo_radius <= depthRight - 1)) | \
+                ((halo_depth + halo_radius >= depthLeft) & (halo_depth - halo_radius <= 1))
+        else:
+            mask = (halo_depth + halo_radius >= depthLeft) & (halo_depth - halo_radius <= depthRight)
+
+        halo_x = halo_x[mask]
+        halo_y = halo_y[mask]
+        halo_radius = halo_radius[mask]
+        del mask
+    del halo_depth
+
+    all_halo_x = na.array([])
+    all_halo_y = na.array([])
+    all_halo_radius = na.array([])
+
+    # Tile halos of width box fraction is greater than one.
+    # Copy original into offset positions to make tiles.
+    for x in range(int(na.ceil(slice['WidthBoxFraction']))):
+        for y in range(int(na.ceil(slice['WidthBoxFraction']))):
+            all_halo_x = na.concatenate([all_halo_x,halo_x+x])
+            all_halo_y = na.concatenate([all_halo_y,halo_y+y])
+            all_halo_radius = na.concatenate([all_halo_radius,halo_radius])
+
+    del halo_x,halo_y,halo_radius
+
+    # Shift centers laterally.
+    offset = copy.deepcopy(slice['ProjectionCenter'])
+    del offset[slice['ProjectionAxis']]
+
+    # Shift x and y positions.
+    all_halo_x -= offset[0]
+    all_halo_y -= offset[1]
+
+    # Wrap off-edge centers back around to other side (periodic boundary conditions).
+    all_halo_x[all_halo_x < 0] += na.ceil(slice['WidthBoxFraction'])
+    all_halo_y[all_halo_y < 0] += na.ceil(slice['WidthBoxFraction'])
+
+    # After shifting, some centers have fractional coverage on both sides of the box.
+    # Find those centers and make copies to be placed on the other side.
+
+    # Centers hanging off the right edge.
+    add_x_right = all_halo_x + all_halo_radius > na.ceil(slice['WidthBoxFraction'])
+    add_x_halo_x = all_halo_x[add_x_right]
+    add_x_halo_x -= na.ceil(slice['WidthBoxFraction'])
+    add_x_halo_y = all_halo_y[add_x_right]
+    add_x_halo_radius = all_halo_radius[add_x_right]
+    del add_x_right
+
+    # Centers hanging off the left edge.
+    add_x_left = all_halo_x - all_halo_radius < 0
+    add2_x_halo_x = all_halo_x[add_x_left]
+    add2_x_halo_x += na.ceil(slice['WidthBoxFraction'])
+    add2_x_halo_y = all_halo_y[add_x_left]
+    add2_x_halo_radius = all_halo_radius[add_x_left]
+    del add_x_left
+
+    # Centers hanging off the top edge.
+    add_y_right = all_halo_y + all_halo_radius > na.ceil(slice['WidthBoxFraction'])
+    add_y_halo_x = all_halo_x[add_y_right]
+    add_y_halo_y = all_halo_y[add_y_right]
+    add_y_halo_y -= na.ceil(slice['WidthBoxFraction'])
+    add_y_halo_radius = all_halo_radius[add_y_right]
+    del add_y_right
+
+    # Centers hanging off the bottom edge.
+    add_y_left = all_halo_y - all_halo_radius < 0
+    add2_y_halo_x = all_halo_x[add_y_left]
+    add2_y_halo_y = all_halo_y[add_y_left]
+    add2_y_halo_y += na.ceil(slice['WidthBoxFraction'])
+    add2_y_halo_radius = all_halo_radius[add_y_left]
+    del add_y_left
+
+    # Add the hanging centers back to the projection data.
+    all_halo_x = na.concatenate([all_halo_x,add_x_halo_x,add2_x_halo_x,add_y_halo_x,add2_y_halo_x])
+    all_halo_y = na.concatenate([all_halo_y,add_x_halo_y,add2_x_halo_y,add_y_halo_y,add2_y_halo_y])
+    all_halo_radius = na.concatenate([all_halo_radius,add_x_halo_radius,add2_x_halo_radius,
+                                      add_y_halo_radius,add2_y_halo_radius])
+
+    del add_x_halo_x,add_x_halo_y,add_x_halo_radius
+    del add2_x_halo_x,add2_x_halo_y,add2_x_halo_radius
+    del add_y_halo_x,add_y_halo_y,add_y_halo_radius
+    del add2_y_halo_x,add2_y_halo_y,add2_y_halo_radius
+
+    # Cut edges to proper width.
+    cut_mask = (all_halo_x - all_halo_radius < slice['WidthBoxFraction']) & \
+        (all_halo_y - all_halo_radius < slice['WidthBoxFraction'])
+    all_halo_x = all_halo_x[cut_mask]
+    all_halo_y = all_halo_y[cut_mask]
+    all_halo_radius = all_halo_radius[cut_mask]
+    del cut_mask
+
+    # Make boolean mask and cut out halos.
+    dx = slice['WidthBoxFraction'] / pixels
+    x = [(q + 0.5) * dx for q in range(pixels)]
+    haloMask = na.ones(shape=(pixels,pixels),dtype=bool)
+
+    # Cut out any pixel that has any part at all in the circle.
+    for q in range(len(all_halo_radius)):
+        dif_xIndex = na.array(int(all_halo_x[q]/dx) - na.array(range(pixels))) != 0
+        dif_yIndex = na.array(int(all_halo_y[q]/dx) - na.array(range(pixels))) != 0
+
+        xDistance = (na.abs(x - all_halo_x[q]) - (0.5 * dx)) * dif_xIndex
+        yDistance = (na.abs(x - all_halo_y[q]) - (0.5 * dx)) * dif_yIndex
+
+        distance = na.array([na.sqrt(w**2 + xDistance**2) for w in yDistance])
+        haloMask *= (distance >= all_halo_radius[q])
+
+    return haloMask

Modified: trunk/yt/extensions/lightcone/LightCone.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightCone.py	(original)
+++ trunk/yt/extensions/lightcone/LightCone.py	Mon Dec  8 21:12:29 2008
@@ -25,12 +25,13 @@
 
 from yt.extensions.lightcone import *
 from yt.logger import lagosLogger as mylog
+from Common_nVolume import *
+from HaloMask import *
 import copy
 import os
 import numpy as na
 import random as rand
 import tables as h5
-from Common_nVolume import *
 
 class LightCone(object):
     def __init__(self,EnzoParameterFile,LightConeParameterFile,verbose=True):
@@ -46,6 +47,7 @@
         self.masterSolution = [] # kept to compare with recycled solutions
         self.projectionStack = []
         self.projectionWeightFieldStack = []
+        self.haloMask = []
 
         # Parameters for recycling light cone solutions.
         self.recycleSolution = False
@@ -58,6 +60,10 @@
         self._ReadLightConeParameterFile()
         self._ReadEnzoParameterFile()
 
+        # Calculate number of pixels.
+        self.pixels = int(self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
+                              self.lightConeParameters['ImageResolutionInArcSeconds'])
+
         # Create output directory.
         if (os.path.exists(self.lightConeParameters['OutputDir'])):
             if not(os.path.isdir(self.lightConeParameters['OutputDir'])):
@@ -82,6 +88,9 @@
         # Make sure recycling flag is off.
         self.recycleSolution = False
 
+        # Get rid of old halo mask, if one was there.
+        self.haloMask = []
+
         if seed is not None:
             self.lightConeParameters['RandomSeed'] = int(seed)
 
@@ -158,9 +167,9 @@
             # Simple error check to make sure more than 100% of box depth is never required.
             if (self.lightConeSolution[q]['DepthBoxFraction'] > 1.0):
                 if self.verbose: mylog.error("Warning: box fraction required to go from z = %f to %f is %f" % (self.lightConeSolution[q]['redshift'],z_next,
-                                                                self.lightConeSolution[q]['DepthBoxFraction']))
+                                                                                                               self.lightConeSolution[q]['DepthBoxFraction']))
                 if self.verbose: mylog.error("Full box delta z is %f, but it is %f to the next data dump." % (self.lightConeSolution[q]['deltaz'],
-                                                                                       self.lightConeSolution[q]['redshift']-z_next))
+                                                                                                              self.lightConeSolution[q]['redshift']-z_next))
 
             # Calculate fraction of box required for width corresponding to requested image size.
             scale = co.AngularScale_1arcsec_kpc(self.lightConeParameters['FinalRedshift'],self.lightConeSolution[q]['redshift'])
@@ -182,7 +191,25 @@
         del self.timeOutputs
         del co
 
-    def ProjectLightCone(self,field,weight_field=None,**kwargs):
+    def GetHaloMask(self,HaloMaskParameterFile,mask_file=None,**kwargs):
+        "Gets a halo mask from a file or makes a new one."
+
+        # Check if file already exists.
+        if (mask_file is not None) and os.path.exists(mask_file):
+            input = h5.openFile(mask_file,'r')
+            self.haloMask = input.root.haloMask.read()
+            input.close()
+
+        # Otherwise, make a halo mask.
+        else:
+            haloMaskCube = MakeLightConeHaloMask(self,HaloMaskParameterFile,mask_file=mask_file,**kwargs)
+            # Collapse cube into final mask.
+            self.haloMask = na.ones(shape=(self.pixels,self.pixels),dtype=bool)
+            for mask in haloMaskCube:
+                self.haloMask *= mask
+            del haloMaskCube
+
+    def ProjectLightCone(self,field,weight_field=None,apply_halo_mask=False,**kwargs):
         "Create projections for light cone, then add them together."
 
         # Clear projection stack.
@@ -194,14 +221,11 @@
         if not(self.lightConeParameters['OutputDir'].endswith("/")):
                  self.lightConeParameters['OutputDir'] += "/"
 
-        pixels = int(self.lightConeParameters['FieldOfViewInArcMinutes'] * 60.0 / \
-            self.lightConeParameters['ImageResolutionInArcSeconds'])
-
         for q,output in enumerate(self.lightConeSolution):
             name = "%s%s_%04d_%04d" % (self.lightConeParameters['OutputDir'],self.lightConeParameters['OutputPrefix'],
                                        q,len(self.lightConeSolution))
             output['object'] = lagos.EnzoStaticOutput(output['filename'])
-            frb = LightConeProjection(output,field,pixels,weight_field=weight_field,
+            frb = LightConeProjection(output,field,self.pixels,weight_field=weight_field,
                                       save_image=self.lightConeParameters['SaveLightConeSlices'],
                                       name=name,**kwargs)
             if (weight_field is not None):
@@ -229,6 +253,14 @@
         # but replace the data with the full light cone projection data.
         frb.data[field] = lightConeProjection
 
+        # Apply halo mask.
+        if apply_halo_mask:
+            if len(self.haloMask) > 0:
+                mylog.info("Applying halo mask.")
+                frb.data[field] *= self.haloMask
+            else:
+                mylog.error("No halo mask loaded, call GetHaloMask.")
+
         # Make a plot collection for the light cone projection.
         center = [0.5 * (self.lightConeSolution[-1]['object'].parameters['DomainLeftEdge'][w] + 
                          self.lightConeSolution[-1]['object'].parameters['DomainRightEdge'][w])
@@ -254,6 +286,9 @@
         True, then it will recycle.  Otherwise, it will create a completely new solution.
         """
 
+        # Get rid of old halo mask, if one was there.
+        self.haloMask = []
+
         if recycle:
             if self.verbose: mylog.info("Recycling solution made with %s with new seed %s." % (self.lightConeParameters['RandomSeed'],
                                                                               newSeed))

Modified: trunk/yt/extensions/lightcone/LightConeProjection.py
==============================================================================
--- trunk/yt/extensions/lightcone/LightConeProjection.py	(original)
+++ trunk/yt/extensions/lightcone/LightConeProjection.py	Mon Dec  8 21:12:29 2008
@@ -30,7 +30,12 @@
 import numpy as na
 import copy
 
-def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name="",**kwargs):
+#### Note: There is an assumption that the box width is 1 in every direction.  This doesn't have to be, but 
+####       I don't have time to fix it now.  All that is required is a multiplication by the box width in 
+####       the direction in question wherever DepthBoxFraction or WidthBoxFraction appears.
+####       To be fixed before turning 30.  - Britton
+
+def LightConeProjection(lightConeSlice,field,pixels,weight_field=None,save_image=False,name="",field_cuts=None,**kwargs):
     "Create a single projection to be added into the light cone stack."
 
     # Use some projection parameters to seed random number generator to make unique node name.
@@ -53,41 +58,64 @@
 
     # 1. The Depth Problem
     # Excise a region from the box with the specified depth and make projection with that.
-    region_center = [0.5,0.5,0.5]
-
-    region_left = [0.0,0.0,0.0]
-    region_left[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] - \
-        0.5 * lightConeSlice['DepthBoxFraction']
+    region_left = copy.deepcopy(dataset_object.parameters['DomainLeftEdge'])
+    region_right = copy.deepcopy(dataset_object.parameters['DomainRightEdge'])
+    region_center = [0.5 * (region_right[q] - region_left[q]) for q in range(len(region_left))]
+
+    # Use coordinate field cut in line of sight to cut projection to proper depth.
+    if (lightConeSlice['DepthBoxFraction'] < 1):
+        axis = ('x','y','z')[lightConeSlice['ProjectionAxis']]
+        depthLeft = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] - 0.5 * lightConeSlice['DepthBoxFraction']
+        depthRight = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] + 0.5 * lightConeSlice['DepthBoxFraction']
+        if (depthLeft < 0):
+            cut_mask = "((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= 0) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= %f)) | ((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= 1))" % \
+                (axis,axis,axis,axis,depthRight,axis,axis,(depthLeft+1),axis,axis)
+        elif (depthRight > 1):
+            cut_mask = "((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= 0) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= %f)) | ((grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"d%s\"] <= 1))" % \
+                (axis,axis,axis,axis,(depthRight-1),axis,axis,depthLeft,axis,axis)
+        else:
+            cut_mask = "(grid[\"%s\"] + 0.5*grid[\"d%s\"] >= %f) & (grid[\"%s\"] - 0.5*grid[\"%s\"] <= %f)" % (axis,axis,depthLeft,axis,axis,depthRight)
+
+        if field_cuts is None:
+           field_cuts = []
+        field_cuts.append(cut_mask)
+
+##### Periodic region.
+#     region_left[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] - \
+#         0.5 * lightConeSlice['DepthBoxFraction']
 
-    region_right = [1.0,1.0,1.0]
-    region_right[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] + \
-        0.5 * lightConeSlice['DepthBoxFraction']
+#     region_right[lightConeSlice['ProjectionAxis']] = lightConeSlice['ProjectionCenter'][lightConeSlice['ProjectionAxis']] + \
+#         0.5 * lightConeSlice['DepthBoxFraction']
 
-    periodic_region = dataset_object.h.periodic_region(region_center,region_left,region_right)
+#     periodic_region = dataset_object.h.periodic_region(region_center,region_left,region_right)
 
     # Make plot collection.
     pc = raven.PlotCollection(dataset_object,center=region_center)
 
     # Make projection.
-    pc.add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,source=periodic_region,use_colorbar=True,
+    pc.add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,field_cuts=field_cuts,use_colorbar=True,
                       node_name=node_name,**kwargs)
 
+##### Periodic region.
+#     pc.add_projection(field,lightConeSlice['ProjectionAxis'],weight_field=weight_field,source=periodic_region,use_colorbar=True,
+#                       node_name=node_name,**kwargs)
+
     # Serialize projection data.
     pc.plots[0].data._serialize(node_name,force=True)
 
     # Delete region, maybe save some ram.
-    del periodic_region
+    #del periodic_region
 
     # 2. The Tile Problem
     # Tile projection to specified width.
 
     # Original projection data.
-    original_px = pc.plots[0].data['px']
-    original_py = pc.plots[0].data['py']
-    original_pdx = pc.plots[0].data['pdx']
-    original_pdy = pc.plots[0].data['pdy']
-    original_field = pc.plots[0].data[field]
-    original_weight_field = pc.plots[0].data['weight_field']
+    original_px = copy.deepcopy(pc.plots[0].data['px'])
+    original_py = copy.deepcopy(pc.plots[0].data['py'])
+    original_pdx = copy.deepcopy(pc.plots[0].data['pdx'])
+    original_pdy = copy.deepcopy(pc.plots[0].data['pdy'])
+    original_field = copy.deepcopy(pc.plots[0].data[field])
+    original_weight_field = copy.deepcopy(pc.plots[0].data['weight_field'])
 
     # Copy original into offset positions to make tiles.
     for x in range(int(na.ceil(lightConeSlice['WidthBoxFraction']))):
@@ -127,40 +155,48 @@
     # Find those cells and make copies to be placed on the other side.
 
     # Cells hanging off the right edge.
-    add_x_px = pc.plots[0]['px'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+    add_x_right = pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])
+    add_x_px = pc.plots[0]['px'][add_x_right]
     add_x_px -= na.ceil(lightConeSlice['WidthBoxFraction'])
-    add_x_py = pc.plots[0]['py'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_x_pdx = pc.plots[0]['pdx'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_x_pdy = pc.plots[0]['pdy'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_x_field = pc.plots[0][field][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_x_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['px'] + 0.5 * pc.plots[0]['pdx'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+    add_x_py = pc.plots[0]['py'][add_x_right]
+    add_x_pdx = pc.plots[0]['pdx'][add_x_right]
+    add_x_pdy = pc.plots[0]['pdy'][add_x_right]
+    add_x_field = pc.plots[0][field][add_x_right]
+    add_x_weight_field = pc.plots[0]['weight_field'][add_x_right]
+    del add_x_right
 
     # Cells hanging off the left edge.
-    add2_x_px = pc.plots[0]['px'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+    add_x_left = pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0
+    add2_x_px = pc.plots[0]['px'][add_x_left]
     add2_x_px += na.ceil(lightConeSlice['WidthBoxFraction'])
-    add2_x_py = pc.plots[0]['py'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
-    add2_x_pdx = pc.plots[0]['pdx'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
-    add2_x_pdy = pc.plots[0]['pdy'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
-    add2_x_field = pc.plots[0][field][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
-    add2_x_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['px'] - 0.5 * pc.plots[0]['pdx'] < 0]
+    add2_x_py = pc.plots[0]['py'][add_x_left]
+    add2_x_pdx = pc.plots[0]['pdx'][add_x_left]
+    add2_x_pdy = pc.plots[0]['pdy'][add_x_left]
+    add2_x_field = pc.plots[0][field][add_x_left]
+    add2_x_weight_field = pc.plots[0]['weight_field'][add_x_left]
+    del add_x_left
 
     # Cells hanging off the top edge.
-    add_y_px = pc.plots[0]['px'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_y_py = pc.plots[0]['py'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+    add_y_right = pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])
+    add_y_px = pc.plots[0]['px'][add_y_right]
+    add_y_py = pc.plots[0]['py'][add_y_right]
     add_y_py -= na.ceil(lightConeSlice['WidthBoxFraction'])
-    add_y_pdx = pc.plots[0]['pdx'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_y_pdy = pc.plots[0]['pdy'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_y_field = pc.plots[0][field][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
-    add_y_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['py'] + 0.5 * pc.plots[0]['pdy'] > na.ceil(lightConeSlice['WidthBoxFraction'])]
+    add_y_pdx = pc.plots[0]['pdx'][add_y_right]
+    add_y_pdy = pc.plots[0]['pdy'][add_y_right]
+    add_y_field = pc.plots[0][field][add_y_right]
+    add_y_weight_field = pc.plots[0]['weight_field'][add_y_right]
+    del add_y_right
 
     # Cells hanging off the bottom edge.
-    add2_y_px = pc.plots[0]['px'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
-    add2_y_py = pc.plots[0]['py'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+    add_y_left = pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0
+    add2_y_px = pc.plots[0]['px'][add_y_left]
+    add2_y_py = pc.plots[0]['py'][add_y_left]
     add2_y_py += na.ceil(lightConeSlice['WidthBoxFraction'])
-    add2_y_pdx = pc.plots[0]['pdx'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
-    add2_y_pdy = pc.plots[0]['pdy'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
-    add2_y_field = pc.plots[0][field][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
-    add2_y_weight_field = pc.plots[0]['weight_field'][pc.plots[0]['py'] - 0.5 * pc.plots[0]['pdy'] < 0]
+    add2_y_pdx = pc.plots[0]['pdx'][add_y_left]
+    add2_y_pdy = pc.plots[0]['pdy'][add_y_left]
+    add2_y_field = pc.plots[0][field][add_y_left]
+    add2_y_weight_field = pc.plots[0]['weight_field'][add_y_left]
+    del add_y_left
 
     # Add the hanging cells back to the projection data.
     pc.plots[0].data['px'] = na.concatenate([pc.plots[0]['px'],add_x_px,add_y_px,add2_x_px,add2_y_px])
@@ -168,8 +204,8 @@
     pc.plots[0].data['pdx'] = na.concatenate([pc.plots[0]['pdx'],add_x_pdx,add_y_pdx,add2_x_pdx,add2_y_pdx])
     pc.plots[0].data['pdy'] = na.concatenate([pc.plots[0]['pdy'],add_x_pdy,add_y_pdy,add2_x_pdy,add2_y_pdy])
     pc.plots[0].data[field] = na.concatenate([pc.plots[0][field],add_x_field,add_y_field,add2_x_field,add2_y_field])
-    pc.plots[0].data['weight_field'] = na.concatenate([pc.plots[0]['weight_field'],
-                                                               add_x_weight_field,add_y_weight_field,add2_x_weight_field,add2_y_weight_field])
+    pc.plots[0].data['weight_field'] = na.concatenate([pc.plots[0]['weight_field'],add_x_weight_field,add_y_weight_field,
+                                                       add2_x_weight_field,add2_y_weight_field])
 
     # Delete original copies of hanging cells.
     del add_x_px,add_y_px,add2_x_px,add2_y_px
@@ -191,6 +227,6 @@
     # Create fixed resolution buffer to return back to the light cone object.
     # These buffers will be stacked together to make the light cone.
     frb = raven.FixedResolutionBuffer(pc.plots[0].data,(0,lightConeSlice['WidthBoxFraction'],0,lightConeSlice['WidthBoxFraction']),
-                                      (pixels,pixels),antialias=True)
+                                      (pixels,pixels),antialias=False)
 
     return frb

Modified: trunk/yt/extensions/lightcone/__init__.py
==============================================================================
--- trunk/yt/extensions/lightcone/__init__.py	(original)
+++ trunk/yt/extensions/lightcone/__init__.py	Mon Dec  8 21:12:29 2008
@@ -27,5 +27,5 @@
 
 from LightConeProjection import *
 from LightCone import *
+from HaloMask import *
 from UniqueSolution import *
-



More information about the yt-svn mailing list