Commit 010e5453 authored by sumpfralle's avatar sumpfralle

fixed code-style issues


git-svn-id: https://pycam.svn.sourceforge.net/svnroot/pycam/trunk@492 bbaffbd6-741e-11dd-a85d-61de82d9cad9
parent 107ff4dc
...@@ -34,7 +34,7 @@ class Iterator: ...@@ -34,7 +34,7 @@ class Iterator:
return v return v
def insertBefore(self, v): def insertBefore(self, v):
self.seq.insert(self.ind-1, v) self.seq.insert(self.ind - 1, v)
self.ind += 1 self.ind += 1
def insert(self, v): def insert(self, v):
...@@ -42,15 +42,15 @@ class Iterator: ...@@ -42,15 +42,15 @@ class Iterator:
self.ind += 1 self.ind += 1
def replace(self, v, w): def replace(self, v, w):
for i in range(0,len(self.seq)): for i in range(len(self.seq)):
if self.seq[i]==v: if self.seq[i] == v:
self.seq[i]=w self.seq[i] = w
def remove(self, v): def remove(self, v):
for i in range(0,len(self.seq)): for i in range(len(self.seq)):
if self.seq[i]==v: if self.seq[i] == v:
del self.seq[i] del self.seq[i]
if i<self.ind: if i < self.ind:
self.ind -= 1 self.ind -= 1
return return
...@@ -64,13 +64,13 @@ class Iterator: ...@@ -64,13 +64,13 @@ class Iterator:
return Iterator(self.seq, self.ind) return Iterator(self.seq, self.ind)
def peek(self, i=0): def peek(self, i=0):
if self.ind+i >= len(self.seq): if self.ind + i >= len(self.seq):
return None return None
else: else:
return self.seq[self.ind+i] return self.seq[self.ind + i]
def remains(self): def remains(self):
return len(self.seq)-self.ind return len(self.seq) - self.ind
class CyclicIterator: class CyclicIterator:
def __init__(self, seq, start=0): def __init__(self, seq, start=0):
...@@ -89,8 +89,8 @@ class CyclicIterator: ...@@ -89,8 +89,8 @@ class CyclicIterator:
return CyclicIterator(self.seq, self.ind) return CyclicIterator(self.seq, self.ind)
def peek(self, i=0): def peek(self, i=0):
idx = self.ind+i idx = self.ind + i
while idx>=len(self.seq): while idx >= len(self.seq):
idx -= len(self.seq) idx -= len(self.seq)
return self.seq[idx] return self.seq[idx]
...@@ -100,10 +100,10 @@ if __name__ == "__main__": ...@@ -100,10 +100,10 @@ if __name__ == "__main__":
i = Iterator(l) i = Iterator(l)
print i.peek() print i.peek()
while True: while True:
v = i.next() val = i.next()
if v == None: if val == None:
break break
if v == 4: if val == 4:
i.insertBefore(3) i.insertBefore(3)
i.insert(5) i.insert(5)
...@@ -118,14 +118,14 @@ if __name__ == "__main__": ...@@ -118,14 +118,14 @@ if __name__ == "__main__":
print "remains=", i.remains() print "remains=", i.remains()
print "l=", l print "l=", l
sum = 0 sum_value = 0
i = CyclicIterator(l) i = CyclicIterator(l)
print "cycle :", print "cycle :",
while sum<30: while sum_value < 30:
v = i.next() val = i.next()
print v, print val,
sum += v sum_value += val
print "=", sum print "=", sum_value
i = Iterator(l) i = Iterator(l)
print "l=", l print "l=", l
...@@ -137,3 +137,4 @@ if __name__ == "__main__": ...@@ -137,3 +137,4 @@ if __name__ == "__main__":
i.remove(4) i.remove(4)
print "remove(4) : ", i.peek() print "remove(4) : ", i.peek()
print "l=", l print "l=", l
...@@ -34,7 +34,7 @@ def get_logger(suffix=None): ...@@ -34,7 +34,7 @@ def get_logger(suffix=None):
def init_logger(log, logfilename=None): def init_logger(log, logfilename=None):
if logfilename: if logfilename:
datetime_format = "%(asctime)s - %(name)s - %(levelname)s - %(message)s" datetime_format = "%(asctime)s - %(name)s - %(levelname)s - %(message)s"
logfile_hander = logging.FileHandler(logfilename) logfile_handler = logging.FileHandler(logfilename)
logfile_handler.setFormatter(datetime_format) logfile_handler.setFormatter(datetime_format)
log.addHandler(logfile_handler) log.addHandler(logfile_handler)
console_output = logging.StreamHandler() console_output = logging.StreamHandler()
......
...@@ -20,75 +20,76 @@ You should have received a copy of the GNU General Public License ...@@ -20,75 +20,76 @@ You should have received a copy of the GNU General Public License
along with PyCAM. If not, see <http://www.gnu.org/licenses/>. along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
""" """
from math import * import math
from math import sqrt
# see BRL-CAD/src/libbn/poly.c # see BRL-CAD/src/libbn/poly.c
EPSILON=1e-4 EPSILON = 1e-4
SMALL=1e-4 SMALL = 1e-4
INV_2=0.5 INV_2 = 0.5
INV_3=1.0/3.0 INV_3 = 1.0 / 3.0
INV_4=0.25 INV_4 = 0.25
INV_27=1.0/27.0 INV_27 = 1.0 / 27.0
SQRT3=sqrt(3.0) SQRT3 = sqrt(3.0)
PI_DIV_3=pi/3.0 PI_DIV_3 = math.pi / 3.0
def near_zero(x, epsilon=EPSILON): def near_zero(x, epsilon=EPSILON):
if x>-epsilon and x<epsilon: if abs(x) < epsilon:
return True return True
else: else:
return False return False
def cuberoot(x): def cuberoot(x):
if x>=0: if x >= 0:
return pow(x,INV_3) return pow(x, INV_3)
else: else:
return -pow(-x,INV_3) return -pow(-x, INV_3)
def poly1_roots(a,b): def poly1_roots(a, b):
if near_zero(a): if near_zero(a):
return None return None
return (-b/a,) return (-b / a, )
def poly2_roots(a,b,c): def poly2_roots(a, b, c):
d = b*b-4*a*c d = b * b - 4 * a * c
if d < 0: if d < 0:
return None return None
if near_zero(a): if near_zero(a):
return poly1_roots(b,c) return poly1_roots(b, c)
if d == 0: if d == 0:
return (-b/(2*a), ) return (-b / (2 * a), )
q = sqrt(d) q = sqrt(d)
if a<0: if a < 0:
return ((-b+q)/(2*a),(-b-q)/(2*a)) return ((-b + q) / (2 * a), (-b - q) / (2 * a))
else: else:
return ((-b-q)/(2*a),(-b+q)/(2*a)) return ((-b - q) / (2 * a), (-b + q) / (2 * a))
def poly3_roots(a,b,c,d): def poly3_roots(a, b, c, d):
if near_zero(a): if near_zero(a):
return poly2_roots(b,c,d) return poly2_roots(b, c, d)
c1=b/a c1 = b / a
c2=c/a c2 = c / a
c3=d/a c3 = d / a
c1_3=c1*INV_3 c1_3 = c1 * INV_3
a = c2-c1*c1_3 a = c2 -c1 * c1_3
b = (2*c1*c1*c1-9*c1*c2+27*c3)*INV_27 b = (2 * c1 * c1 * c1 - 9 * c1 * c2 + 27 * c3) * INV_27
delta = a*a delta = a * a
delta = b*b*INV_4+delta*a*INV_27 delta = b * b * INV_4 + delta * a * INV_27
if delta>0: if delta > 0:
r_delta = sqrt(delta) r_delta = sqrt(delta)
A = -INV_2*b + r_delta A = -INV_2 * b + r_delta
B = -INV_2*b - r_delta B = -INV_2 * b - r_delta
A = cuberoot(A) A = cuberoot(A)
B = cuberoot(B) B = cuberoot(B)
return (A+B-c1_3, ) return (A + B - c1_3, )
elif delta==0: elif delta == 0:
b_2 = -b*INV_2 b_2 = -b * INV_2
s = cuberoot(b_2) s = cuberoot(b_2)
return (2*s-c1_3,-s-c1_3,-s-c1_3,) return (2 * s - c1_3, -s - c1_3, -s - c1_3, )
else: else:
if a>0: if a > 0:
fact = 0 fact = 0
phi = 0 phi = 0
cs_phi = 1.0 cs_phi = 1.0
...@@ -96,75 +97,66 @@ def poly3_roots(a,b,c,d): ...@@ -96,75 +97,66 @@ def poly3_roots(a,b,c,d):
else: else:
a *= -INV_3 a *= -INV_3
fact = sqrt(a) fact = sqrt(a)
f = -b*INV_2/(a*fact) f = -b * INV_2 / (a * fact)
if f >= 1.0: if f >= 1.0:
phi = 0 phi = 0
cs_phi = 1.0 cs_phi = 1.0
sn_phi_s3 = 0.0 sn_phi_s3 = 0.0
elif f <= -1.0: elif f <= -1.0:
phi = PI_DIV_3 phi = PI_DIV_3
cs_phi = cos(phi) cs_phi = math.cos(phi)
sn_phi_s3 = sin(phi)*SQRT3 sn_phi_s3 = math.sin(phi) * SQRT3
else: else:
phi = acos(f)*INV_3 phi = math.acos(f) * INV_3
cs_phi = cos(phi) cs_phi = math.cos(phi)
sn_phi_s3 = sin(phi)*SQRT3 sn_phi_s3 = math.sin(phi) * SQRT3
r1 = 2*fact*cs_phi r1 = 2 * fact * cs_phi
r2 = fact*( sn_phi_s3 - cs_phi) r2 = fact * ( sn_phi_s3 - cs_phi)
r3 = fact*(-sn_phi_s3 - cs_phi) r3 = fact * (-sn_phi_s3 - cs_phi)
return (r1-c1_3, r2-c1_3, r3-c1_3) return (r1 - c1_3, r2 - c1_3, r3 - c1_3)
def max3(a,b,c): def poly4_roots(a, b, c, d, e):
if a>b: if a == 0:
max_ab = a return poly3_roots(b, c, d, e)
else: c1 = float(b) / a
max_ab = b c2 = float(c) / a
if c>max_ab: c3 = float(d) / a
return c c4 = float(e) / a
else: roots3 = poly3_roots(1.0, -c2, c3 * c1 - 4 * c4, -c3 * c3 - c4 * c1 * c1 \
return max_ab + 4 * c4 * c2)
def poly4_roots(a,b,c,d,e):
if a==0:
return poly3_roots(b,c,d,e)
c1 = float(b)/a
c2 = float(c)/a
c3 = float(d)/a
c4 = float(e)/a
roots3 = poly3_roots(1.0, -c2, c3*c1-4*c4, -c3*c3-c4*c1*c1+4*c4*c2)
if not roots3: if not roots3:
return None return None
if len(roots3)==1: if len(roots3) == 1:
U = roots3[0] U = roots3[0]
else: else:
U = max3(roots3[0],roots3[1],roots3[2]) U = max(roots3[0], roots3[1], roots3[2])
p = c1*c1*INV_4+U-c2 p = c1 * c1 * INV_4 + U - c2
U *= INV_2 U *= INV_2
q = U*U-c4 q = U * U - c4
if p<0: if p < 0:
if p<-SMALL: if p < -SMALL:
return None return None
p = 0 p = 0
else: else:
p = sqrt(p) p = sqrt(p)
if q<0: if q < 0:
if q<-SMALL: if q < -SMALL:
return None return None
q = 0 q = 0
else: else:
q = sqrt(q) q = sqrt(q)
quad1 = [1.0,c1*INV_2-p,0] quad1 = [1.0, c1 * INV_2 - p, 0]
quad2 = [1.0,c1*INV_2+p,0] quad2 = [1.0, c1 * INV_2 + p, 0]
q1 = U - q q1 = U - q
q2 = U + q q2 = U + q
p = quad1[1]*q2+quad2[1]*q1-c3 p = quad1[1] * q2 + quad2[1] * q1 - c3
if near_zero(p): if near_zero(p):
quad1[2] = q1 quad1[2] = q1
quad2[2] = q2 quad2[2] = q2
else: else:
q = quad1[1]*q1+quad2[1]*q2-c3 q = quad1[1] * q1 + quad2[1] * q2 - c3
if near_zero(q): if near_zero(q):
quad1[2] = q2 quad1[2] = q2
quad2[2] = q1 quad2[2] = q1
...@@ -175,7 +167,7 @@ def poly4_roots(a,b,c,d,e): ...@@ -175,7 +167,7 @@ def poly4_roots(a,b,c,d,e):
roots1 = poly2_roots(quad1[0], quad1[1], quad1[2]) roots1 = poly2_roots(quad1[0], quad1[1], quad1[2])
roots2 = poly2_roots(quad2[0], quad2[1], quad2[2]) roots2 = poly2_roots(quad2[0], quad2[1], quad2[2])
if roots1 and roots2: if roots1 and roots2:
return roots1+roots2 return roots1 + roots2
elif roots1: elif roots1:
return roots1 return roots1
elif roots2: elif roots2:
...@@ -183,65 +175,65 @@ def poly4_roots(a,b,c,d,e): ...@@ -183,65 +175,65 @@ def poly4_roots(a,b,c,d,e):
else: else:
return None return None
def test_poly1(a,b): def test_poly1(a, b):
roots = poly1_roots(a,b) roots = poly1_roots(a, b)
print a, "*x+",b,"=0 ", roots print a, "*x+", b, "=0 ", roots
if roots: if roots:
for r in roots: for r in roots:
f = a*r+b f = a * r + b
if not near_zero(f): if not near_zero(f):
print "ERROR:", print "ERROR:",
print " f(",r,") =", f print " f(%f)=%f" % (r, f)
def test_poly2(a,b,c): def test_poly2(a, b, c):
roots = poly2_roots(a,b,c) roots = poly2_roots(a, b, c)
print a, "*x^2+",b,"*x+",c,"=0 ", roots print a, "*x^2+", b, "*x+", c, "=0 ", roots
if roots: if roots:
for r in roots: for r in roots:
f = a*r*r+b*r+c f = a * r * r + b * r + c
if not near_zero(f): if not near_zero(f):
print "ERROR:", print "ERROR:",
print " f(",r,") = ", f print " f(%f)=%f" % (r, f)
def test_poly3(a,b,c,d): def test_poly3(a, b, c, d):
roots = poly3_roots(a,b,c,d) roots = poly3_roots(a, b, c, d)
print a, "*x^3+",b,"*x^2+",c,"*x+",d,"=0 ", roots print a, "*x^3+", b, "*x^2+", c, "*x+", d, "=0 ", roots
if roots: if roots:
for r in roots: for r in roots:
f = a*r*r*r+b*r*r+c*r+d f = a * r * r * r + b * r * r + c * r + d
if not near_zero(f): if not near_zero(f):
print "ERROR:", print "ERROR:",
print " f(",r,") = ", f print " f(%f)=%f" % (r, f)
def test_poly4(a,b,c,d,e): def test_poly4(a, b, c, d, e):
roots = poly4_roots(a,b,c,d,e) roots = poly4_roots(a, b, c, d, e)
print "f(x)=%g*x**4%+g*x**3%+g*x**2%+g*x%+g" % (a,b,c,d,e) print "f(x)=%g*x**4%+g*x**3%+g*x**2%+g*x%+g" % (a, b, c, d, e)
print "roots:", roots print "roots:", roots
if roots: if roots:
for r in roots: for r in roots:
f = a*r*r*r*r+b*r*r*r+c*r*r+d*r+e f = a * r * r * r * r + b * r * r * r + c * r * r + d * r + e
if not near_zero(f, epsilon=SMALL): if not near_zero(f, epsilon=SMALL):
print "ERROR:", print "ERROR:",
print " f(",r,") = ", f print " f(%f)=%f" % (r, f)
return roots return roots
if __name__ == "__main__": if __name__ == "__main__":
test_poly1(1,2) test_poly1(1, 2)
test_poly2(1,2,0) test_poly2(1, 2, 0)
test_poly2(1,2,1) test_poly2(1, 2, 1)
test_poly2(1,2,2) test_poly2(1, 2, 2)
test_poly3(1,0,0,0) test_poly3(1, 0, 0, 0)
test_poly3(1,0,0,-1) test_poly3(1, 0, 0, -1)
test_poly3(1,-1,0,0) test_poly3(1, -1, 0, 0)
test_poly3(1,0,-2,0) test_poly3(1, 0, -2, 0)
test_poly3(1,0,-2,1) test_poly3(1, 0, -2, 1)
test_poly4(1,0,0,0,0) test_poly4(1, 0, 0, 0, 0)
test_poly4(1,0,0,0,-1) test_poly4(1, 0, 0, 0, -1)
test_poly4(1,0,-2,0,1) test_poly4(1, 0, -2, 0, 1)
test_poly4(1,-10,35,-50,+24) test_poly4(1, -10, 35, -50, +24)
test_poly4(1,0,6,-60,36) test_poly4(1, 0, 6, -60, 36)
test_poly4(1,-25,235.895,-995.565,1585.25); test_poly4(1, -25, 235.895, -995.565, 1585.25)
...@@ -23,42 +23,41 @@ along with PyCAM. If not, see <http://www.gnu.org/licenses/>. ...@@ -23,42 +23,41 @@ along with PyCAM. If not, see <http://www.gnu.org/licenses/>.
def find_root_subdivide(f, x0, x1, tolerance, scale): def find_root_subdivide(f, x0, x1, tolerance, scale):
ymin = 0 ymin = 0
xmin = 0 xmin = 0
imin = 0 while x1 - x0 > tolerance:
while x1-x0>tolerance: for i in range(scale):
for i in range(0,scale): x = x1 + (i / scale) * (x1 - x0)
x = x1 + (i/scale)*(x1-x0)
y = f(x) y = f(x)
abs_y = abs(y) abs_y = abs(y)
if i==0: if i == 0:
ymin = abs_y ymin = abs_y
xmin = x xmin = x
imin = 0
else: else:
if abs_y<ymin: if abs_y < ymin:
ymin = abs_y ymin = abs_y
xmin = x xmin = x
imin = i x0 = xmin - 1 / scale
x0 = xmin - 1/scale x1 = xmin + 1 / scale
x1 = xmin + 1/scale
scale /= 10 scale /= 10
return xmin return xmin
def find_root_newton_raphson(f, df, x0, tolerance, maxiter): def find_root_newton_raphson(f, df, x0, tolerance, maxiter):
x = x0 x = x0
iter = 0 iter_count = 0
while iter<maxiter: while iter_count < maxiter:
y = f(x) y = f(x)
if y == 0: if y == 0:
return x return x
dy = df(x) dy = df(x)
if dy == 0: if dy == 0:
return None return None
dx = y/dy dx = y / dy
x = x - dx x = x - dx
if dx < tolerance: if dx < tolerance:
break break
iter += 1 iter_count += 1
return x return x
def find_root(f, df=None, x0=0, x1=1, tolerance=0.001): def find_root(f, df=None, x0=0, x1=1, tolerance=0.001):
return find_root_subdivide(f=f,x0=x0,x1=x1,tolerance=tolerance, scale=10.0) return find_root_subdivide(f=f, x0=x0, x1=x1, tolerance=tolerance,
scale=10.0)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment