suzanne.soy/geosolver
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

matfunc.py - preparing forchange from vector ->matfunc.

Browse Source 
intersection.py - added tests (preparing for change vector->matfunc)
clsolver2D.py - added stub for HH2S
This commit is contained in:
kwikrick 2011-08-16 20:32:20 +00:00
parent 69ee114ddc
commit 9d6002102f
5 changed files with 605 additions and 102 deletions

102
TODO.txt
View File

@ -63,3 +63,105 @@ def diamond_3d():
or:
FixConstraint(v2=[-5.0, 5.0, 0.0]) not satisfied
FixConstraint(v1=[0.0, 0.0, 0.0]) not satisfied
- when solvin DAD type clusters, sometimes raises:
StandardError: cluster rigid#157388428(['v1', 'v2', 'v3']) already in clsolver
Debug output:
Solving....
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: a = v3
triplet2dad: c = v1
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: a = v3
triplet2dad: c = v1
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: start
triplet2dad: b = v2
triplet2dad: b in rigids
Traceback (most recent call last):
File "test.py", line 971, in <module>
t2d()
File "test.py", line 968, in test2d
test(dad_problem())
File "test.py", line 780, in test
solver = GeometricSolver(problem, use_prototype)
File "/home/rick/lib/python/geosolver/geometric.py", line 294, in __init__
self._add_constraint(con)
File "/home/rick/lib/python/geosolver/geometric.py", line 486, in _add_constraint
self.dr.add(hog)
File "/home/rick/lib/python/geosolver/clsolver.py", line 79, in add
self._process_new()
File "/home/rick/lib/python/geosolver/clsolver.py", line 363, in _process_new
self._add_method_complete(method)
File "/home/rick/lib/python/geosolver/clsolver.py", line 487, in _add_method_complete
self._add_cluster(output)
File "/home/rick/lib/python/geosolver/clsolver.py", line 236, in _add_cluster
raise StandardError, "cluster %s already in clsolver"%(str(newcluster))
StandardError: cluster rigid#157388428(['v1', 'v2', 'v3']) already in clsolver

412
geosolver/clsolver2D.py
View File

@ -16,7 +16,7 @@ class ClusterSolver2D(ClusterSolver):
def __init__(self):
"""Instantiate a ClusterSolver2D"""
ClusterSolver.__init__(self, [CheckAR, MergePR, MergeRR, DeriveDDD, DeriveDAD])
ClusterSolver.__init__(self, [CheckAR, MergePR, MergeRR, DeriveDDD, DeriveDAD, DeriveADD, DeriveHH2S])
# ----------------------------------------------
@ -169,7 +169,7 @@ class DeriveDDD(ClusterMethod):
#patterngraph = _pattern()
def _incremental_matcher(solver):
triplets = Triplets(solver, Rigids(solver))
triplets = DistanceTriplets(solver, Rigids(solver))
matcher = incremental.Map(triplet2ddd, triplets)
return matcher
@ -244,32 +244,33 @@ class DeriveDAD(ClusterMethod):
def isdad(triplet):
dad = triplet2dad(triplet)
return isinstance(DeriveDAD, dad)
return isinstance(dad, DeriveDAD)
def triplet2dad(triplet):
print "triplet2dad: start"
hogs = filter(lambda c: isinstance(Hog, c), triplet)
rigids= filter(lambda c: isinstance(Rigid, c), triplet)
hogs = filter(lambda c: isinstance(c, Hedgehog), triplet)
rigids= filter(lambda c: isinstance(c, Rigid), triplet)
if not(len(hogs)==1 and len(rigids)==2): return None
hog = hogs[0]
r1 = rigids[0]
r2 = rigids[2]
b = hog.apex;
r2 = rigids[1]
b = hog.cvar;
print "triplet2dad: b = ", b
if not(b in r1.vars): return None
if not(b in r2.vars): return None
print "triplet2dad: b in rigids"
p1s = r1.vars.intersection(hog.vars)
p2s = r2.vars.intersection(hog.vars)
p1s = r1.vars.intersection(hog.xvars)
p2s = r2.vars.intersection(hog.xvars)
if not(len(p1s) == 1): return None
if not(len(p2s) == 1): return None
a = p1s[0]
b = p1s[1]
a = list(p1s)[0]
c = list(p2s)[0]
print "triplet2dad: a = ", a
print "triplet2dad: c = ", c
if a==c: return None
return DeriveDAD( {"$d_ab":r1, "$a_abc":hog, "$d_bc":r2, "$a":a, "$b":b, "$c":c })
# end def
triplets = Triplets(solver, solver.top_level())
triplets = ConnectedTriplets(solver, solver.top_level())
matchtriplets = incremental.Filter(lambda triplet: isdad(triplet), triplets)
matcher = incremental.Map(triplet2dad, matchtriplets)
return matcher
@ -296,6 +297,107 @@ class DeriveDAD(ClusterMethod):
solutions = solve_dad(v1,v2,v3,d12,a123,d23)
return solutions
class DeriveADD(ClusterMethod):
"""Represents a merging of two distances and an angle"""
def __init__(self, map):
# check inputs
self.a_cab = map["$a_cab"]
self.d_ab = map["$d_ab"]
self.d_bc = map["$d_bc"]
self.a = map["$a"]
self.b = map["$b"]
self.c = map["$c"]
# create output
out = Rigid([self.a,self.b,self.c])
# set method parameters
self._inputs = [self.a_cab, self.d_ab, self.d_bc]
self._outputs = [out]
ClusterMethod.__init__(self)
# do not remove input clusters (because root not considered here)
self.noremove = True
#def _pattern():
# pattern = [["rigid","$d_ab",["$a", "$b"]],
# ["hedgehog", "$a_abc",["$b", "$a", "$c"]],
# ["rigid", "$d_bc",["$b","$c"]]]
# return pattern2graph(pattern)
#pattern = staticmethod(_pattern)
#patterngraph = _pattern()
def _incremental_matcher(solver):
def isadd(triplet):
dad = triplet2add(triplet)
return isinstance(dad, DeriveADD)
def triplet2add(triplet):
print "triplet2add: start"
hogs = filter(lambda c: isinstance(c, Hedgehog), triplet)
rigids= filter(lambda c: isinstance(c, Rigid), triplet)
if not(len(hogs)==1 and len(rigids)==2): return None
hog = hogs[0]
print "triplet2add: hog = ",hog
r1 = rigids[0]
r2 = rigids[1]
a = hog.cvar;
print "triplet2add: a = ", a
if a in r1.vars and not(a in r2.vars):
d_ab = r1
d_bc = r2
elif a in r2.vars and not(a in r1.vars):
d_ab = r2
d_bc = r1
else:
return None
print "d_ab:",d_ab
print "d_bc:",d_bc
pbs = d_ab.vars.intersection(hog.xvars)
if not(len(pbs) == 1): return None
b = list(pbs)[0]
print "triplet2add: b = ", b
pcs = d_bc.vars.intersection(hog.xvars).difference([b])
if not(len(pcs) == 1): return None
c = list(pcs)[0]
print "triplet2add: c = ", c
return DeriveADD( {"$a_cab":hog, "$d_ab":d_ab, "$d_bc":d_bc, "$a":a, "$b":b, "$c":c })
# end def
triplets = ConnectedTriplets(solver, solver.top_level())
matchtriplets = incremental.Filter(lambda triplet: isadd(triplet), triplets)
matcher = incremental.Map(triplet2add, matchtriplets)
return matcher
incremental_matcher = staticmethod(_incremental_matcher)
def __str__(self):
s = "DeriveADD("+str(self._inputs[0])+"+"+str(self._inputs[1])+"+"+str(self._inputs[2])+"->"+str(self._outputs[0])+")"
s += "[" + self.status_str()+"]"
return s
def multi_execute(self, inmap):
diag_print("DeriveADD.multi_execute called","clmethods")
c312 = inmap[self.a_cab]
c12 = inmap[self.d_ab]
c23 = inmap[self.d_bc]
v1 = self.a
v2 = self.b
v3 = self.c
a312 = angle_3p(c312.get(v3),c312.get(v1),c312.get(v2))
d12 = distance_2p(c12.get(v1),c12.get(v2))
d23 = distance_2p(c23.get(v2),c23.get(v3))
solutions = solve_add(v1,v2,v3,a312,d12,d23)
return solutions
def prototype_constraints(self):
constraints = []
constraints.append(SelectionConstraint(fnot(is_obtuse),[self.a,self.c,self.b]))
constraints.append(SelectionConstraint(fnot(is_acute),[self.a,self.c,self.b]))
return constraints
class CheckAR(ClusterMethod):
"""Represents the overconstrained merging a hedgehog and a rigid that completely overlaps it."""
def __init__(self, map):
@ -352,6 +454,63 @@ class CheckAR(ClusterMethod):
# all checks passed, return rigid configuration
return [rigid]
class DeriveHH2S(ClusterMethod):
"""Derive a scalable from two angles"""
def __init__(self, map):
# check inputs
self.a_cab = map["$a_cab"]
self.a_abc = map["$a_abc"]
self.a = map["$a"]
self.b = map["$b"]
self.c = map["$c"]
# create output
out = Balloon([self.a,self.b,self.c])
# set method parameters
self._inputs = [self.a_cab, self.a_abc]
self._outputs = [out]
ClusterMethod.__init__(self)
# do not remove input clusters (because root not considered here)
self.noremove = True
def __str__(self):
s = "DeriveHH2S("+str(self._inputs[0])+"+"+str(self._inputs[1])+"->"+str(self._outputs[0])+")"
s += "[" + self.status_str()+"]"
return s
def _incremental_matcher(solver):
def pair_is_hh2s(pair):
method = pair_to_hh2s(pair)
return isinstance(method, DeriveHH2S)
def pair_to_hh2s(pair):
print "pair_to_hhs2s: start"
print "not implemented!"
return None
#return DeriveHH2S( {"$a_cab":a_cab, "$a_abc":a_abc, "$a":a, "$b":b, "$c":c })
# end def
hogs = Hogs(solver)
pairs = ConnectedPairs1(solver, hogs)
matchingpairs = incremental.Filter(lambda pair: pair_is_hh2s(pair), pairs)
matcher = incremental.Map(pair_to_hh2s, matchingpairs)
return matcher
incremental_matcher = staticmethod(_incremental_matcher)
def multi_execute(self, inmap):
diag_print("DeriveHH2S.multi_execute called","clmethods")
c312 = inmap[self.a_cab]
c123 = inmap[self.a_abc]
v1 = self.a
v2 = self.b
v3 = self.c
a312 = angle_3p(c312.get(v3),c312.get(v1),c312.get(v2))
d12 = 1.0
a123 = angle_3p(c123.get(v1),c123.get(v2),c123.get(v3))
solutions = solve_ada_3D(v1,v2,v3,a312,d12,a123)
return solutions
# ---------------------------------------------------------
# ------- functions to determine configurations ----------
@ -385,11 +544,123 @@ def solve_dad(v1,v2,v3,d12,a123,d23):
solutions.append(solution)
return solutions
def solve_add(a,b,c, a_cab, d_ab, d_bc):
"""returns a list of Configurations of v1,v2,v3 such that distance v1-v2=d12 etc.
v<x>: name of point variables
d<xy>: numeric distance values
a<xyz>: numeric angle in radians
"""
diag_print("solve_dad: %s %s %s %f %f %f"%(a,b,c,a_cab,d_ab,d_bc),"clmethods")
p_a = vector.vector([0.0,0.0])
p_b = vector.vector([d_ab,0.0])
dir = vector.vector([math.cos(-a_cab),math.sin(-a_cab)])
solutions = cr_int(p_b, d_bc, p_a, dir)
rval = []
for p_c in solutions:
map = {a:p_a, b:p_b, c:p_c}
rval.append(Configuration(map))
return rval
# -------------------------------------
# --------- incremental sets ----------
# -------------------------------------
class Rigids(incremental.Filter):
def __init__(self, solver):
self._solver = solver
incremental.Filter.__init__(self, lambda c: isinstance(c, Rigid), self._solver.top_level())
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Rigids("+repr(self._solver)+")"
class Hogs(incremental.Filter):
def __init__(self, solver):
self._solver = solver
incremental.Filter.__init__(self, lambda c: isinstance(c, Hedgehog), self._solver.top_level())
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Hogs("+repr(self._solver)+")"
class Balloons(incremental.Filter):
def __init__(self, solver):
self._solver = solver
incremental.Filter.__init__(self, lambda c: isinstance(c, Balloon), self._solver.top_level())
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Balloons("+repr(self._solver)+")"
class Points(incremental.Filter):
def __init__(self, solver):
self._solver = solver
rigids = Rigids(solver)
incremental.Filter.__init__(self, lambda c: len(c.vars)==1, rigids)
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Points("+repr(self._solver)+")"
class Distances(incremental.Filter):
def __init__(self, solver):
self._solver = solver
rigids = Rigids(solver)
incremental.Filter.__init__(self, lambda c: len(c.vars)==2, rigids)
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Distances("+repr(self._solver)+")"
class ConnectedPairs1(incremental.IncrementalSet):
"""Incremental set of all pairs of connected clusters in 1 incremental set"""
@ -467,69 +738,10 @@ class ConnectedPairs2(incremental.IncrementalSet):
def __hash__(self):
return hash((self._solver, self._incrset1, self._incrset2))
class Rigids(incremental.Filter):
def __init__(self, solver):
self._solver = solver
incremental.Filter.__init__(self, lambda c: isinstance(c, Rigid), self._solver.top_level())
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Rigids("+repr(self._solver)+")"
class Points(incremental.Filter):
def __init__(self, solver):
self._solver = solver
rigids = Rigids(solver)
incremental.Filter.__init__(self, lambda c: len(c.vars)==1, rigids)
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Points("+repr(self._solver)+")"
class Distances(incremental.Filter):
def __init__(self, solver):
self._solver = solver
rigids = Rigids(solver)
incremental.Filter.__init__(self, lambda c: len(c.vars)==2, rigids)
def __hash__(self):
return hash((self.__class__, self._solver))
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver
else:
return False
def __repr__(self):
return "Distances("+repr(self._solver)+")"
class Triplets(incremental.IncrementalSet):
class DistanceTriplets(incremental.IncrementalSet):
def __init__(self, solver, incrset):
"""Creates an incremental set of all tripltes of connected clusters in incrset, according to solver"""
"""Creates an incremental set of all tripltes of 1-connected clusters in incrset, according to solver"""
self._solver = solver
self._incrset = incrset
incremental.IncrementalSet.__init__(self, [incrset])
@ -574,6 +786,56 @@ class Triplets(incremental.IncrementalSet):
return hash((self.__class__, self._solver, self._incrset))
def __repr__(self):
return "Triplets("+repr(self._solver)+","+repr(self._incset)+")"
return "DistanceTriplets("+repr(self._solver)+","+repr(self._incset)+")"
class ConnectedTriplets(incremental.IncrementalSet):
def __init__(self, solver, incrset):
"""Creates an incremental set of all triplets of connected clusters in incrset, according to solver"""
self._solver = solver
self._incrset = incrset
incremental.IncrementalSet.__init__(self, [incrset])
return
def _receive_add(self,source, obj):
connected = set()
for var in obj.vars:
dependend = self._solver.find_dependend(var)
dependend = filter(lambda x: x in self._incrset, dependend)
dependend = filter(lambda x: len(x.vars.intersection(obj.vars))>=1, dependend)
connected.update(dependend)
if obj in connected:
connected.remove(obj)
obj1 = obj
if len(connected) >= 2:
l = list(connected)
for i in range(len(l)):
obj2 = l[i]
shared12 = obj1.vars.intersection(obj2.vars)
for j in range(i):
obj3 = l[j]
shared23 = obj2.vars.intersection(obj3.vars)
if len(shared23)>=1:
shared13 = obj1.vars.intersection(obj3.vars)
if len(shared13)>=1:
self._add(frozenset((obj1,obj2,obj3)))
def _receive_remove(self,source, obj):
for frozen in list(self):
if obj in frozen:
self._remove(frozen)
def __eq__(self, other):
if isinstance(other, self.__class__):
return self._solver == other._solver and self._incrset == other._incrset
else:
return False
def __hash__(self):
return hash((self.__class__, self._solver, self._incrset))
def __repr__(self):
return "ConnectedTriplets("+repr(self._solver)+","+repr(self._incset)+")"

