[Yt-svn] yt-commit r1616 - in trunk/yt: . _amr_utils extensions/volume_rendering lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Thu Feb 11 11:40:17 PST 2010
Author: mturk
Date: Thu Feb 11 11:40:12 2010
New Revision: 1616
URL: http://yt.enzotools.org/changeset/1616
Log:
Backporting from hg: new "instinfo" command, vast improvements and speedups to
the volume renderer, bug fix for longdouble on some arches.
Modified:
trunk/yt/_amr_utils/FixedInterpolator.c
trunk/yt/_amr_utils/VolumeIntegrator.pyx
trunk/yt/commands.py
trunk/yt/extensions/volume_rendering/software_sampler.py
trunk/yt/lagos/HDF5LightReader.c
Modified: trunk/yt/_amr_utils/FixedInterpolator.c
==============================================================================
--- trunk/yt/_amr_utils/FixedInterpolator.c (original)
+++ trunk/yt/_amr_utils/FixedInterpolator.c Thu Feb 11 11:40:12 2010
@@ -97,8 +97,12 @@
dp[i] = backup;
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;
+ }
}
Modified: trunk/yt/_amr_utils/VolumeIntegrator.pyx
==============================================================================
--- trunk/yt/_amr_utils/VolumeIntegrator.pyx (original)
+++ trunk/yt/_amr_utils/VolumeIntegrator.pyx Thu Feb 11 11:40:12 2010
@@ -45,7 +45,9 @@
return f1
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,
np.float64_t a, np.float64_t b):
@@ -76,7 +78,7 @@
cdef np.float64_t *vs[4]
cdef int nbins
cdef public int ns
- cdef np.float64_t dbin
+ cdef np.float64_t dbin, idbin
cdef np.float64_t light_color[3]
cdef np.float64_t light_dir[3]
cdef int use_light
@@ -96,6 +98,7 @@
self.x_bounds[1] = tf_obj.x_bounds[1]
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]
self.light_color[2] = tf_obj.light_color[2]
@@ -108,44 +111,40 @@
for i in range(3): self.light_dir[i] /= normval
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:
cdef public object avp_pos, avp_dir, acenter, aimage
@@ -181,9 +180,12 @@
self.pdx = (self.bounds[1] - self.bounds[0])/self.nv
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):
cx += self.center[i] * self.x_vec[i]
@@ -215,6 +217,7 @@
cdef np.float64_t left_edge[3]
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 dims[3]
@@ -234,6 +237,7 @@
self.right_edge[i] = right_edge[i]
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
@@ -293,7 +297,7 @@
cdef int cur_ind[3], step[3], x, y, i, n, flat_ind, hit, direction
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):
if (v_dir[i] < 0):
@@ -327,9 +331,12 @@
intersect[i] = v_pos[i] + intersect_t * v_dir[i]
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
if cur_ind[i] < 0 or cur_ind[i] >= self.dims[i]: return 0
@@ -353,37 +360,39 @@
hit += 1
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]
if enter_t > 1.0: break
return hit
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
cdef void sample_values(self,
np.float64_t v_pos[3],
np.float64_t v_dir[3],
@@ -393,15 +402,21 @@
np.float64_t *rgba,
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) / (tf.ns-1) # five samples, so divide by four
- for dti in range(tf.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)
tf.eval_transfer(dt, dv, rgba, grad)
Modified: trunk/yt/commands.py
==============================================================================
--- trunk/yt/commands.py (original)
+++ trunk/yt/commands.py Thu Feb 11 11:40:12 2010
@@ -27,7 +27,7 @@
from yt.funcs import *
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(
axis = dict(short="-a", long="--axis",
@@ -205,6 +205,87 @@
func(self, subcmd, opts, arg)
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,6 +301,37 @@
"""
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'])
@check_args
def do_hop(self, subcmd, opts, arg):
Modified: trunk/yt/extensions/volume_rendering/software_sampler.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/software_sampler.py (original)
+++ trunk/yt/extensions/volume_rendering/software_sampler.py Thu Feb 11 11:40:12 2010
@@ -74,7 +74,6 @@
norm_vec = cp._norm_vec * (2.0*W*na.sqrt(3))
hit = 0
tnow = time.time()
- every = na.ceil(len(partitioned_grids) / 100.0)
vp = VectorPlane(vectors, norm_vec, back_center,
(xp0, xp1, yp0, yp1), image, cp._x_vec, cp._y_vec)
@@ -89,11 +88,14 @@
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()
return partitioned_grids[ind], image, vectors, norm_vec, pos
Modified: trunk/yt/lagos/HDF5LightReader.c
==============================================================================
--- trunk/yt/lagos/HDF5LightReader.c (original)
+++ trunk/yt/lagos/HDF5LightReader.c Thu Feb 11 11:40:12 2010
@@ -34,6 +34,10 @@
#include "numpy/ndarrayobject.h"
+#ifndef npy_float128
+#define npy_float128 npy_longdouble
+#endif
+
static PyObject *_hdf5ReadError;
herr_t iterate_dataset(hid_t loc_id, const char *name, void *nodelist);
More information about the yt-svn
mailing list