[yt-svn] commit/yt: 4 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Dec 5 11:02:34 PST 2016
4 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/b850dff755db/
Changeset: b850dff755db
Branch: yt
User: xarthisius
Date: 2016-12-02 22:06:46+00:00
Summary: Use explicit multiplication or sqrt instead of pow(..., 2/0.5)
Affected #: 13 files
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -311,7 +311,7 @@
else:
self.mk[ii[2], ii[1], ii[0], offset] = mk + (fields[0] - mk) / k
self.qk[ii[2], ii[1], ii[0], offset] = \
- qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
+ qk + (k - 1.0) * (fields[0] - mk) * (fields[0] - mk) / k
self.i[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -1711,10 +1711,13 @@
cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil:
# this is the sphere selection
cdef int i
- cdef np.float64_t dist2 = 0
+ cdef np.float64_t dist, dist2_max, dist2 = 0
for i in range(3):
- dist2 += self.difference(pos[i], self.center[i], i)**2
- if dist2 <= (self.mag[0]+radius)**2: return 1
+ dist = self.difference(pos[i], self.center[i], i)
+ dist2 += dist * dist
+ dist2_max = (self.mag[0] + radius) * (self.mag[0] + radius)
+ if dist2 <= dist2_max:
+ return 1
return 0
@cython.boundscheck(False)
@@ -1724,7 +1727,7 @@
np.float64_t right_edge[3]) nogil:
# This is the sphere selection
cdef int i
- cdef np.float64_t box_center, relcenter, closest, dist, edge
+ cdef np.float64_t box_center, relcenter, closest, dist, edge, dist_max
if left_edge[0] <= self.center[0] <= right_edge[0] and \
left_edge[1] <= self.center[1] <= right_edge[1] and \
left_edge[2] <= self.center[2] <= right_edge[2]:
@@ -1737,7 +1740,9 @@
edge = right_edge[i] - left_edge[i]
closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0)
dist += closest * closest
- if dist <= self.mag[0]**2: return 1
+ dist_max = self.mag[0] * self.mag[0]
+ if dist <= dist_max:
+ return 1
return 0
def _hash_vals(self):
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/alt_ray_tracers.pyx
--- a/yt/utilities/lib/alt_ray_tracers.pyx
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -117,12 +117,12 @@
zleft = left_edges[:,1]
zright = right_edges[:,1]
- a = (dpcart[0]**2) + (dpcart[1]**2)
- b = (2*dpcart[0]*p1cart[0]) + (2*dpcart[1]*p1cart[1])
- cleft = ((p1cart[0]**2) + (p1cart[1]**2)) - rleft**2
- cright = ((p1cart[0]**2) + (p1cart[1]**2)) - rright**2
- twoa = 2*a
- bsqrd = b**2
+ a = dpcart[0] * dpcart[0] + dpcart[1] * dpcart[1]
+ b = 2 * dpcart[0] * p1cart[0] + 2 * dpcart[1] * p1cart[1]
+ cleft = p1cart[0] * p1cart[0] + p1cart[1] * p1cart[1] - rleft * rleft
+ cright = p1cart[0] * p1cart[0] + p1cart[1] * p1cart[1] - rright * rright
+ twoa = 2 * a
+ bsqrd = b * b
# Compute positive and negative times and associated masks
I = np.intp(left_edges.shape[0])
@@ -203,17 +203,17 @@
tsect, tinds = np.unique(tsect[tmask], return_index=True)
inds = inds[tmask][tinds]
xyz = dsect[tmask][tinds]
- s = np.sqrt(((xyz - p1cart)**2).sum(axis=1))
+ s = np.sqrt(((xyz - p1cart) * (xyz - p1cart)).sum(axis=1))
s, sinds = np.unique(s, return_index=True)
inds = inds[sinds]
xyz = xyz[sinds]
- t = s/np.sqrt((dpcart**2).sum())
+ t = s/np.sqrt((dpcart*dpcart).sum())
sinds = s.argsort()
s = s[sinds]
t = t[sinds]
inds = inds[sinds]
xyz = xyz[sinds]
- rztheta = np.concatenate([np.sqrt(xyz[:,0]**2 + xyz[:,1]**2)[:,np.newaxis],
+ rztheta = np.concatenate([np.sqrt(xyz[:,0] * xyz[:,0] + xyz[:,1] * xyz[:,1])[:,np.newaxis],
xyz[:,2:3],
np.arctan2(xyz[:,1], xyz[:,0])[:,np.newaxis]], axis=1)
return t, s, rztheta, inds
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/cosmology_time.pyx
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -1,12 +1,14 @@
cimport numpy as np
import numpy as np
-
+from libc.math cimport sqrt
cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0):
- return ( aexp_tau**3 * (O_mat_0 + O_vac_0*aexp_tau**3 + O_k_0*aexp_tau) )**0.5
+ double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
+ return sqrt( aexp_tau3 * (O_mat_0 + O_vac_0*aexp_tau3 + O_k_0*aexp_tau) )
cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
- return ( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t**3 + O_k_0*aexp_t) )**0.5
+ double aexp_t3 = aexp_t * aexp_t * aexp_t
+ return sqrt( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t3 + O_k_0*aexp_t) )
cdef step_cosmo(double alpha,double tau,double aexp_tau,double t,double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -445,18 +445,18 @@
+ rp*sm*tp*( r - s + t - 2.0)*vals[5] \
+ rp*sp*tp*( r + s + t - 2.0)*vals[6] \
+ rm*sp*tp*(-r + s + t - 2.0)*vals[7] \
- + 2.0*(1.0 - r**2)*sm*tm*vals[8] \
- + 2.0*rp*(1.0 - s**2)*tm*vals[9] \
- + 2.0*(1.0 - r**2)*sp*tm*vals[10] \
- + 2.0*rm*(1.0 - s**2)*tm*vals[11] \
- + 2.0*rm*sm*(1.0 - t**2)*vals[12] \
- + 2.0*rp*sm*(1.0 - t**2)*vals[13] \
- + 2.0*rp*sp*(1.0 - t**2)*vals[14] \
- + 2.0*rm*sp*(1.0 - t**2)*vals[15] \
- + 2.0*(1.0 - r**2)*sm*tp*vals[16] \
- + 2.0*rp*(1.0 - s**2)*tp*vals[17] \
- + 2.0*(1.0 - r**2)*sp*tp*vals[18] \
- + 2.0*rm*(1.0 - s**2)*tp*vals[19]
+ + 2.0*(1.0 - r*r)*sm*tm*vals[8] \
+ + 2.0*rp*(1.0 - s*s)*tm*vals[9] \
+ + 2.0*(1.0 - r*r)*sp*tm*vals[10] \
+ + 2.0*rm*(1.0 - s*s)*tm*vals[11] \
+ + 2.0*rm*sm*(1.0 - t*t)*vals[12] \
+ + 2.0*rp*sm*(1.0 - t*t)*vals[13] \
+ + 2.0*rp*sp*(1.0 - t*t)*vals[14] \
+ + 2.0*rm*sp*(1.0 - t*t)*vals[15] \
+ + 2.0*(1.0 - r*r)*sm*tp*vals[16] \
+ + 2.0*rp*(1.0 - s*s)*tp*vals[17] \
+ + 2.0*(1.0 - r*r)*sp*tp*vals[18] \
+ + 2.0*rm*(1.0 - s*s)*tp*vals[19]
return 0.125*F
@cython.boundscheck(False)
@@ -517,18 +517,18 @@
+ rp*sm*tp*( r - s + t - 2.0)*vertices[15 + i] \
+ rp*sp*tp*( r + s + t - 2.0)*vertices[18 + i] \
+ rm*sp*tp*(-r + s + t - 2.0)*vertices[21 + i] \
- + 2.0*(1.0 - r**2)*sm*tm*vertices[24 + i] \
- + 2.0*rp*(1.0 - s**2)*tm*vertices[27 + i] \
- + 2.0*(1.0 - r**2)*sp*tm*vertices[30 + i] \
- + 2.0*rm*(1.0 - s**2)*tm*vertices[33 + i] \
- + 2.0*rm*sm*(1.0 - t**2)*vertices[36 + i] \
- + 2.0*rp*sm*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*rp*sp*(1.0 - t**2)*vertices[42 + i] \
- + 2.0*rm*sp*(1.0 - t**2)*vertices[45 + i] \
- + 2.0*(1.0 - r**2)*sm*tp*vertices[48 + i] \
- + 2.0*rp*(1.0 - s**2)*tp*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*sp*tp*vertices[54 + i] \
- + 2.0*rm*(1.0 - s**2)*tp*vertices[57 + i] \
+ + 2.0*(1.0 - r*r)*sm*tm*vertices[24 + i] \
+ + 2.0*rp*(1.0 - s*s)*tm*vertices[27 + i] \
+ + 2.0*(1.0 - r*r)*sp*tm*vertices[30 + i] \
+ + 2.0*rm*(1.0 - s*s)*tm*vertices[33 + i] \
+ + 2.0*rm*sm*(1.0 - t*t)*vertices[36 + i] \
+ + 2.0*rp*sm*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*rp*sp*(1.0 - t*t)*vertices[42 + i] \
+ + 2.0*rm*sp*(1.0 - t*t)*vertices[45 + i] \
+ + 2.0*(1.0 - r*r)*sm*tp*vertices[48 + i] \
+ + 2.0*rp*(1.0 - s*s)*tp*vertices[51 + i] \
+ + 2.0*(1.0 - r*r)*sp*tp*vertices[54 + i] \
+ + 2.0*rm*(1.0 - s*s)*tp*vertices[57 + i] \
- 8.0*phys_x[i]
@@ -566,17 +566,17 @@
+ (sp*tp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (sp*tp*(r - s - t + 2.0) - rm*sp*tp)*vertices[21 + i] \
- 4.0*r*sm*tm*vertices[24 + i] \
- + 2.0*(1.0 - s**2)*tm*vertices[27 + i] \
+ + 2.0*(1.0 - s*s)*tm*vertices[27 + i] \
- 4.0*r*sp*tm*vertices[30 + i] \
- - 2.0*(1.0 - s**2)*tm*vertices[33 + i] \
- - 2.0*sm*(1.0 - t**2)*vertices[36 + i] \
- + 2.0*sm*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*sp*(1.0 - t**2)*vertices[42 + i] \
- - 2.0*sp*(1.0 - t**2)*vertices[45 + i] \
+ - 2.0*(1.0 - s*s)*tm*vertices[33 + i] \
+ - 2.0*sm*(1.0 - t*t)*vertices[36 + i] \
+ + 2.0*sm*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*sp*(1.0 - t*t)*vertices[42 + i] \
+ - 2.0*sp*(1.0 - t*t)*vertices[45 + i] \
- 4.0*r*sm*tp*vertices[48 + i] \
- + 2.0*(1.0 - s**2)*tp*vertices[51 + i] \
+ + 2.0*(1.0 - s*s)*tp*vertices[51 + i] \
- 4.0*r*sp*tp*vertices[54 + i] \
- - 2.0*(1.0 - s**2)*tp*vertices[57 + i]
+ - 2.0*(1.0 - s*s)*tp*vertices[57 + i]
scol[i] = ( rm*tm*(r + s + t + 2.0) - rm*sm*tm)*vertices[0 + i] \
+ (-rp*tm*(r - s - t - 2.0) - rp*sm*tm)*vertices[3 + i] \
+ ( rp*tm*(r + s - t - 2.0) + rp*sp*tm)*vertices[6 + i] \
@@ -585,17 +585,17 @@
+ (-rp*tp*(r - s + t - 2.0) - rp*sm*tp)*vertices[15 + i] \
+ ( rp*tp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (-rm*tp*(r - s - t + 2.0) + rm*sp*tp)*vertices[21 + i] \
- - 2.0*(1.0 - r**2)*tm*vertices[24 + i] \
+ - 2.0*(1.0 - r*r)*tm*vertices[24 + i] \
- 4.0*rp*s*tm*vertices[27 + i] \
- + 2.0*(1.0 - r**2)*tm*vertices[30 + i] \
+ + 2.0*(1.0 - r*r)*tm*vertices[30 + i] \
- 4.0*rm*s*tm*vertices[33 + i] \
- - 2.0*rm*(1.0 - t**2)*vertices[36 + i] \
- - 2.0*rp*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*rp*(1.0 - t**2)*vertices[42 + i] \
- + 2.0*rm*(1.0 - t**2)*vertices[45 + i] \
- - 2.0*(1.0 - r**2)*tp*vertices[48 + i] \
+ - 2.0*rm*(1.0 - t*t)*vertices[36 + i] \
+ - 2.0*rp*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*rp*(1.0 - t*t)*vertices[42 + i] \
+ + 2.0*rm*(1.0 - t*t)*vertices[45 + i] \
+ - 2.0*(1.0 - r*r)*tp*vertices[48 + i] \
- 4.0*rp*s*tp*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*tp*vertices[54 + i] \
+ + 2.0*(1.0 - r*r)*tp*vertices[54 + i] \
- 4.0*rm*s*tp*vertices[57 + i]
tcol[i] = ( rm*sm*(r + s + t + 2.0) - rm*sm*tm)*vertices[0 + i] \
+ (-rp*sm*(r - s - t - 2.0) - rp*sm*tm)*vertices[3 + i] \
@@ -605,18 +605,18 @@
+ ( rp*sm*(r - s + t - 2.0) + rp*sm*tp)*vertices[15 + i] \
+ ( rp*sp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (-rm*sp*(r - s - t + 2.0) + rm*sp*tp)*vertices[21 + i] \
- - 2.0*(1.0 - r**2)*sm*vertices[24 + i] \
- - 2.0*rp*(1.0 - s**2)*vertices[27 + i] \
- - 2.0*(1.0 - r**2)*sp*vertices[30 + i] \
- - 2.0*rm*(1.0 - s**2)*vertices[33 + i] \
+ - 2.0*(1.0 - r*r)*sm*vertices[24 + i] \
+ - 2.0*rp*(1.0 - s*s)*vertices[27 + i] \
+ - 2.0*(1.0 - r*r)*sp*vertices[30 + i] \
+ - 2.0*rm*(1.0 - s*s)*vertices[33 + i] \
- 4.0*rm*sm*t*vertices[36 + i] \
- 4.0*rp*sm*t*vertices[39 + i] \
- 4.0*rp*sp*t*vertices[42 + i] \
- 4.0*rm*sp*t*vertices[45 + i] \
- + 2.0*(1.0 - r**2)*sm*vertices[48 + i] \
- + 2.0*rp*(1.0 - s**2)*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*sp*vertices[54 + i] \
- + 2.0*rm*(1.0 - s**2)*vertices[57 + i]
+ + 2.0*(1.0 - r*r)*sm*vertices[48 + i] \
+ + 2.0*rp*(1.0 - s*s)*vertices[51 + i] \
+ + 2.0*(1.0 - r*r)*sp*vertices[54 + i] \
+ + 2.0*rm*(1.0 - s*s)*vertices[57 + i]
cdef class W1Sampler3D(NonlinearSolveSampler3D):
@@ -815,15 +815,19 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
- cdef double phi0, phi1, phi2, phi3, phi4, phi5
+ cdef double phi0, phi1, phi2, phi3, phi4, phi5, c0sq, c1sq, c0c1
- phi0 = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
- 2 * coord[1]**2 + 4 * coord[0] * coord[1]
- phi1 = -coord[0] + 2 * coord[0]**2
- phi2 = -coord[1] + 2 * coord[1]**2
- phi3 = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1]
- phi4 = 4 * coord[0] * coord[1]
- phi5 = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1]
+ c0sq = coord[0] * coord[0]
+ c1sq = coord[1] * coord[1]
+ c0c1 = coord[0] * coord[1]
+
+ phi0 = 1 - 3 * coord[0] + 2 * c0sq - 3 * coord[1] + \
+ 2 * c1sq + 4 * c0c1
+ phi1 = -coord[0] + 2 * c0sq
+ phi2 = -coord[1] + 2 * c1sq
+ phi3 = 4 * coord[0] - 4 * c0sq - 4 * c0c1
+ phi4 = 4 * c0c1
+ phi5 = 4 * coord[1] - 4 * c1sq - 4 * c0c1
return vals[0]*phi0 + vals[1]*phi1 + vals[2]*phi2 + vals[3]*phi3 + \
vals[4]*phi4 + vals[5]*phi5
@@ -861,22 +865,26 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double[10] phi
+ cdef double coordsq[3]
cdef int i
cdef double return_value = 0
- phi[0] = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
- 2 * coord[1]**2 - 3 * coord[2] + 2 * coord[2]**2 + \
+ for i in range(3):
+ coordsq[i] = coord[i] * coord[i]
+
+ phi[0] = 1 - 3 * coord[0] + 2 * coordsq[0] - 3 * coord[1] + \
+ 2 * coordsq[1] - 3 * coord[2] + 2 * coordsq[2] + \
4 * coord[0] * coord[1] + 4 * coord[0] * coord[2] + \
4 * coord[1] * coord[2]
- phi[1] = -coord[0] + 2 * coord[0]**2
- phi[2] = -coord[1] + 2 * coord[1]**2
- phi[3] = -coord[2] + 2 * coord[2]**2
- phi[4] = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1] - \
+ phi[1] = -coord[0] + 2 * coordsq[0]
+ phi[2] = -coord[1] + 2 * coordsq[1]
+ phi[3] = -coord[2] + 2 * coordsq[2]
+ phi[4] = 4 * coord[0] - 4 * coordsq[0] - 4 * coord[0] * coord[1] - \
4 * coord[0] * coord[2]
phi[5] = 4 * coord[0] * coord[1]
- phi[6] = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1] - \
+ phi[6] = 4 * coord[1] - 4 * coordsq[1] - 4 * coord[0] * coord[1] - \
4 * coord[1] * coord[2]
- phi[7] = 4 * coord[2] - 4 * coord[2]**2 - 4 * coord[2] * coord[0] - \
+ phi[7] = 4 * coord[2] - 4 * coordsq[2] - 4 * coord[2] * coord[0] - \
4 * coord[2] * coord[1]
phi[8] = 4 * coord[0] * coord[2]
phi[9] = 4 * coord[1] * coord[2]
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -19,7 +19,7 @@
#cimport healpix_interface
from libc.stdlib cimport malloc, calloc, free, abs
from libc.math cimport exp, floor, log2, \
- fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
+ fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI, sqrt
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, i64clip
from field_interpolation_tables cimport \
FieldInterpolationTable, FIT_initialize_table, FIT_eval_transfer,\
@@ -317,7 +317,7 @@
# zb = (x*x/8.0 + y*y/2.0 - 1.0)
# if zb > 0: continue
# z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- # z = z**0.5
+ # z = sqrt(z)
# # Longitude
# phi = (2.0*atan(z*x/(2.0 * (2.0*z*z-1.0))) + pi)
# # Latitude
@@ -352,7 +352,7 @@
px = (2.0 * (nimi*nx + i)) / resolution - 1.0
for j in range(ny):
py = (2.0 * (nimj*ny + j)) / resolution - 1.0
- r = (px*px + py*py)**0.5
+ r = sqrt(px*px + py*py)
if r > 1.01:
vp[i,j,0] = vp[i,j,1] = vp[i,j,2] = 0.0
continue
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/image_samplers.pyx
--- a/yt/utilities/lib/image_samplers.pyx
+++ b/yt/utilities/lib/image_samplers.pyx
@@ -412,7 +412,9 @@
self.vra.n_samples = n_samples
self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
- light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
+ light_dir /= np.sqrt(light_dir[0] * light_dir[0] +
+ light_dir[1] * light_dir[1] +
+ light_dir[2] * light_dir[2])
for i in range(3):
self.vra.light_dir[i] = light_dir[i]
for i in range(4):
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/lenses.pyx
--- a/yt/utilities/lib/lenses.pyx
+++ b/yt/utilities/lib/lenses.pyx
@@ -140,8 +140,9 @@
if acos(sight_angle_cos) < 0.5 * M_PI and sight_angle_cos != 0.0:
sight_length = cam_width[2] / sight_angle_cos
else:
- sight_length = sqrt(cam_width[0]**2 + cam_width[1]**2)
- sight_length = sight_length / sqrt(1.0 - sight_angle_cos**2)
+ sight_length = sqrt(cam_width[0] * cam_width[0] +
+ cam_width[1] * cam_width[1])
+ sight_length /= sqrt(1.0 - sight_angle_cos * sight_angle_cos)
fma(sight_length, sight_vector, cam_pos, pos1)
subtract(pos1, sight_center, pos1)
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/marching_cubes.pyx
--- a/yt/utilities/lib/marching_cubes.pyx
+++ b/yt/utilities/lib/marching_cubes.pyx
@@ -18,6 +18,7 @@
import numpy as np
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip
from libc.stdlib cimport malloc, free, abs
+from libc.math cimport sqrt
from fixed_interpolator cimport *
cdef extern from "marching_cubes.h":
@@ -357,7 +358,7 @@
for n in range(3):
temp += normal[n]*normal[n]
# Take the negative, to ensure it points inwardly
- temp = -(temp**0.5)
+ temp = -sqrt(temp)
# Dump this somewhere for now
temp = wval * (fv[0] * normal[0] +
fv[1] * normal[1] +
@@ -368,15 +369,15 @@
for n in range(3):
fv[n] = 0.0
for n in range(3):
- fv[0] += (current.p[0][n] - current.p[2][n])**2.0
- fv[1] += (current.p[1][n] - current.p[0][n])**2.0
- fv[2] += (current.p[2][n] - current.p[1][n])**2.0
+ fv[0] += (current.p[0][n] - current.p[2][n]) * (current.p[0][n] - current.p[2][n])
+ fv[1] += (current.p[1][n] - current.p[0][n]) * (current.p[1][n] - current.p[0][n])
+ fv[2] += (current.p[2][n] - current.p[1][n]) * (current.p[2][n] - current.p[1][n])
s = 0.0
for n in range(3):
- fv[n] = fv[n]**0.5
+ fv[n] = sqrt(fv[n])
s += 0.5 * fv[n]
area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))
- area = area**0.5
+ area = sqrt(area)
flux += temp*area
last = current
if current.next == NULL: break
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -47,7 +47,7 @@
np.ndarray[np.float64_t, ndim=2] qresult,
np.ndarray[np.uint8_t, ndim=1, cast=True] used):
cdef int n, fi, bin
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -57,12 +57,13 @@
wresult[bin] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin,fi]
# qresult has to have the previous wresult
- qresult[bin,fi] += (oldwr * wval * (bval - mresult[bin,fi])**2) / \
+ qresult[bin,fi] += (oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin,fi] += wval * (bval - mresult[bin,fi]) / wresult[bin]
+ mresult[bin,fi] += wval * bval_mresult / wresult[bin]
used[bin] = 1
return
@@ -79,7 +80,7 @@
np.ndarray[np.float64_t, ndim=3] qresult,
np.ndarray[np.uint8_t, ndim=2, cast=True] used):
cdef int n, fi, bin_x, bin_y
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -90,12 +91,13 @@
wresult[bin_x,bin_y] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin_x,bin_y,fi]
# qresult has to have the previous wresult
- qresult[bin_x,bin_y,fi] += (oldwr * wval * (bval - mresult[bin_x,bin_y,fi])**2) / \
+ qresult[bin_x,bin_y,fi] += oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin_x,bin_y,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin_x,bin_y,fi] += wval * (bval - mresult[bin_x,bin_y,fi]) / wresult[bin_x,bin_y]
+ mresult[bin_x,bin_y,fi] += wval * bval_mresult / wresult[bin_x,bin_y]
used[bin_x,bin_y] = 1
return
@@ -113,7 +115,7 @@
np.ndarray[np.float64_t, ndim=4] qresult,
np.ndarray[np.uint8_t, ndim=3, cast=True] used):
cdef int n, fi, bin_x, bin_y, bin_z
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -125,14 +127,14 @@
wresult[bin_x,bin_y,bin_z] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin_x,bin_y,bin_z,fi]
# qresult has to have the previous wresult
qresult[bin_x,bin_y,bin_z,fi] += \
- (oldwr * wval * (bval - mresult[bin_x,bin_y,bin_z,fi])**2) / \
+ oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin_x,bin_y,bin_z,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin_x,bin_y,bin_z,fi] += wval * \
- (bval - mresult[bin_x,bin_y,bin_z,fi]) / \
+ mresult[bin_x,bin_y,bin_z,fi] += wval * bval_mresult / \
wresult[bin_x,bin_y,bin_z]
used[bin_x,bin_y,bin_z] = 1
return
@@ -149,16 +151,17 @@
np.ndarray[np.float64_t, ndim=1] qresult,
np.ndarray[np.float64_t, ndim=1] used):
cdef int n, bin
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bin = bins_x[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bin] += (wresult[bin] * wval * (bval - mresult[bin])**2) / \
+ bval_mresult = bval - mresult[bin]
+ qresult[bin] += wresult[bin] * wval * bval_mresult * bval_mresult / \
(wresult[bin] + wval)
wresult[bin] += wval
bresult[bin] += wval*bval
- mresult[bin] += wval * (bval - mresult[bin]) / wresult[bin]
+ mresult[bin] += wval * bval_mresult / wresult[bin]
used[bin] = 1
return
@@ -175,17 +178,18 @@
np.ndarray[np.float64_t, ndim=2] qresult,
np.ndarray[np.float64_t, ndim=2] used):
cdef int n, bini, binj
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bini = bins_x[n]
binj = bins_y[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bini, binj] += (wresult[bini, binj] * wval * (bval - mresult[bini, binj])**2) / \
+ bval_mresult = bval - mresult[bini, binj]
+ qresult[bini, binj] += (wresult[bini, binj] * wval * bval_mresult * bval_mresult / \
(wresult[bini, binj] + wval)
wresult[bini, binj] += wval
bresult[bini, binj] += wval*bval
- mresult[bini, binj] += wval * (bval - mresult[bini, binj]) / wresult[bini, binj]
+ mresult[bini, binj] += wval * bval_mresult / wresult[bini, binj]
used[bini, binj] = 1
return
@@ -203,18 +207,19 @@
np.ndarray[np.float64_t, ndim=3] qresult,
np.ndarray[np.float64_t, ndim=3] used):
cdef int n, bini, binj, bink
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bini = bins_x[n]
binj = bins_y[n]
bink = bins_z[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bini, binj, bink] += (wresult[bini, binj, bink] * wval * (bval - mresult[bini, binj, bink])**2) / \
+ bval_mresult = bval - mresult[bini, binj, bink]
+ qresult[bini, binj, bink] += wresult[bini, binj, bink] * wval * bval_mresult * bval_mresult / \
(wresult[bini, binj, bink] + wval)
wresult[bini, binj, bink] += wval
bresult[bini, binj, bink] += wval*bval
- mresult[bini, binj, bink] += wval * (bval - mresult[bini, binj, bink]) / wresult[bini, binj, bink]
+ mresult[bini, binj, bink] += wval * bval_mresult / wresult[bini, binj, bink]
used[bini, binj, bink] = 1
return
@@ -344,8 +349,8 @@
z1 = zs[j+1]
dx = abs(x1-x0)
dy = abs(y1-y0)
- dzx = (z1-z0) / (dx**2 + dy**2) * dx
- dzy = (z1-z0) / (dx**2 + dy**2) * dy
+ dzx = (z1-z0) / (dx * dx + dy * dy) * dx
+ dzy = (z1-z0) / (dx * dx + dy * dy) * dy
err = dx - dy
if crop == 1 and (dx > nx/2.0 or dy > ny/2.0):
continue
@@ -1016,7 +1021,7 @@
i = 0
n_q = mass.size
pbar = get_pbar("Calculating potential for %d cells" % n_q,
- 0.5 * (n_q**2 - n_q))
+ 0.5 * (n_q * n_q - n_q))
for q_outer in range(n_q - 1):
this_potential = 0.
mass_o = mass[q_outer]
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -352,8 +352,8 @@
for i in range(8):
rbounds[0] = fmin(rbounds[0], corners[i])
rbounds[1] = fmax(rbounds[1], corners[i])
- rbounds[0] = rbounds[0]**0.5
- rbounds[1] = rbounds[1]**0.5
+ rbounds[0] = math.sqrt(rbounds[0])
+ rbounds[1] = math.sqrt(rbounds[1])
# If we include the origin in either direction, we need to have radius of
# zero as our lower bound.
if x0 < 0 and x1 > 0:
@@ -479,8 +479,8 @@
y = (-1.0 + j * dy)*s2
zb = (x*x/8.0 + y*y/2.0 - 1.0)
if zb > 0: continue
- z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- z = z**0.5
+ z = (1.0 - (x * 0.25) * (x * 0.25) - (y * 0.5) * (y * 0.5))
+ z = math.sqrt(z)
# Longitude
theta0 = 2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0)))
# Latitude
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/lib/points_in_volume.pyx
--- a/yt/utilities/lib/points_in_volume.pyx
+++ b/yt/utilities/lib/points_in_volume.pyx
@@ -17,6 +17,7 @@
import numpy as np
cimport numpy as np
cimport cython
+from libc.math cimport sqrt
cdef extern from "math.h":
double fabs(double x)
@@ -123,7 +124,7 @@
cdef np.float64_t norm = 0.0
for i in range(3):
norm += vec[i]*vec[i]
- norm = norm**0.5
+ norm = sqrt(norm)
for i in range(3):
vec[i] /= norm
diff -r bb7f23a27456bed189cbd20b0b5cdc968a2f31e0 -r b850dff755db3483deebebe018e1d8639fc1f8e9 yt/utilities/spatial/ckdtree.pyx
--- a/yt/utilities/spatial/ckdtree.pyx
+++ b/yt/utilities/spatial/ckdtree.pyx
@@ -3,6 +3,7 @@
import numpy as np
cimport numpy as np
cimport libc.stdlib as stdlib
+from libc.math cimport sqrt
cimport cython
import kdtree
@@ -700,7 +701,7 @@
for j in range(num_neighbors):
pj = tags_temp[j]
r2 = dist_temp[j] * ih2
- rs = 2.0 - (r2**0.5)
+ rs = 2.0 - sqrt(r2)
if (r2 < 1.0):
rs = (1.0 - 0.75*rs*r2)
else:
https://bitbucket.org/yt_analysis/yt/commits/f17df4b963f9/
Changeset: f17df4b963f9
Branch: yt
User: xarthisius
Date: 2016-12-02 22:20:13+00:00
Summary: Fix syntax error
Affected #: 1 file
diff -r b850dff755db3483deebebe018e1d8639fc1f8e9 -r f17df4b963f97baf40f73b7b246207c431ee25fd yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -59,7 +59,7 @@
bval = bsource[n,fi]
bval_mresult = bval - mresult[bin,fi]
# qresult has to have the previous wresult
- qresult[bin,fi] += (oldwr * wval * bval_mresult * bval_mresult / \
+ qresult[bin,fi] += oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin,fi] += wval*bval
# mresult needs the new wresult
@@ -185,7 +185,7 @@
bval = bsource[n]
wval = wsource[n]
bval_mresult = bval - mresult[bini, binj]
- qresult[bini, binj] += (wresult[bini, binj] * wval * bval_mresult * bval_mresult / \
+ qresult[bini, binj] += wresult[bini, binj] * wval * bval_mresult * bval_mresult / \
(wresult[bini, binj] + wval)
wresult[bini, binj] += wval
bresult[bini, binj] += wval*bval
https://bitbucket.org/yt_analysis/yt/commits/d13b12c9b1ad/
Changeset: d13b12c9b1ad
Branch: yt
User: xarthisius
Date: 2016-12-03 00:25:47+00:00
Summary: Add missing cdef statements
Affected #: 1 file
diff -r f17df4b963f97baf40f73b7b246207c431ee25fd -r d13b12c9b1ad7a31a8d5727f9f67e04f617e2742 yt/utilities/lib/cosmology_time.pyx
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -3,11 +3,11 @@
from libc.math cimport sqrt
cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0):
- double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
+ cdef double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
return sqrt( aexp_tau3 * (O_mat_0 + O_vac_0*aexp_tau3 + O_k_0*aexp_tau) )
cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
- double aexp_t3 = aexp_t * aexp_t * aexp_t
+ cdef double aexp_t3 = aexp_t * aexp_t * aexp_t
return sqrt( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t3 + O_k_0*aexp_t) )
https://bitbucket.org/yt_analysis/yt/commits/0874610be596/
Changeset: 0874610be596
Branch: yt
User: ngoldbaum
Date: 2016-12-05 19:02:07+00:00
Summary: Merged in xarthisius/yt (pull request #2454)
Use explicit multiplication or sqrt instead of pow(..., 2/0.5)
Affected #: 13 files
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -311,7 +311,7 @@
else:
self.mk[ii[2], ii[1], ii[0], offset] = mk + (fields[0] - mk) / k
self.qk[ii[2], ii[1], ii[0], offset] = \
- qk + (k - 1.0) * (fields[0] - mk)**2.0 / k
+ qk + (k - 1.0) * (fields[0] - mk) * (fields[0] - mk) / k
self.i[ii[2], ii[1], ii[0], offset] += 1
def finalize(self):
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -1711,10 +1711,13 @@
cdef int select_sphere(self, np.float64_t pos[3], np.float64_t radius) nogil:
# this is the sphere selection
cdef int i
- cdef np.float64_t dist2 = 0
+ cdef np.float64_t dist, dist2_max, dist2 = 0
for i in range(3):
- dist2 += self.difference(pos[i], self.center[i], i)**2
- if dist2 <= (self.mag[0]+radius)**2: return 1
+ dist = self.difference(pos[i], self.center[i], i)
+ dist2 += dist * dist
+ dist2_max = (self.mag[0] + radius) * (self.mag[0] + radius)
+ if dist2 <= dist2_max:
+ return 1
return 0
@cython.boundscheck(False)
@@ -1724,7 +1727,7 @@
np.float64_t right_edge[3]) nogil:
# This is the sphere selection
cdef int i
- cdef np.float64_t box_center, relcenter, closest, dist, edge
+ cdef np.float64_t box_center, relcenter, closest, dist, edge, dist_max
if left_edge[0] <= self.center[0] <= right_edge[0] and \
left_edge[1] <= self.center[1] <= right_edge[1] and \
left_edge[2] <= self.center[2] <= right_edge[2]:
@@ -1737,7 +1740,9 @@
edge = right_edge[i] - left_edge[i]
closest = relcenter - fclip(relcenter, -edge/2.0, edge/2.0)
dist += closest * closest
- if dist <= self.mag[0]**2: return 1
+ dist_max = self.mag[0] * self.mag[0]
+ if dist <= dist_max:
+ return 1
return 0
def _hash_vals(self):
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/alt_ray_tracers.pyx
--- a/yt/utilities/lib/alt_ray_tracers.pyx
+++ b/yt/utilities/lib/alt_ray_tracers.pyx
@@ -117,12 +117,12 @@
zleft = left_edges[:,1]
zright = right_edges[:,1]
- a = (dpcart[0]**2) + (dpcart[1]**2)
- b = (2*dpcart[0]*p1cart[0]) + (2*dpcart[1]*p1cart[1])
- cleft = ((p1cart[0]**2) + (p1cart[1]**2)) - rleft**2
- cright = ((p1cart[0]**2) + (p1cart[1]**2)) - rright**2
- twoa = 2*a
- bsqrd = b**2
+ a = dpcart[0] * dpcart[0] + dpcart[1] * dpcart[1]
+ b = 2 * dpcart[0] * p1cart[0] + 2 * dpcart[1] * p1cart[1]
+ cleft = p1cart[0] * p1cart[0] + p1cart[1] * p1cart[1] - rleft * rleft
+ cright = p1cart[0] * p1cart[0] + p1cart[1] * p1cart[1] - rright * rright
+ twoa = 2 * a
+ bsqrd = b * b
# Compute positive and negative times and associated masks
I = np.intp(left_edges.shape[0])
@@ -203,17 +203,17 @@
tsect, tinds = np.unique(tsect[tmask], return_index=True)
inds = inds[tmask][tinds]
xyz = dsect[tmask][tinds]
- s = np.sqrt(((xyz - p1cart)**2).sum(axis=1))
+ s = np.sqrt(((xyz - p1cart) * (xyz - p1cart)).sum(axis=1))
s, sinds = np.unique(s, return_index=True)
inds = inds[sinds]
xyz = xyz[sinds]
- t = s/np.sqrt((dpcart**2).sum())
+ t = s/np.sqrt((dpcart*dpcart).sum())
sinds = s.argsort()
s = s[sinds]
t = t[sinds]
inds = inds[sinds]
xyz = xyz[sinds]
- rztheta = np.concatenate([np.sqrt(xyz[:,0]**2 + xyz[:,1]**2)[:,np.newaxis],
+ rztheta = np.concatenate([np.sqrt(xyz[:,0] * xyz[:,0] + xyz[:,1] * xyz[:,1])[:,np.newaxis],
xyz[:,2:3],
np.arctan2(xyz[:,1], xyz[:,0])[:,np.newaxis]], axis=1)
return t, s, rztheta, inds
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/cosmology_time.pyx
--- a/yt/utilities/lib/cosmology_time.pyx
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -1,12 +1,14 @@
cimport numpy as np
import numpy as np
-
+from libc.math cimport sqrt
cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0):
- return ( aexp_tau**3 * (O_mat_0 + O_vac_0*aexp_tau**3 + O_k_0*aexp_tau) )**0.5
+ cdef double aexp_tau3 = aexp_tau * aexp_tau * aexp_tau
+ return sqrt( aexp_tau3 * (O_mat_0 + O_vac_0*aexp_tau3 + O_k_0*aexp_tau) )
cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
- return ( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t**3 + O_k_0*aexp_t) )**0.5
+ cdef double aexp_t3 = aexp_t * aexp_t * aexp_t
+ return sqrt( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t3 + O_k_0*aexp_t) )
cdef step_cosmo(double alpha,double tau,double aexp_tau,double t,double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -445,18 +445,18 @@
+ rp*sm*tp*( r - s + t - 2.0)*vals[5] \
+ rp*sp*tp*( r + s + t - 2.0)*vals[6] \
+ rm*sp*tp*(-r + s + t - 2.0)*vals[7] \
- + 2.0*(1.0 - r**2)*sm*tm*vals[8] \
- + 2.0*rp*(1.0 - s**2)*tm*vals[9] \
- + 2.0*(1.0 - r**2)*sp*tm*vals[10] \
- + 2.0*rm*(1.0 - s**2)*tm*vals[11] \
- + 2.0*rm*sm*(1.0 - t**2)*vals[12] \
- + 2.0*rp*sm*(1.0 - t**2)*vals[13] \
- + 2.0*rp*sp*(1.0 - t**2)*vals[14] \
- + 2.0*rm*sp*(1.0 - t**2)*vals[15] \
- + 2.0*(1.0 - r**2)*sm*tp*vals[16] \
- + 2.0*rp*(1.0 - s**2)*tp*vals[17] \
- + 2.0*(1.0 - r**2)*sp*tp*vals[18] \
- + 2.0*rm*(1.0 - s**2)*tp*vals[19]
+ + 2.0*(1.0 - r*r)*sm*tm*vals[8] \
+ + 2.0*rp*(1.0 - s*s)*tm*vals[9] \
+ + 2.0*(1.0 - r*r)*sp*tm*vals[10] \
+ + 2.0*rm*(1.0 - s*s)*tm*vals[11] \
+ + 2.0*rm*sm*(1.0 - t*t)*vals[12] \
+ + 2.0*rp*sm*(1.0 - t*t)*vals[13] \
+ + 2.0*rp*sp*(1.0 - t*t)*vals[14] \
+ + 2.0*rm*sp*(1.0 - t*t)*vals[15] \
+ + 2.0*(1.0 - r*r)*sm*tp*vals[16] \
+ + 2.0*rp*(1.0 - s*s)*tp*vals[17] \
+ + 2.0*(1.0 - r*r)*sp*tp*vals[18] \
+ + 2.0*rm*(1.0 - s*s)*tp*vals[19]
return 0.125*F
@cython.boundscheck(False)
@@ -517,18 +517,18 @@
+ rp*sm*tp*( r - s + t - 2.0)*vertices[15 + i] \
+ rp*sp*tp*( r + s + t - 2.0)*vertices[18 + i] \
+ rm*sp*tp*(-r + s + t - 2.0)*vertices[21 + i] \
- + 2.0*(1.0 - r**2)*sm*tm*vertices[24 + i] \
- + 2.0*rp*(1.0 - s**2)*tm*vertices[27 + i] \
- + 2.0*(1.0 - r**2)*sp*tm*vertices[30 + i] \
- + 2.0*rm*(1.0 - s**2)*tm*vertices[33 + i] \
- + 2.0*rm*sm*(1.0 - t**2)*vertices[36 + i] \
- + 2.0*rp*sm*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*rp*sp*(1.0 - t**2)*vertices[42 + i] \
- + 2.0*rm*sp*(1.0 - t**2)*vertices[45 + i] \
- + 2.0*(1.0 - r**2)*sm*tp*vertices[48 + i] \
- + 2.0*rp*(1.0 - s**2)*tp*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*sp*tp*vertices[54 + i] \
- + 2.0*rm*(1.0 - s**2)*tp*vertices[57 + i] \
+ + 2.0*(1.0 - r*r)*sm*tm*vertices[24 + i] \
+ + 2.0*rp*(1.0 - s*s)*tm*vertices[27 + i] \
+ + 2.0*(1.0 - r*r)*sp*tm*vertices[30 + i] \
+ + 2.0*rm*(1.0 - s*s)*tm*vertices[33 + i] \
+ + 2.0*rm*sm*(1.0 - t*t)*vertices[36 + i] \
+ + 2.0*rp*sm*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*rp*sp*(1.0 - t*t)*vertices[42 + i] \
+ + 2.0*rm*sp*(1.0 - t*t)*vertices[45 + i] \
+ + 2.0*(1.0 - r*r)*sm*tp*vertices[48 + i] \
+ + 2.0*rp*(1.0 - s*s)*tp*vertices[51 + i] \
+ + 2.0*(1.0 - r*r)*sp*tp*vertices[54 + i] \
+ + 2.0*rm*(1.0 - s*s)*tp*vertices[57 + i] \
- 8.0*phys_x[i]
@@ -566,17 +566,17 @@
+ (sp*tp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (sp*tp*(r - s - t + 2.0) - rm*sp*tp)*vertices[21 + i] \
- 4.0*r*sm*tm*vertices[24 + i] \
- + 2.0*(1.0 - s**2)*tm*vertices[27 + i] \
+ + 2.0*(1.0 - s*s)*tm*vertices[27 + i] \
- 4.0*r*sp*tm*vertices[30 + i] \
- - 2.0*(1.0 - s**2)*tm*vertices[33 + i] \
- - 2.0*sm*(1.0 - t**2)*vertices[36 + i] \
- + 2.0*sm*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*sp*(1.0 - t**2)*vertices[42 + i] \
- - 2.0*sp*(1.0 - t**2)*vertices[45 + i] \
+ - 2.0*(1.0 - s*s)*tm*vertices[33 + i] \
+ - 2.0*sm*(1.0 - t*t)*vertices[36 + i] \
+ + 2.0*sm*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*sp*(1.0 - t*t)*vertices[42 + i] \
+ - 2.0*sp*(1.0 - t*t)*vertices[45 + i] \
- 4.0*r*sm*tp*vertices[48 + i] \
- + 2.0*(1.0 - s**2)*tp*vertices[51 + i] \
+ + 2.0*(1.0 - s*s)*tp*vertices[51 + i] \
- 4.0*r*sp*tp*vertices[54 + i] \
- - 2.0*(1.0 - s**2)*tp*vertices[57 + i]
+ - 2.0*(1.0 - s*s)*tp*vertices[57 + i]
scol[i] = ( rm*tm*(r + s + t + 2.0) - rm*sm*tm)*vertices[0 + i] \
+ (-rp*tm*(r - s - t - 2.0) - rp*sm*tm)*vertices[3 + i] \
+ ( rp*tm*(r + s - t - 2.0) + rp*sp*tm)*vertices[6 + i] \
@@ -585,17 +585,17 @@
+ (-rp*tp*(r - s + t - 2.0) - rp*sm*tp)*vertices[15 + i] \
+ ( rp*tp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (-rm*tp*(r - s - t + 2.0) + rm*sp*tp)*vertices[21 + i] \
- - 2.0*(1.0 - r**2)*tm*vertices[24 + i] \
+ - 2.0*(1.0 - r*r)*tm*vertices[24 + i] \
- 4.0*rp*s*tm*vertices[27 + i] \
- + 2.0*(1.0 - r**2)*tm*vertices[30 + i] \
+ + 2.0*(1.0 - r*r)*tm*vertices[30 + i] \
- 4.0*rm*s*tm*vertices[33 + i] \
- - 2.0*rm*(1.0 - t**2)*vertices[36 + i] \
- - 2.0*rp*(1.0 - t**2)*vertices[39 + i] \
- + 2.0*rp*(1.0 - t**2)*vertices[42 + i] \
- + 2.0*rm*(1.0 - t**2)*vertices[45 + i] \
- - 2.0*(1.0 - r**2)*tp*vertices[48 + i] \
+ - 2.0*rm*(1.0 - t*t)*vertices[36 + i] \
+ - 2.0*rp*(1.0 - t*t)*vertices[39 + i] \
+ + 2.0*rp*(1.0 - t*t)*vertices[42 + i] \
+ + 2.0*rm*(1.0 - t*t)*vertices[45 + i] \
+ - 2.0*(1.0 - r*r)*tp*vertices[48 + i] \
- 4.0*rp*s*tp*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*tp*vertices[54 + i] \
+ + 2.0*(1.0 - r*r)*tp*vertices[54 + i] \
- 4.0*rm*s*tp*vertices[57 + i]
tcol[i] = ( rm*sm*(r + s + t + 2.0) - rm*sm*tm)*vertices[0 + i] \
+ (-rp*sm*(r - s - t - 2.0) - rp*sm*tm)*vertices[3 + i] \
@@ -605,18 +605,18 @@
+ ( rp*sm*(r - s + t - 2.0) + rp*sm*tp)*vertices[15 + i] \
+ ( rp*sp*(r + s + t - 2.0) + rp*sp*tp)*vertices[18 + i] \
+ (-rm*sp*(r - s - t + 2.0) + rm*sp*tp)*vertices[21 + i] \
- - 2.0*(1.0 - r**2)*sm*vertices[24 + i] \
- - 2.0*rp*(1.0 - s**2)*vertices[27 + i] \
- - 2.0*(1.0 - r**2)*sp*vertices[30 + i] \
- - 2.0*rm*(1.0 - s**2)*vertices[33 + i] \
+ - 2.0*(1.0 - r*r)*sm*vertices[24 + i] \
+ - 2.0*rp*(1.0 - s*s)*vertices[27 + i] \
+ - 2.0*(1.0 - r*r)*sp*vertices[30 + i] \
+ - 2.0*rm*(1.0 - s*s)*vertices[33 + i] \
- 4.0*rm*sm*t*vertices[36 + i] \
- 4.0*rp*sm*t*vertices[39 + i] \
- 4.0*rp*sp*t*vertices[42 + i] \
- 4.0*rm*sp*t*vertices[45 + i] \
- + 2.0*(1.0 - r**2)*sm*vertices[48 + i] \
- + 2.0*rp*(1.0 - s**2)*vertices[51 + i] \
- + 2.0*(1.0 - r**2)*sp*vertices[54 + i] \
- + 2.0*rm*(1.0 - s**2)*vertices[57 + i]
+ + 2.0*(1.0 - r*r)*sm*vertices[48 + i] \
+ + 2.0*rp*(1.0 - s*s)*vertices[51 + i] \
+ + 2.0*(1.0 - r*r)*sp*vertices[54 + i] \
+ + 2.0*rm*(1.0 - s*s)*vertices[57 + i]
cdef class W1Sampler3D(NonlinearSolveSampler3D):
@@ -815,15 +815,19 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
- cdef double phi0, phi1, phi2, phi3, phi4, phi5
+ cdef double phi0, phi1, phi2, phi3, phi4, phi5, c0sq, c1sq, c0c1
- phi0 = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
- 2 * coord[1]**2 + 4 * coord[0] * coord[1]
- phi1 = -coord[0] + 2 * coord[0]**2
- phi2 = -coord[1] + 2 * coord[1]**2
- phi3 = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1]
- phi4 = 4 * coord[0] * coord[1]
- phi5 = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1]
+ c0sq = coord[0] * coord[0]
+ c1sq = coord[1] * coord[1]
+ c0c1 = coord[0] * coord[1]
+
+ phi0 = 1 - 3 * coord[0] + 2 * c0sq - 3 * coord[1] + \
+ 2 * c1sq + 4 * c0c1
+ phi1 = -coord[0] + 2 * c0sq
+ phi2 = -coord[1] + 2 * c1sq
+ phi3 = 4 * coord[0] - 4 * c0sq - 4 * c0c1
+ phi4 = 4 * c0c1
+ phi5 = 4 * coord[1] - 4 * c1sq - 4 * c0c1
return vals[0]*phi0 + vals[1]*phi1 + vals[2]*phi2 + vals[3]*phi3 + \
vals[4]*phi4 + vals[5]*phi5
@@ -861,22 +865,26 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double[10] phi
+ cdef double coordsq[3]
cdef int i
cdef double return_value = 0
- phi[0] = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
- 2 * coord[1]**2 - 3 * coord[2] + 2 * coord[2]**2 + \
+ for i in range(3):
+ coordsq[i] = coord[i] * coord[i]
+
+ phi[0] = 1 - 3 * coord[0] + 2 * coordsq[0] - 3 * coord[1] + \
+ 2 * coordsq[1] - 3 * coord[2] + 2 * coordsq[2] + \
4 * coord[0] * coord[1] + 4 * coord[0] * coord[2] + \
4 * coord[1] * coord[2]
- phi[1] = -coord[0] + 2 * coord[0]**2
- phi[2] = -coord[1] + 2 * coord[1]**2
- phi[3] = -coord[2] + 2 * coord[2]**2
- phi[4] = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1] - \
+ phi[1] = -coord[0] + 2 * coordsq[0]
+ phi[2] = -coord[1] + 2 * coordsq[1]
+ phi[3] = -coord[2] + 2 * coordsq[2]
+ phi[4] = 4 * coord[0] - 4 * coordsq[0] - 4 * coord[0] * coord[1] - \
4 * coord[0] * coord[2]
phi[5] = 4 * coord[0] * coord[1]
- phi[6] = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1] - \
+ phi[6] = 4 * coord[1] - 4 * coordsq[1] - 4 * coord[0] * coord[1] - \
4 * coord[1] * coord[2]
- phi[7] = 4 * coord[2] - 4 * coord[2]**2 - 4 * coord[2] * coord[0] - \
+ phi[7] = 4 * coord[2] - 4 * coordsq[2] - 4 * coord[2] * coord[0] - \
4 * coord[2] * coord[1]
phi[8] = 4 * coord[0] * coord[2]
phi[9] = 4 * coord[1] * coord[2]
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -19,7 +19,7 @@
#cimport healpix_interface
from libc.stdlib cimport malloc, calloc, free, abs
from libc.math cimport exp, floor, log2, \
- fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
+ fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI, sqrt
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, i64clip
from field_interpolation_tables cimport \
FieldInterpolationTable, FIT_initialize_table, FIT_eval_transfer,\
@@ -317,7 +317,7 @@
# zb = (x*x/8.0 + y*y/2.0 - 1.0)
# if zb > 0: continue
# z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- # z = z**0.5
+ # z = sqrt(z)
# # Longitude
# phi = (2.0*atan(z*x/(2.0 * (2.0*z*z-1.0))) + pi)
# # Latitude
@@ -352,7 +352,7 @@
px = (2.0 * (nimi*nx + i)) / resolution - 1.0
for j in range(ny):
py = (2.0 * (nimj*ny + j)) / resolution - 1.0
- r = (px*px + py*py)**0.5
+ r = sqrt(px*px + py*py)
if r > 1.01:
vp[i,j,0] = vp[i,j,1] = vp[i,j,2] = 0.0
continue
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/image_samplers.pyx
--- a/yt/utilities/lib/image_samplers.pyx
+++ b/yt/utilities/lib/image_samplers.pyx
@@ -412,7 +412,9 @@
self.vra.n_samples = n_samples
self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
- light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
+ light_dir /= np.sqrt(light_dir[0] * light_dir[0] +
+ light_dir[1] * light_dir[1] +
+ light_dir[2] * light_dir[2])
for i in range(3):
self.vra.light_dir[i] = light_dir[i]
for i in range(4):
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/lenses.pyx
--- a/yt/utilities/lib/lenses.pyx
+++ b/yt/utilities/lib/lenses.pyx
@@ -140,8 +140,9 @@
if acos(sight_angle_cos) < 0.5 * M_PI and sight_angle_cos != 0.0:
sight_length = cam_width[2] / sight_angle_cos
else:
- sight_length = sqrt(cam_width[0]**2 + cam_width[1]**2)
- sight_length = sight_length / sqrt(1.0 - sight_angle_cos**2)
+ sight_length = sqrt(cam_width[0] * cam_width[0] +
+ cam_width[1] * cam_width[1])
+ sight_length /= sqrt(1.0 - sight_angle_cos * sight_angle_cos)
fma(sight_length, sight_vector, cam_pos, pos1)
subtract(pos1, sight_center, pos1)
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/marching_cubes.pyx
--- a/yt/utilities/lib/marching_cubes.pyx
+++ b/yt/utilities/lib/marching_cubes.pyx
@@ -18,6 +18,7 @@
import numpy as np
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip
from libc.stdlib cimport malloc, free, abs
+from libc.math cimport sqrt
from fixed_interpolator cimport *
cdef extern from "marching_cubes.h":
@@ -357,7 +358,7 @@
for n in range(3):
temp += normal[n]*normal[n]
# Take the negative, to ensure it points inwardly
- temp = -(temp**0.5)
+ temp = -sqrt(temp)
# Dump this somewhere for now
temp = wval * (fv[0] * normal[0] +
fv[1] * normal[1] +
@@ -368,15 +369,15 @@
for n in range(3):
fv[n] = 0.0
for n in range(3):
- fv[0] += (current.p[0][n] - current.p[2][n])**2.0
- fv[1] += (current.p[1][n] - current.p[0][n])**2.0
- fv[2] += (current.p[2][n] - current.p[1][n])**2.0
+ fv[0] += (current.p[0][n] - current.p[2][n]) * (current.p[0][n] - current.p[2][n])
+ fv[1] += (current.p[1][n] - current.p[0][n]) * (current.p[1][n] - current.p[0][n])
+ fv[2] += (current.p[2][n] - current.p[1][n]) * (current.p[2][n] - current.p[1][n])
s = 0.0
for n in range(3):
- fv[n] = fv[n]**0.5
+ fv[n] = sqrt(fv[n])
s += 0.5 * fv[n]
area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))
- area = area**0.5
+ area = sqrt(area)
flux += temp*area
last = current
if current.next == NULL: break
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -47,7 +47,7 @@
np.ndarray[np.float64_t, ndim=2] qresult,
np.ndarray[np.uint8_t, ndim=1, cast=True] used):
cdef int n, fi, bin
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -57,12 +57,13 @@
wresult[bin] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin,fi]
# qresult has to have the previous wresult
- qresult[bin,fi] += (oldwr * wval * (bval - mresult[bin,fi])**2) / \
+ qresult[bin,fi] += oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin,fi] += wval * (bval - mresult[bin,fi]) / wresult[bin]
+ mresult[bin,fi] += wval * bval_mresult / wresult[bin]
used[bin] = 1
return
@@ -79,7 +80,7 @@
np.ndarray[np.float64_t, ndim=3] qresult,
np.ndarray[np.uint8_t, ndim=2, cast=True] used):
cdef int n, fi, bin_x, bin_y
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -90,12 +91,13 @@
wresult[bin_x,bin_y] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin_x,bin_y,fi]
# qresult has to have the previous wresult
- qresult[bin_x,bin_y,fi] += (oldwr * wval * (bval - mresult[bin_x,bin_y,fi])**2) / \
+ qresult[bin_x,bin_y,fi] += oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin_x,bin_y,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin_x,bin_y,fi] += wval * (bval - mresult[bin_x,bin_y,fi]) / wresult[bin_x,bin_y]
+ mresult[bin_x,bin_y,fi] += wval * bval_mresult / wresult[bin_x,bin_y]
used[bin_x,bin_y] = 1
return
@@ -113,7 +115,7 @@
np.ndarray[np.float64_t, ndim=4] qresult,
np.ndarray[np.uint8_t, ndim=3, cast=True] used):
cdef int n, fi, bin_x, bin_y, bin_z
- cdef np.float64_t wval, bval, oldwr
+ cdef np.float64_t wval, bval, oldwr, bval_mresult
cdef int nb = bins_x.shape[0]
cdef int nf = bsource.shape[1]
for n in range(nb):
@@ -125,14 +127,14 @@
wresult[bin_x,bin_y,bin_z] += wval
for fi in range(nf):
bval = bsource[n,fi]
+ bval_mresult = bval - mresult[bin_x,bin_y,bin_z,fi]
# qresult has to have the previous wresult
qresult[bin_x,bin_y,bin_z,fi] += \
- (oldwr * wval * (bval - mresult[bin_x,bin_y,bin_z,fi])**2) / \
+ oldwr * wval * bval_mresult * bval_mresult / \
(oldwr + wval)
bresult[bin_x,bin_y,bin_z,fi] += wval*bval
# mresult needs the new wresult
- mresult[bin_x,bin_y,bin_z,fi] += wval * \
- (bval - mresult[bin_x,bin_y,bin_z,fi]) / \
+ mresult[bin_x,bin_y,bin_z,fi] += wval * bval_mresult / \
wresult[bin_x,bin_y,bin_z]
used[bin_x,bin_y,bin_z] = 1
return
@@ -149,16 +151,17 @@
np.ndarray[np.float64_t, ndim=1] qresult,
np.ndarray[np.float64_t, ndim=1] used):
cdef int n, bin
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bin = bins_x[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bin] += (wresult[bin] * wval * (bval - mresult[bin])**2) / \
+ bval_mresult = bval - mresult[bin]
+ qresult[bin] += wresult[bin] * wval * bval_mresult * bval_mresult / \
(wresult[bin] + wval)
wresult[bin] += wval
bresult[bin] += wval*bval
- mresult[bin] += wval * (bval - mresult[bin]) / wresult[bin]
+ mresult[bin] += wval * bval_mresult / wresult[bin]
used[bin] = 1
return
@@ -175,17 +178,18 @@
np.ndarray[np.float64_t, ndim=2] qresult,
np.ndarray[np.float64_t, ndim=2] used):
cdef int n, bini, binj
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bini = bins_x[n]
binj = bins_y[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bini, binj] += (wresult[bini, binj] * wval * (bval - mresult[bini, binj])**2) / \
+ bval_mresult = bval - mresult[bini, binj]
+ qresult[bini, binj] += wresult[bini, binj] * wval * bval_mresult * bval_mresult / \
(wresult[bini, binj] + wval)
wresult[bini, binj] += wval
bresult[bini, binj] += wval*bval
- mresult[bini, binj] += wval * (bval - mresult[bini, binj]) / wresult[bini, binj]
+ mresult[bini, binj] += wval * bval_mresult / wresult[bini, binj]
used[bini, binj] = 1
return
@@ -203,18 +207,19 @@
np.ndarray[np.float64_t, ndim=3] qresult,
np.ndarray[np.float64_t, ndim=3] used):
cdef int n, bini, binj, bink
- cdef np.float64_t wval, bval
+ cdef np.float64_t wval, bval, bval_mresult
for n in range(bins_x.shape[0]):
bini = bins_x[n]
binj = bins_y[n]
bink = bins_z[n]
bval = bsource[n]
wval = wsource[n]
- qresult[bini, binj, bink] += (wresult[bini, binj, bink] * wval * (bval - mresult[bini, binj, bink])**2) / \
+ bval_mresult = bval - mresult[bini, binj, bink]
+ qresult[bini, binj, bink] += wresult[bini, binj, bink] * wval * bval_mresult * bval_mresult / \
(wresult[bini, binj, bink] + wval)
wresult[bini, binj, bink] += wval
bresult[bini, binj, bink] += wval*bval
- mresult[bini, binj, bink] += wval * (bval - mresult[bini, binj, bink]) / wresult[bini, binj, bink]
+ mresult[bini, binj, bink] += wval * bval_mresult / wresult[bini, binj, bink]
used[bini, binj, bink] = 1
return
@@ -344,8 +349,8 @@
z1 = zs[j+1]
dx = abs(x1-x0)
dy = abs(y1-y0)
- dzx = (z1-z0) / (dx**2 + dy**2) * dx
- dzy = (z1-z0) / (dx**2 + dy**2) * dy
+ dzx = (z1-z0) / (dx * dx + dy * dy) * dx
+ dzy = (z1-z0) / (dx * dx + dy * dy) * dy
err = dx - dy
if crop == 1 and (dx > nx/2.0 or dy > ny/2.0):
continue
@@ -1016,7 +1021,7 @@
i = 0
n_q = mass.size
pbar = get_pbar("Calculating potential for %d cells" % n_q,
- 0.5 * (n_q**2 - n_q))
+ 0.5 * (n_q * n_q - n_q))
for q_outer in range(n_q - 1):
this_potential = 0.
mass_o = mass[q_outer]
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -352,8 +352,8 @@
for i in range(8):
rbounds[0] = fmin(rbounds[0], corners[i])
rbounds[1] = fmax(rbounds[1], corners[i])
- rbounds[0] = rbounds[0]**0.5
- rbounds[1] = rbounds[1]**0.5
+ rbounds[0] = math.sqrt(rbounds[0])
+ rbounds[1] = math.sqrt(rbounds[1])
# If we include the origin in either direction, we need to have radius of
# zero as our lower bound.
if x0 < 0 and x1 > 0:
@@ -479,8 +479,8 @@
y = (-1.0 + j * dy)*s2
zb = (x*x/8.0 + y*y/2.0 - 1.0)
if zb > 0: continue
- z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- z = z**0.5
+ z = (1.0 - (x * 0.25) * (x * 0.25) - (y * 0.5) * (y * 0.5))
+ z = math.sqrt(z)
# Longitude
theta0 = 2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0)))
# Latitude
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/lib/points_in_volume.pyx
--- a/yt/utilities/lib/points_in_volume.pyx
+++ b/yt/utilities/lib/points_in_volume.pyx
@@ -17,6 +17,7 @@
import numpy as np
cimport numpy as np
cimport cython
+from libc.math cimport sqrt
cdef extern from "math.h":
double fabs(double x)
@@ -123,7 +124,7 @@
cdef np.float64_t norm = 0.0
for i in range(3):
norm += vec[i]*vec[i]
- norm = norm**0.5
+ norm = sqrt(norm)
for i in range(3):
vec[i] /= norm
diff -r 00e84c8e7bfa68ca8bb45a8ec4331acfe5c8d2a1 -r 0874610be59627d1973a759d88c3d1636369fd4e yt/utilities/spatial/ckdtree.pyx
--- a/yt/utilities/spatial/ckdtree.pyx
+++ b/yt/utilities/spatial/ckdtree.pyx
@@ -3,6 +3,7 @@
import numpy as np
cimport numpy as np
cimport libc.stdlib as stdlib
+from libc.math cimport sqrt
cimport cython
import kdtree
@@ -700,7 +701,7 @@
for j in range(num_neighbors):
pj = tags_temp[j]
r2 = dist_temp[j] * ih2
- rs = 2.0 - (r2**0.5)
+ rs = 2.0 - sqrt(r2)
if (r2 < 1.0):
rs = (1.0 - 0.75*rs*r2)
else:
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list