173
geosolver/intersections.py
View File

@ -1,5 +1,6 @@
import vector
import math
import random
from matfunc import Mat,Vec,eye
from tolerance import *
from diagnostic import *
@ -15,6 +16,14 @@ def sign(x):
else:
return 0.0
def sign2(x):
"""Returns 1 if x>0, else -1 (even if x=0)"""
if tol_gt(x,0.0):
return 1.0
else:
return -1.0
# -------- 3D intersections ---------------
@ -118,15 +127,14 @@ def cl_int(p1,r,p2,v):
E = r*r*d2 - D*D
if tol_gt(d2, 0) and tol_gt(E, 0):
sE = math.sqrt(E)
x1 = p1[0] + (D * v[1] + sign(v[1])*v[0]*sE) / d2
x2 = p1[0] + (D * v[1] - sign(v[1])*v[0]*sE) / d2
x1 = p1[0] + (D * v[1] + sign2(v[1])*v[0]*sE) / d2
x2 = p1[0] + (D * v[1] - sign2(v[1])*v[0]*sE) / d2
y1 = p1[1] + (-D * v[0] + abs(v[1])*sE) / d2
y2 = p1[1] + (-D * v[0] - abs(v[1])*sE) / d2
return [vector.vector([x1,y1]), vector.vector([x2,y2])]
elif tol_eq(E, 0):
x1 = p1[0] + D * v[1] / d2
y1 = p1[1] + -D * v[0] / d2
# return [vector.vector([x1,y1]), vector.vector([x1,y1])]
return [vector.vector([x1,y1])]
else:
return []
@ -608,8 +616,118 @@ def test_sss_int():
# print sat
return sat
def test1():
#diag_select(".*")
def test_cc_int():
"""Generates two random circles, computes the intersection,
and verifies that the number of intersections and the
positions of the intersections are correct.
Returns True or False"""
# gen two random circles p1,r2 and p2, r2
p1 = vector.randvec(2, 0.0, 10.0,1.0)
r1 = random.uniform(0, 10.0)
p2 = vector.randvec(2, 0.0, 10.0,1.0)
r2 = random.uniform(0, 10.0)
# 33% change that r2=abs(r1-|p1-p2|)
if random.random() < 0.33:
r2 = abs(r1-vector.norm(p1-p2))
# do interesection
diag_print("problem:"+str(p1)+","+str(r1)+","+str(p2)+","+str(r2),"test_cc_int")
sols = cc_int(p1, r1, p2, r2)
diag_print("solutions:"+str(map(str, sols)),"test_cc_int")
# test number of intersections
if len(sols) == 0:
if not tol_gt(vector.norm(p2-p1),r1 + r2) and not tol_lt(vector.norm(p2-p1),abs(r1 - r2)) and not tol_eq(vector.norm(p1-p2),0):
diag_print("number of solutions 0 is wrong","test_cc_int")
return False
elif len(sols) == 1:
if not tol_eq(vector.norm(p2-p1),r1 + r2) and not tol_eq(vector.norm(p2-p1),abs(r1-r2)):
diag_print("number of solutions 1 is wrong","test_cc_int")
return False
elif len(sols) == 2:
if not tol_lt(vector.norm(p2-p1),r1 + r2) and not tol_gt(vector.norm(p2-p1),abs(r1 - r2)):
diag_prin("number of solutions 2 is wrong")
return False
else:
diag_print("number of solutions > 2 is wrong","test_cc_int")
return False
# test intersection coords
for p3 in sols:
if not tol_eq(vector.norm(p3-p1), r1):
diag_print("solution not on circle 1","test_cc_int")
return False
if not tol_eq(vector.norm(p3-p2), r2):
diag_print("solution not on circle 2","test_cc_int")
return False
diag_print("OK","test_cc_int")
return True
def test_cl_int():
"""Generates random circle and line, computes the intersection,
and verifies that the number of intersections and the
positions of the intersections are correct.
Returns True or False"""
# 3 random points
p1 = vector.randvec(2, 0.0, 10.0,1.0)
p2 = vector.randvec(2, 0.0, 10.0,1.0)
p3 = vector.randvec(2, 0.0, 10.0,1.0)
# prevent div by zero / no line direction
if tol_eq(vector.norm(p1-p2),0):
p2 = p1 + p3 * 0.1
# line (o,v): origin p1, direction p1-p2
o = p1
v = (p2 - p1) / vector.norm(p2 - p1)
# cirlce (c, r): centered in p3, radius p3-p2 + rx
c = p3
r0 = vector.norm(p3-p2)
# cases rx = 0, rx > 0, rx < 0
case = random.choice([1,2,3])
if case==1:
r = r0 #should have one intersection (unles r0 = 0)
elif case==2:
r = random.random() * r0 # should have no ints (unless r0=0)
elif case==3:
r = r0 + random.random() * r0 # should have 2 ints (unless r0=0)
# do interesection
diag_print("problem:"+str(c)+","+str(r)+","+str(o)+","+str(v),"test_cl_int")
sols = cl_int(c,r,o,v)
diag_print("solutions:"+str(map(str, sols)),"test_cl_int")
# distance from point on line closest to circle center
l = vector.dot(c-o, v) / vector.norm(v)
p = o + v * l / vector.norm(v)
d = vector.norm(p-c)
diag_print("distance center to line="+str(d),"test_cl_int")
# test number of intersections
if len(sols) == 0:
if not tol_gt(d, r):
diag_print("wrong number of solutions: 0", "test_cl_int")
return False
elif len(sols) == 1:
if not tol_eq(d, r):
diag_print("wrong number of solutions: 1", "test_cl_int")
return False
elif len(sols) == 2:
if not tol_lt(d, r):
diag_print("wrong number of solutions: 2", "test_cl_int")
return False
else:
diag_print("wrong number of solutions: >2", "test_cl_int")
# test coordinates of intersection
for s in sols:
# s on line (o,v)
if not is_colinear(s, o, o+v):
diag_print("solution not on line", "test_cl_int")
return False
# s on circle c, r
if not tol_eq(vector.norm(s-c), r):
diag_print("solution not on circle", "test_cl_int")
return False
return True
def test_intersections():
sat = True
for i in range(0,100):
sat = sat and test_ll_int()
@ -619,7 +737,7 @@ def test1():
if sat:
print "ll_int() passed"
else:
print "ll_int() failed"
print "ll_int() FAILED"
sat = True
for i in range(0,100):
@ -627,27 +745,26 @@ def test1():
if not sat:
print "rr_int() failed"
return
if sat:
print "rr_int() passed"
else:
print "rr_int() failed"
print "rr_int() FAILED"
#sat = True
#for i in range(0,100):
# sat = sat and test_cc_int()
#if sat:
# print "cc_int() passed"
#else:
# print "cc_int() failed"
sat = True
for i in range(0,100):
sat = sat and test_cc_int()
if sat:
print "cc_int() passed"
else:
print "cc_int() FAILED"
#sat = True
#for i in range(0,100):
# sat = sat and test_cl_int()
#if sat:
# print "cl_int() passed"
#else:
# print "cl_int() failed"
sat = True
for i in range(0,100):
sat = sat and test_cl_int()
if sat:
print "cl_int() passed"
else:
print "cl_int() FAILED"
sat = True
for i in range(0,100):
@ -655,8 +772,10 @@ def test1():
if sat:
print "sss_int() passed"
else:
print "sss_int() failed"
print "sss_int() FAILED"
print "Note: did not test degenerate cases for sss_int"
def test_angles():
print "2D angles"
for i in xrange(9):
a = i * 45 * math.pi / 180
@ -701,7 +820,11 @@ def test_hcs_degen():
print make_hcs_3d(p1,p2,p3)
if __name__ == '__main__':
#test1()
#diag_select("test_cc_int")
#diag_select("test_cl_int")
#diag_select(".*")
test_intersections()
#test_angles()
#test_perp_3d()
#test_sss_degen()
test_hcs_degen()
#test_hcs_degen()

