Changeset 1618

Timestamp:

02/11/10 11:47:43 (6 months ago)

Author:

mturk

Message:

Bug fixes from trunk and instinfo command.

Location:

branches/yt-1.6/yt

Files:

• commands.py (modified) (3 diffs)
• extensions/HaloMassFcn.py (modified) (2 diffs)
• extensions/volume_rendering/software_sampler.py (modified) (3 diffs)
• lagos/EnzoDefs.py (modified) (1 diff)
• lagos/EnzoFields.py (modified) (1 diff)
• lagos/HDF5LightReader.c (modified) (1 diff)
• lagos/HierarchyType.py (modified) (1 diff)
• lagos/parallelHOP/parallelHOP.py (modified) (2 diffs)
• lagos/Profiles.py (modified) (1 diff)
• _amr_utils/FixedInterpolator.c (modified) (1 diff)
• _amr_utils/VolumeIntegrator.pyx (modified) (13 diffs)

branches/yt-1.6/yt/commands.py

@@ -28 +28 @@
 from yt.recipes import _fix_pf
 import yt.cmdln as cmdln
-import optparse, os, os.path, math, sys
+import optparse, os, os.path, math, sys, time, subprocess
 
 _common_options = dict(
@@ -206 +206 @@
     return arg_iterate
 
+def _get_vcs_type(path):
+    if os.path.exists(os.path.join(path, ".hg")):
+        return "hg"
+    if os.path.exists(os.path.join(path, ".svn")):
+        return "svn"
+    return None
+
+def _get_svn_version(path):
+    p = subprocess.Popen(["svn", "info", path], stdout = subprocess.PIPE,
+                                               stderr = subprocess.STDOUT)
+    stdout, stderr = p.communicate()
+    return stdout
+
+def _update_svn(path):
+    f = open(os.path.join(path, "yt_updater.log"), "a")
+    f.write("\n\nUPDATE PROCESS: %s\n\n" % (time.asctime()))
+    p = subprocess.Popen(["svn", "up", path], stdout = subprocess.PIPE,
+                                              stderr = subprocess.STDOUT)
+    stdout, stderr = p.communicate()
+    f.write(stdout)
+    f.write("\n\n")
+    if p.returncode:
+        print "BROKEN: See %s" % (os.path.join(path, "yt_updater.log"))
+        sys.exit(1)
+    f.write("Rebuilding modules\n\n")
+    p = subprocess.Popen([sys.executable, "setup.py", "build_ext", "-i"], cwd=path,
+                        stdout = subprocess.PIPE, stderr = subprocess.STDOUT)
+    stdout, stderr = p.communicate()
+    f.write(stdout)
+    f.write("\n\n")
+    if p.returncode:
+        print "BROKEN: See %s" % (os.path.join(path, "yt_updater.log"))
+        sys.exit(1)
+    f.write("Successful!\n")
+
+def _update_hg(path):
+    from mercurial import hg, ui, commands
+    f = open(os.path.join(path, "yt_updater.log"), "a")
+    u = ui.ui()
+    u.pushbuffer()
+    config_fn = os.path.join(path, ".hg", "hgrc")
+    print "Reading configuration from ", config_fn
+    u.readconfig(config_fn)
+    repo = hg.repository(u, path)
+    commands.pull(u, repo)
+    f.write(u.popbuffer())
+    f.write("\n\n")
+    u.pushbuffer()
+    commands.identify(u, repo)
+    if "+" in u.popbuffer():
+        print "Can't rebuild modules by myself."
+        print "You will have to do this yourself.  Here's a sample commands:"
+        print
+        print "    $ cd %s" % (path)
+        print "    $ hg up"
+        print "    $ %s setup.py develop" % (sys.executable)
+        sys.exit(1)
+    f.write("Rebuilding modules\n\n")
+    p = subprocess.Popen([sys.executable, "setup.py", "build_ext", "-i"], cwd=path,
+                        stdout = subprocess.PIPE, stderr = subprocess.STDOUT)
+    stdout, stderr = p.communicate()
+    f.write(stdout)
+    f.write("\n\n")
+    if p.returncode:
+        print "BROKEN: See %s" % (os.path.join(path, "yt_updater.log"))
+        sys.exit(1)
+    f.write("Successful!\n")
+
+def _get_hg_version(path):
+    from mercurial import hg, ui, commands
+    u = ui.ui()
+    u.pushbuffer()
+    repo = hg.repository(u, path)
+    commands.identify(u, repo)
+    return u.popbuffer()
+
+_vcs_identifier = dict(svn = _get_svn_version,
+                        hg = _get_hg_version)
+_vcs_updater = dict(svn = _update_svn,
+                     hg = _update_hg)
+
 class YTCommands(cmdln.Cmdln):
     name="yt"
@@ -220 +301 @@
         """
         self.cmdloop()
+
+    @cmdln.option("-u", "--update-source", action="store_true",
+                  default = False,
+                  help="Update the yt installation, if able")
+    def do_instinfo(self, subcmd, opts):
+        """
+        ${cmd_name}: Get some information about the yt installation
+
+        ${cmd_usage}
+        ${cmd_option_list}
+        """
+        import pkg_resources
+        yt_provider = pkg_resources.get_provider("yt")
+        path = os.path.dirname(yt_provider.module_path)
+        print
+        print "yt module located at:"
+        print "    %s" % (path)
+        if "site-packages" not in path:
+            vc_type = _get_vcs_type(path)
+            vstring = _vcs_identifier[vc_type](path)
+            print
+            print "The current version of the code is:"
+            print
+            print "---"
+            print vstring
+            print "---"
+            print
+            print "This installation CAN be automatically updated."
+            if opts.update_source:  
+                _vcs_updater[vc_type](path)
+            print "Updated successfully."
 
     @add_cmd_options(['outputfn','bn','thresh','dm_only','skip'])

branches/yt-1.6/yt/extensions/HaloMassFcn.py

@@ -179 +179 @@
         mylog.info("Reading halo masses from %s" % self.halo_file)
         f = open(self.halo_file,'r')
-        line = f.readline() # burn the top header line.
         line = f.readline()
+        while line[0] == '#':
+            line = f.readline()
         self.haloes = []
         while line:
@@ -204 +205 @@
             if i == (self.num_sigma_bins - 3): break
 
-        self.dis = dis / self.pf['CosmologyComovingBoxSize']**3.0
+        self.dis = dis  / self.pf['CosmologyComovingBoxSize']**3.0 * self.hubble0**3.0
 
     def sigmaM(self):

branches/yt-1.6/yt/extensions/volume_rendering/software_sampler.py

@@ -29 +29 @@
 
 def direct_ray_cast(pf, L, center, W, Nvec, tf, 
-                    partitioned_grids = None, field = 'Density', log_field = True, whole_box=False):
+                    partitioned_grids = None, field = 'Density',
+                    log_field = True, whole_box=False,
+                    nsamples = 5):
     center = na.array(center, dtype='float64')
 
@@ -73 +75 @@
     hit = 0
     tnow = time.time()
-    every = na.ceil(len(partitioned_grids) / 100.0)
 
     vp = VectorPlane(vectors, norm_vec, back_center,
@@ -86 +87 @@
     
     tfp = TransferFunctionProxy(tf)
+    tfp.ns = nsamples
 
-    pbar = get_pbar("Ray casting ", len(partitioned_grids))
+
+    total_cells = sum(na.prod(g.my_data.shape) for g in partitioned_grids)
+    pbar = get_pbar("Ray casting ", total_cells)
+    total_cells = 0
     for i,g in enumerate(partitioned_grids[ind]):
-        if (i % every) == 0: 
-            pbar.update(i)
+        pbar.update(total_cells)
         pos = g.cast_plane(tfp, vp)
+        total_cells += na.prod(g.my_data.shape)
     pbar.finish()
 

branches/yt-1.6/yt/lagos/EnzoDefs.py

@@ -86 +86 @@
                   'pc'    : 1e6,
                   'au'    : 2.063e11,
-                  'rsun'  : 2.2167e13,
+                  'rsun'  : 4.43664e13,
                   'cm'    : 3.0857e24,
                   'miles' : 1.917e19}

branches/yt-1.6/yt/lagos/EnzoFields.py

@@ -343 +343 @@
           projection_conversion="1")
 
+def _IsStarParticle(field, data):
+    is_star = (data['creation_time'] > 0).astype('float64')
+    return is_star
+add_field('IsStarParticle', function=_IsStarParticle,
+          particle_type = True)
+
 #
 # Now we do overrides for 2D fields

branches/yt-1.6/yt/lagos/HDF5LightReader.c

@@ -34 +34 @@
 
 #include "numpy/ndarrayobject.h"
+
+#ifndef npy_float128
+#define npy_float128 npy_longdouble
+#endif
 
 

branches/yt-1.6/yt/lagos/HierarchyType.py

@@ -250 +250 @@
 
         full_name = "%s/%s" % (node, name)
-        return self._data_file[full_name][:]
+        try:
+            return self._data_file[full_name][:]
+        except TypeError:
+            return self._data_file[full_name]
 
     def _close_data_file(self):

branches/yt-1.6/yt/lagos/parallelHOP/parallelHOP.py

@@ -105 +105 @@
                 # max_padding, which is more likely the case with load-balancing
                 # turned on.
-                dx = na.min( na.abs(my_vertex[0] - vertex[0]), \
-                    self.period[0] - na.abs(my_vertex[0] - vertex[0]))
-                dy = na.min( na.abs(my_vertex[1] - vertex[1]), \
-                    self.period[1] - na.abs(my_vertex[1] - vertex[1]))
-                dz = na.min( na.abs(my_vertex[2] - vertex[2]), \
-                    self.period[2] - na.abs(my_vertex[2] - vertex[2]))
+                dx = min( na.fabs(my_vertex[0] - vertex[0]), \
+                    self.period[0] - na.fabs(my_vertex[0] - vertex[0]))
+                dy = min( na.fabs(my_vertex[1] - vertex[1]), \
+                    self.period[1] - na.fabs(my_vertex[1] - vertex[1]))
+                dz = min( na.fabs(my_vertex[2] - vertex[2]), \
+                    self.period[2] - na.fabs(my_vertex[2] - vertex[2]))
                 d = na.sqrt(dx*dx + dy*dy + dz*dz)
                 if d <= self.max_padding:
@@ -1375 +1375 @@
         if calc:
             com = self.CoM[subchain]
-            rad = na.abs(com - loc)
+            rad = na.fabs(com - loc)
             dist = (na.minimum(rad, self.period - rad)**2.).sum(axis=1)
             dist = dist[sort]

branches/yt-1.6/yt/lagos/Profiles.py

@@ -167 +167 @@
                 # in the way.
                 if field in self.pf.field_info \
-                    and self.pf.field_info[field].particle_type \
-                    and not self._data_source._is_fully_enclosed(source):
-                    pointI = self._data_source._get_particle_indices(source)
+                    and self.pf.field_info[field].particle_type:
+                    if not self._data_source._is_fully_enclosed(source):
+                        pointI = self._data_source._get_particle_indices(source)
                 else:
                     pointI = self._data_source._get_point_indices(source)

branches/yt-1.6/yt/_amr_utils/FixedInterpolator.c

@@ -98 +98 @@
       normval += grad[i]*grad[i];
     }
-    normval = sqrt(normval);
-    for (i = 0; i < 3; i++) grad[i] /= -normval;
-    //fprintf(stderr, "Normval: %0.3lf %0.3lf %0.3lf %0.3lf\n",
-    //        normval, grad[0], grad[1], grad[2]);
+    if (normval != 0.0){
+      normval = sqrt(normval);
+      for (i = 0; i < 3; i++) grad[i] /= -normval;
+      //fprintf(stderr, "Normval: %0.3lf %0.3lf %0.3lf %0.3lf\n",
+      //        normval, grad[0], grad[1], grad[2]);
+    }else{
+      grad[0]=grad[1]=grad[2]=0.0;
+    }
 }

branches/yt-1.6/yt/_amr_utils/VolumeIntegrator.pyx

@@ -46 +46 @@
 
 cdef inline int iclip(int i, int a, int b):
-    return imin(imax(i, a), b)
+    if i < a: return a
+    if i > b: return b
+    return i
 
 cdef inline np.float64_t fclip(np.float64_t f,
@@ -76 +78 @@
     cdef np.float64_t *vs[4]
     cdef int nbins
-    cdef np.float64_t dbin
+    cdef public int ns
+    cdef np.float64_t dbin, idbin
     cdef np.float64_t light_color[3]
     cdef np.float64_t light_dir[3]
@@ -96 +99 @@
         self.nbins = tf_obj.nbins
         self.dbin = (self.x_bounds[1] - self.x_bounds[0])/self.nbins
+        self.idbin = 1.0/self.dbin
         self.light_color[0] = tf_obj.light_color[0]
         self.light_color[1] = tf_obj.light_color[1]
@@ -108 +112 @@
         self.use_light = tf_obj.use_light
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cdef void interpolate(self, np.float64_t dv, np.float64_t *trgba):
+        cdef int bin_id, channel
+        cdef np.float64_t bv, dy, dd, tf
+        bin_id = <int> ((dv - self.x_bounds[0]) * self.idbin)
+        # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+        dd = dv-(self.x_bounds[0] + bin_id * self.dbin) # x - x0
+        for channel in range(4):
+            bv = self.vs[channel][bin_id] # This is x0
+            dy = self.vs[channel][bin_id+1]-bv # dy
+                # This is our final value for transfer function on the entering face
+            trgba[channel] = bv+dd*dy*self.idbin
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
     cdef void eval_transfer(self, np.float64_t dt, np.float64_t dv,
                                     np.float64_t *rgba, np.float64_t *grad):
         cdef int i
-        cdef int bin_id
-        cdef np.float64_t tf, trgba[4], bv, dx, dy, dd, ta, dot_prod
-        dx = self.dbin
-
+        cdef np.float64_t ta, tf, trgba[4], dot_prod
+        self.interpolate(dv, trgba) 
         # get source alpha first
         # First locate our points
-        bin_id = iclip(<int> floor((dv - self.x_bounds[0]) / dx),
-                        0, self.nbins-2)
-            # Recall that linear interpolation is y0 + (x-x0) * dx/dy
-        bv = self.vs[3][bin_id] # This is x0
-        dy = self.vs[3][bin_id+1]-bv # dy
-        dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
-            # This is our final value for transfer function on the entering face
-        tf = bv+dd*(dy/dx) 
-        ta = tf  # Store the source alpha
         dot_prod = 0.0
-        for i in range(3):
-            dot_prod += self.light_dir[i] * grad[i]
-        #print dot_prod, grad[0], grad[1], grad[2]
-        dot_prod = fmax(0.0, dot_prod)
-        for i in range(3):
-            # Recall that linear interpolation is y0 + (x-x0) * dx/dy
-            bv = self.vs[i][bin_id] # This is x0
-            dy = self.vs[i][bin_id+1]-bv # dy
-            dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
-            # This is our final value for transfer function on the entering face
-            tf = bv+dd*(dy/dx) + dot_prod * self.light_color[i]
-            # alpha blending
-            rgba[i] += (1. - rgba[3])*ta*tf*dt
-        #update alpha
-        rgba[3] += (1. - rgba[3])*ta*dt
-        # We should really do some alpha blending.
-        # Front to back blending is defined as:
-        #  dst.rgb = dst.rgb + (1 - dst.a) * src.a * src.rgb
-        #  dst.a   = dst.a   + (1 - dst.a) * src.a     
+        if self.use_light:
+            for i in range(3):
+                dot_prod += self.light_dir[i] * grad[i]
+            dot_prod = fmax(0.0, dot_prod)
+            for i in range(3):
+                trgba[i] += dot_prod*self.light_color[i]
+        # alpha blending
+        ta = (1.0 - rgba[3])*dt*trgba[3]
+        for i in range(4):
+            rgba[i] += ta*trgba[i]
 
 cdef class VectorPlane:
@@ -181 +181 @@
         self.pdy = (self.bounds[3] - self.bounds[2])/self.nv
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
     cdef void get_start_stop(self, np.float64_t *ex, int *rv):
         # Extrema need to be re-centered
         cdef np.float64_t cx, cy
+        cdef int i
         cx = cy = 0.0
         for i in range(3):
@@ -215 +218 @@
     cdef np.float64_t right_edge[3]
     cdef np.float64_t dds[3]
+    cdef np.float64_t idds[3]
     cdef public np.float64_t min_dds
-    cdef int ns
     cdef int dims[3]
 
@@ -235 +238 @@
             self.dims[i] = dims[i]
             self.dds[i] = (self.right_edge[i] - self.left_edge[i])/dims[i]
+            self.idds[i] = 1.0/self.dds[i]
         self.my_data = data
         self.data = <np.float64_t*> data.data
@@ -245 +249 @@
         # like http://courses.csusm.edu/cs697exz/ray_box.htm
         cdef int vi, vj, hit, i, ni, nj, nn
-        self.ns = 5 #* (1 + <int> log2(self.dds[0] / self.min_dds))
         cdef int iter[4]
         cdef np.float64_t v_pos[3], v_dir[3], rgba[4], extrema[4]
@@ -295 +298 @@
         cdef np.float64_t intersect_t = 1.0
         cdef np.float64_t intersect[3], tmax[3], tdelta[3]
-        cdef np.float64_t enter_t, dist, alpha, dt
+        cdef np.float64_t enter_t, dist, alpha, dt, exit_t
         cdef np.float64_t tr, tl, temp_x, temp_y, dv
         for i in range(3):
@@ -329 +332 @@
             cur_ind[i] = <int> floor((intersect[i] +
                                       step[i]*1e-8*self.dds[i] -
-                                      self.left_edge[i])/self.dds[i])
+                                      self.left_edge[i])*self.idds[i])
             tmax[i] = (((cur_ind[i]+step[i])*self.dds[i])+
                         self.left_edge[i]-v_pos[i])/v_dir[i]
+            # This deals with the asymmetry in having our indices refer to the
+            # left edge of a cell, but the right edge of the brick being one
+            # extra zone out.
             if cur_ind[i] == self.dims[i] and step[i] < 0:
                 cur_ind[i] = self.dims[i] - 1
@@ -355 +361 @@
             if tmax[0] < tmax[1]:
                 if tmax[0] < tmax[2]:
-                    self.sample_values(v_pos, v_dir, enter_t, tmax[0], cur_ind,
+                    exit_t = fmin(tmax[0], 1.0)
+                    self.sample_values(v_pos, v_dir, enter_t, exit_t, cur_ind,
                                        rgba, tf)
                     cur_ind[0] += step[0]
-                    dt = fmin(tmax[0], 1.0) - enter_t
                     enter_t = tmax[0]
                     tmax[0] += tdelta[0]
                 else:
-                    self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+                    exit_t = fmin(tmax[2], 1.0)
+                    self.sample_values(v_pos, v_dir, enter_t, exit_t, cur_ind,
                                        rgba, tf)
                     cur_ind[2] += step[2]
-                    dt = fmin(tmax[2], 1.0) - enter_t
                     enter_t = tmax[2]
                     tmax[2] += tdelta[2]
             else:
                 if tmax[1] < tmax[2]:
-                    self.sample_values(v_pos, v_dir, enter_t, tmax[1], cur_ind,
+                    exit_t = fmin(tmax[1], 1.0)
+                    self.sample_values(v_pos, v_dir, enter_t, exit_t, cur_ind,
                                        rgba, tf)
                     cur_ind[1] += step[1]
-                    dt = fmin(tmax[1], 1.0) - enter_t
                     enter_t = tmax[1]
                     tmax[1] += tdelta[1]
                 else:
-                    self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+                    exit_t = fmin(tmax[2], 1.0)
+                    self.sample_values(v_pos, v_dir, enter_t, exit_t, cur_ind,
                                        rgba, tf)
                     cur_ind[2] += step[2]
-                    dt = fmin(tmax[2], 1.0) - enter_t
                     enter_t = tmax[2]
                     tmax[2] += tdelta[2]
@@ -386 +392 @@
         return hit
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
     cdef void sample_values(self,
                             np.float64_t v_pos[3],
@@ -395 +403 @@
                             TransferFunctionProxy tf):
         cdef np.float64_t cp[3], dp[3], temp, dt, t, dv
-        cdef np.float64_t grad[3]
+        cdef np.float64_t grad[3], ds[3]
+        grad[0] = grad[1] = grad[2] = 0.0
         cdef int dti, i
-        dt = (exit_t - enter_t) / (self.ns-1) # five samples, so divide by four
-        for dti in range(self.ns - 1):
-            t = enter_t + dt * dti
+        dt = (exit_t - enter_t) / (tf.ns) # 4 samples should be dt=0.25
+        for i in range(3):
+            temp = ci[i] * self.dds[i] + self.left_edge[i]
+            dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i] - temp
+            dp[i] *= self.idds[i]
+            ds[i] = v_dir[i] * self.idds[i] * dt
+        for dti in range(tf.ns): 
             for i in range(3):
-                cp[i] = v_pos[i] + t * v_dir[i]
-                dp[i] = fclip(fmod(cp[i], self.dds[i])/self.dds[i], 0, 1.0)
+                dp[i] += ds[i]
             dv = trilinear_interpolate(self.dims, ci, dp, self.data)
+            if not ((dv > tf.x_bounds[0]) and (dv < tf.x_bounds[1])):
+                continue
             if tf.use_light == 1:
                 eval_gradient(self.dims, ci, dp, self.data, grad)