diff --git a/TODO.txt b/TODO.txt index 2159b91..685ed98 100644 --- a/TODO.txt +++ b/TODO.txt @@ -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 + 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 diff --git a/geosolver/clsolver2D.py b/geosolver/clsolver2D.py index 5b328da..1596b51 100644 --- a/geosolver/clsolver2D.py +++ b/geosolver/clsolver2D.py @@ -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: name of point variables + d: numeric distance values + a: 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)+")" + - diff --git a/geosolver/intersections.py b/geosolver/intersections.py index 51f7928..5fc248a 100644 --- a/geosolver/intersections.py +++ b/geosolver/intersections.py @@ -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() diff --git a/geosolver/matfunc.py b/geosolver/matfunc.py index 9ad7cbe..f4824b1 100644 --- a/geosolver/matfunc.py +++ b/geosolver/matfunc.py @@ -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]) diff --git a/test/test.py b/test/test.py index c58c32e..1b7a081 100644 --- a/test/test.py +++ b/test/test.py @@ -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()