16
geosolver/matfunc.py
View File

@ -5,6 +5,7 @@
Version = 'File MATFUNC.PY, Ver 183, Date 12-Dec-2002,14:33:42'
# Modified by Rick van der Meiden 25-1-2007: included a makehash() and __hash__ in Table.
# Modified by Rick van der Meiden 20090311: added normSquared
# Modified by Rick van der Meiden 20101102: added Vec2 and Vec3 constructor functions
import operator, math, random
NPRE, NPOST = 0, 0 # Disables pre and post condition checks
@ -259,6 +260,8 @@ class LowerTri(Triangular):
x.append( (b[i] - x.dot(self[i][:i])) / self[i][i] )
return x
# ---------- class like constructor functions --------
def Mat( elems ):
'Factory function to create a new matrix.'
m, n = len(elems), len(elems[0])
@ -273,6 +276,13 @@ def Mat( elems ):
return Square(elems)
return LowerTri(elems)
def Vec2( x, y ):
return Vec([x,y])
def Vec3( x, y, z ):
return Vec([x,y,z])
# -------------- computing functions -----------
def funToVec( tgtfun, low=-1, high=1, steps=40, EqualSpacing=0 ):
'''Compute x,y points from evaluating a target function over an interval (low to high)
@ -298,6 +308,8 @@ def ratfit( (xvec, yvec), degree=2 ):
'Solves design matrix for approximating rational polynomial coefficients (a*x**2 + b*x + c)/(d*x**2 + e*x + 1)'
return Mat([[x**n for n in range(degree,-1,-1)]+[-y*x**n for n in range(degree,0,-1)] for x,y in zip(xvec,yvec)]).solve(Vec(yvec))
# ----------- matrix generators -------------
def genmat(m, n, func):
if not n: n=m
return Mat([ [func(i,j) for i in range(n)] for j in range(m) ])
@ -308,7 +320,7 @@ def zeroes(m=1, n=None):
def eye(m=1, n=None):
'Identity matrix with side length m-by-m or m-by-n'
return genmat(m,n, lambda i,j: float(i==j))
return genmat(m,n, lambda i,j: i==j)
def hilb(m=1, n=None):
'Hilbert matrix with side length m-by-m or m-by-n. Elem[i][j]=1/(i+j+1)'
@ -318,6 +330,8 @@ def rand(m=1, n=None):
'Random matrix with side length m-by-m or m-by-n'
return genmat(m,n, lambda i,j: random.random())
# ------------ unit test ----------
if __name__ == '__main__':
import cmath
a = Table([1+2j,2,3,4])

4
test/test.py
View File

@ -965,7 +965,9 @@ def test2d():
#test(ddd_problem())
#test(double_triangle())
#test(triple_double_triangle())
test(dad_problem())
#test(dad_problem())
#test(add_problem())
test(aad_problem())
if __name__ == "__main__":
test2d()