class GeometricProblem (Notifier, Listener):
    """A geometric constraint problem with a prototpe.
    
       A problem consists of geometric variables (just variables for short), a prototype
       for each variable and constraints.

       Variables are of type Point, Line, etc. Alternatively, variables of any other hashable type
        (e.g. strings) are assumed to be points. (Depricated - used for backwards compatibility)
       
       Prototypes are of type vector 
       A point prototype must have length equal to the dimensionality as the problem (D).
       A line prototype must have length 2*D: it represents two points though which the line passes 
       A plane prototype must have length 3*D: it represents three points though which the plane passes 
       
       Supported constraints are instances of ParametricConstraint, FixConstraint, SelectionConstraint, etc.
       
       GeometricProblem listens for changes in constraint parameters and passes
       these changes, and changes in the system of constraints and the prototype, 
       to any other listerers (e.g. GeometricSolver) 
    """
    
    def __init__(self, dimension, use_prototype=True):
        """Initialize a new problem. Must specify dimensionality of the problem (i.e. dimension of points) and 
            wheter to use the prototype for solution selection."""
        Notifier.__init__(self)
        Listener.__init__(self)
        self.dimension = dimension      # dimensionality of points
        self.prototype = {}             # mapping from variables to prototypes
        self.cg = ConstraintGraph()     # constraint graph
        self.use_prototype = use_prototype;     # whether to use prototype for solution selection

    # ----------- prototype --------

    def set_prototype_selection(self, enabled):
        """Enable (True, the default) or disable (False) use of prototype for solution selection"""
        self.use_prototype = enabled
        self.dr.set_prototype_selection(self.use_prototype)
        self.send_notify(("set_prototype_selection", self.use_prototype))

    def get_prototype_selection(self):
        """Return True if prototype selection has been enabled (the default) or False otherwise"""
        return self.use_prototype

    # ---------------- variables -----------

    def add_variable(self, variable, prototype):
        """add a variable with a prototype"""
        prototypevector = vector.vector(prototype)
        # check dimension of prototype
        if isinstance(variable, Point):
            assert len(prototypevector) == self.dimension
        elif isinstance(variable, Line):
            assert len(prototypevector) == 2 * self.dimension
        elif isinstance(variable, Plane):
            assert len(prototypevector) == 3 * self.dimension
        else:
            # assume point
            assert len(prototypevector) == self.dimension
        if variable not in self.prototype:
            self.prototype[variable] = prototypevector
            self.cg.add_variable(variable)
        else:
            raise StandardError, "variable already in problem"

    def has_variable(self, variable):
        """returns True if variable in problem"""
        return variable in self.prototype

    def rem_variable(self, variable):
        """remove a variable (and all constraints incident imposed on it)"""
        if variable in self.prototype:
            del self.prototype[variable] 
            self.cg.rem_variable(variable) 

    def set_prototype(self, variable, prototype):
        """set prototype of variable"""
        prototypevector = vector.vector(prototype)
        if variable in self.prototype:
            self.prototype[variable] = prototypevector
            self.send_notify(("set_prototype", (variable,prototypevector)))
        else:
            raise StandardError, "unknown variable variable"

    def get_prototype(self, variable):
        """get prototype of variable"""
        if variable in self.prototype:
            return self.prototype[variable]
        else:
            raise StandardError, "unknown variable "+str(variable)

    # --------------- points - depricated  -------------
    def add_point(self, variable, prototype):
        """depricated - use add_variable"""
        return self.add_variable(variable, prototype)

    def has_point(self, variable):
        """depricated - use has_variable""" 
        return self.has_variable(variable)

    def set_point(self, variable, prototype):
        """depricated - use set_prototype"""
        return self.set_prototype(variable, prototype)

    def get_point(self, variable):
        """depricated - use get_prototype"""
        return self.get_prototype(variable)
    
    # ----------- constraints --------
 
    def add_constraint(self, con):
        """add a constraint"""
        # check that variables in problem
        for var in con.variables():
                if var not in self.prototype:
                    raise StandardError, "variable %s not in problem"%(var)
        # check that constraint not already in problem
        if isinstance(con, DistanceConstraint):
            if self.get_distance(con.variables()[0],con.variables()[1]):
                raise StandardError, "distance already in problem"
        elif isinstance(con, AngleConstraint):
            if self.get_angle(con.variables()[0],con.variables()[1], con.variables()[2]):
                raise StandardError, "angle already in problem"
        elif isinstance(con, RigidConstraint):
            if self.get_rigid(con.variables()):
                raise StandardError, "rigid already in problem"
        elif isinstance(con, MateConstraint):
            if self.get_mate(con.variables()):
                raise StandardError, "mate constraint already in problem"
        elif isinstance(con, SelectionConstraint):
            pass
        elif isinstance(con, FixConstraint):
            if self.get_fix(con.variables()[0]):
                raise StandardError, "fix already in problem"
        elif isinstance(con, CoincidenceConstraint):
            if self.get_coincidence(con.variables()[0], con.variables()[1]):
                raise StandardError, "coincidence already in problem"
        else:
            raise StandardError, "unsupported constraint type"
        # passed tests, add to poblem
        if isinstance(self, ParametricConstraint):
            con.add_listener(self)
        self.cg.add_constraint(con)

    def get_distance(self, a, b):
        """return the distance constraint on given points, or None"""
        on_a = self.cg.get_constraints_on(a)
        on_b = self.cg.get_constraints_on(b)
        on_ab = filter(lambda c: c in on_a and c in on_b, on_a)
        distances = filter(lambda c: isinstance(c, DistanceConstraint), on_ab)
        if len(distances) > 1:
            raise StandardError, "multiple constraints found"
        elif len(distances) == 1:
            return distances[0]
        else:
            return None
 
    def get_angle(self, a, b, c):
        """return the angle constraint on given points, or None"""
        on_a = self.cg.get_constraints_on(a)
        on_b = self.cg.get_constraints_on(b)
        on_c = self.cg.get_constraints_on(c)
        on_abc = filter(lambda x: x in on_a and x in on_b and x in on_c, on_a)
        angles = filter(lambda x: isinstance(x, AngleConstraint), on_abc)
        candidates = filter(lambda x: x.variables()[1] == b, angles)
        if len(candidates) > 1:
            raise StandardError, "multiple constraints found"
        elif len(candidates) == 1:
            return candidates[0]
        else:
            return None

    def get_fix(self, p):
        """return the fix constraint on given point, or None"""
        on_p = self.cg.get_constraints_on(p)
        fixes = filter(lambda x: isinstance(x, FixConstraint), on_p)
        if len(fixes) > 1:
            raise StandardError, "multiple constraints found"
        elif len(fixes) == 1:
            return fixes[0]
        else:
            return None

    def get_constraints_with_type_on_variables(self, constrainttype, variables):
        candidates = None
        for var in variables:
            if candidates == None:
                candidates = set(filter(lambda c:isinstance(c,constrainttype),self.cg.get_constraints_on(var)))
            else:    
                candidates.intersection_update(filter(lambda c:isinstance(c,constrainttype),self.cg.get_constraints_on(var)))
        return candidates

    def get_unique_constraint(self, constrainttype, variables):
        candidates = self.get_constraints_with_type_on_variables(constrainttype, variables)
        if len(candidates) == 0:
            return None
        elif len(candidates) == 1:
            return list(candidates)[0]
        else: # >= 1
            raise StandardError, "multiple constraints found"
        
    def get_coincidence(self, p, g):
        return self.get_unique_constraint(CoincidenceConstraint, [p,g])

    def get_rigid(self, variables):
        return self.get_unique_constraint(RigidConstraint, variables)

    def get_mate(self, variables):
        return self.get_unique_constraint(MateConstraint, variables)

    def get_coincident_points(self, geometry):
        coincidences = self.get_constraints_with_type_on_variables(CoincidenceConstraint, [geometry])
        points = set()
        for constraint in coincidences:
            points.update(filter(lambda var: isinstance(var, Point) and var != geometry, constraint.variables()))
        return points

    def verify(self, solution):
        """returns true iff all constraints satisfied by given solution. 
           solution is a dictionary mapping variables (names) to values (points)"""
        if solution == None:
            sat = False
        else:
            sat = True
            for con in self.cg.constraints():
                solved = True
                for v in con.variables():
                    if v not in solution:
                        solved = False
                        break
                if not solved:
                    diag_print(str(con)+" not solved", "GeometricProblem.verify")
                    sat = False
                elif not con.satisfied(solution):
                    diag_print(str(con)+" not satisfied", "GeometricProblem.verify")
                    sat = False
        return sat
       
    def rem_point(self, var):
        """remove a point variable from the constraint system"""
        if var in self.prototype:
            self.cg.rem_variable(var)
            del self.prototype[var]
        else:
            raise StandardError, "variable "+str(var)+" not in problem."
            
    def rem_constraint(self, con):
        """remove a constraint from the constraint system"""
        if con in self.cg.constraints():
            if isinstance(con, SelectionConstraint): 
                self.send_notify(("rem_selection_constraint", con))
            self.cg.rem_constraint(con)
        else:
            raise StandardError, "no constraint "+str(con)+" in problem."

    def receive_notify(self, object, notify):
        """When notified of changed constraint parameters, pass on to listeners"""
        if isinstance(object, ParametricConstraint):
            (message, data) = notify
            if message == "set_parameter":
                self.send_notify(("set_parameter",(object,data)))
        #elif object == self.cg:
        #    self.send_notify(notify)

    def __str__(self):
        s = ""
        for v in self.prototype:
            s += str(v) + " = " + str(self.prototype[v]) + "\n"
        for con in self.cg.constraints():
            s += str(con) + "\n"
        s+= "prototype-based selection = " + str(self.use_prototype)
        return s class GeometricSolver (Listener):
    """The GeometricSolver monitors changes in a GeometricProblem and 
       maps any changes to corresponding changes in a GeometricDecomposition
    """

    # public methods
    def __init__(self, problem):
        """Create a new GeometricSolver instance
        
           keyword args
            problem        - the GeometricProblem instance to be monitored for changes
        """
        # init superclasses
        Listener.__init__(self)

        # init variables
        # the problem on which this solver works 
        self.problem = problem
        # shortcut to the dimension of the problem
        self.dimension = problem.dimension
        # shortcut to the constraint graph of the problem
        self.cg = problem.cg  
        # the cluster-solver (or dr for decomposition-recombination) 
        if self.problem.dimension == 2:
            self.dr = ClusterSolver2D()
        elif self.problem.dimension == 3:
            self.dr = ClusterSolver3D()
        else:
            raise StandardError, "Sorry, can't solve problems of dimension < 2 or > 3."
        
        # a map from problem variables and constraints to clusters, and vice versa
        self._map = {}
        
        # enable prototype based selection by default
        self._set_prototype_selection(problem.get_prototype_selection())

        # register 
        self.cg.add_listener(self)
        self.dr.add_listener(self)

        # create an initial fix cluster
        self.fixvars = []
        self.fixcluster = None

        # add variables
        for var in self.cg.variables():
            self._add_variable(var)
       
        # add constraints
        toadd = set(self.cg.constraints())

        # add coincidences first. Prevents re-mapping of primitves and re-solving of problem 
        for con in list(toadd):
            if isinstance(con, CoincidenceConstraint): 
                self._add_constraint(con)
                toadd.remove(con)

        # add selection constraints first. Prevents re-evaluation 
        for con in list(toadd):
            if isinstance(con, SelectionConstraint): 
                self._add_constraint(con)
                toadd.remove(con)

        # add distances first. Nicer decomposition in Rigids
        for con in list(toadd):
            if isinstance(con, DistanceConstraint): 
                self._add_constraint(con)
                toadd.remove(con)

        # add other constraints. 
        for con in toadd:
            self._add_constraint(con)

    def get_constrainedness(self):
        """Depricated. Use get_status instead"""
        return self.get_status()

    def get_result(self):
        """Depricated. Use get_decomposition."""
        return self.get_decomposition()  

    def get_cluster(self):
        """Depricated. Use get_decomposition."""
        return self.get_decomposition()  
    
    def get_decomposition(self):
        """Returns a GeometricDecomposition (the root of a tree of clusters),
         describing the solutions and the decomposition of the problem."""
        # several drcluster can maps to a single geoclusters 
        map = {}   
        geoclusters = []
        # map dr clusters
        for drcluster in filter(lambda c: isinstance(c, Rigid), self.dr.clusters()):
            # create geocluster and map to drcluster (and vice versa)
            geocluster = GeometricDecomposition(drcluster.vars)
            if geocluster not in map:
                map[drcluster] = geocluster
                map[geocluster] = [drcluster]
                geoclusters.append(geocluster)
            else:
                geocluster = map[map[geocluster][0]]
                map[drcluster] = geocluster
                map[geocluster].append(drcluster)
            
        for geocluster in geoclusters:
            # pick newest drcluster 
            drclusters = map[geocluster]
            drcluster = max(drclusters, key=lambda c: c.creationtime)
            # determine solutions
            geocluster.solutions = self._map_cluster_solutions(drcluster)
            # determine incidental underconstrainedness
            underconstrained = False
            configurations = self.dr.get(drcluster)
            if configurations != None:
                for config in configurations:
                    if config.underconstrained:
                        underconstrained = True
            # determine flag
            if drcluster.overconstrained:
                geocluster.flag = GeometricDecomposition.OVERCONSTRAINED
            elif geocluster.solutions == None:
                geocluster.flag = GeometricDecomposition.UNSOLVED
            elif len(geocluster.solutions) == 0:
                geocluster.flag = GeometricDecomposition.INCONSISTENT
            elif underconstrained:
                geocluster.flag = GeometricDecomposition.DEGENERTE
            else:
                geocluster.flag = GeometricDecomposition.OK


        # determine subclusters
        for method in self.dr.methods():
            if not isinstance(method, PrototypeMethod) and not isinstance(method, SelectionMethod):
                for out in method.outputs():
                    if isinstance(out, Rigid):
                        parent = map[out]
                        for inp in method.inputs():
                            if isinstance(inp, Rigid):
                                sub = map[inp]
                                if sub != parent and sub not in parent.subs:
                                    parent.subs.append(sub)
        
        # determine result from top-level clusters 
        top = self.dr.top_level()
        rigids = filter(lambda c: isinstance(c, Rigid), top)
        if len(top) > 1:
            # structurally underconstrained cluster
            result = GeometricDecomposition(self.problem.cg.variables())
            result.flag = GeometricDecomposition.UNDERCONSTRAINED
            for cluster in rigids:
                result.subs.append(map[cluster])
        else:
            if len(rigids) == 1:
                # structurally well constrained, or structurally overconstrained
                result = map[rigids[0]]
            else:
                # no variables in problem?
                result = GeometricDecomposition(self.problem.cg.variables())
                result.variables = []
                result.subs = []
                result.solutions = []
                result.flags = GeometricDecomposition.UNSOLVED
        return result 

    
    def get_solutions(self):
        """Returns a list of Configurations, which will be empty if the
           problem has no solutions. Note: this method is
           cheaper but less informative than get_decomposition. 
        """
        #"""The list and the configurations should not be changed (since they are
        #references to objects in the solver)."""
        # find top level rigid and all its configurations
        rigids = filter(lambda c: isinstance(c, Rigid), self.dr.top_level())
        if len(rigids) != 0:   
            solutions = self._map_cluster_solutions(drcluster)
        else:
            solutions = []
        return solutions

    def _map_cluster_solutions(self, drcluster):
        # map dr-cluster configurations to solutions, i.e. a map from problem variables to values           
        configurations = self.dr.get(drcluster)
        solutions = []
        diag_print("mapping cluster "+str(drcluster)+" #configurations="+str(len(configurations)),"GeometricSolver")
        for configuration in configurations:
            solution = {}
            for var in self.problem.cg.variables():
                if isinstance(var, Point):
                    assert len(self._map[var].vars) == 1
                    point = iter(self._map[var].vars).next()
                    if point in configuration:
                        solution[var] = configuration[point]
                elif isinstance(var, Line):
                    line_rigid = self._map[var]
                    if self.dimension == 2:
                        line_vertex = line_rigid.vertex 
                        line_normal = line_rigid.normal
                        if line_vertex in configuration and line_normal in configuration:
                            p1 = configuration[line_vertex]
                            n = configuration[line_normal]
                            p2 = p1 + perp_2d(n-p1)
                            solution[var] = p1.concatonated(p2)
                        #endif in config
                    elif self.dimension == 3:
                        line_vertex = line_rigid.vertex 
                        line_normal1 = line_rigid.normal1
                        line_normal2 = line_rigid.normal2
                        if line_vertex in configuration and line_normal1 in configuration and line_normal2 in configuration:
                            p1 = configuration[line_vertex]
                            n1 = configuration[line_normal1]
                            n2 = configuration[line_normal2]
                            p2 = p1 + vector.cross(n1-p1, n2-p1)
                            solution[var] = p1.concatonated(p2)
                        #endif in config
                    #endif dimension 
                else:
                    # assume point - depricated
                    assert len(self._map[var].vars) == 1
                    point = iter(self._map[var].vars).next()
                    if point in configuration:
                        solution[var] = configuration[point]
            #for
            solutions.append(solution)
        #for    
        return solutions

    def get_status(self):
        """Returns a symbolic flag, one of:
            GeometricDecomposition.UNDERCONSTRAINED, 
            GeometricDecomposition.OVERCONSTRAINED,
            GeometricDecomposition.OK,
            GeometricDecomposition.UNSOLVED,
            GeometricDecomposition.EMPTY,
            GeometricDecomposition.INCONSISTENT,
            GeometricDecomposition.DEGENERTE.
           Note: this method is cheaper but less informative than get_decomposition. 
        """
        rigids = filter(lambda c: isinstance(c, Rigid), self.dr.top_level())
        if len(rigids) == 0:            
            return GeometricDecomposition.EMPTY
        elif len(rigids) == 1:
            drcluster = rigids[0]
            solutions = self.dr.get(drcluster)
            underconstrained = False
            if solutions == None:
                return GeometricDecomposition.UNSOLVED
            else:
                for solution in solutions:
                    if solution.underconstrained:
                        underconstrained = True
            if drcluster.overconstrained:
                return GeometricDecomposition.OVERCONSTRAINED
            elif len(solutions) == 0:
                return GeometricDecomposition.INCONSISTENT
            elif underconstrained:
                return GeometricDecomposition.DEGENERTE
            else:
                return GeometricDecomposition.OK
        else:
            return GeometricDecomposition.UNDERCONSTRAINED
   

    def receive_notify(self, object, message):
        """Take notice of changes in constraint problem"""
        if object == self.cg:
            (type, data) = message
            if type == "add_constraint":
                self._add_constraint(data)
            elif type == "rem_constraint":
                self._rem_constraint(data)
            elif type == "add_variable":
                self._add_variable(data)
            elif type == "rem_variable":
                self._rem_variable(data)
            elif type == "set_prototype_selection":
                self.set_prototype_selection(data)
            else:
                raise StandardError, "unknown message type"+str(type)
        elif object == self.problem:
            (type, data) = message
            if type == "set_point":
                (variable, point) = data
                self._update_variable(variable) 
            elif type == "set_parameter":
                (constraint, value) = data
                self._update_constraint(constraint)
            else:
                raise StandardError, "unknown message type"+str(type)
        elif object == self.dr:
            pass
        else:
            raise StandardError, "message from unknown source"+str((object, message))
    
    # --------------- internal methods ------------------

    def _set_prototype_selection(self, enabled):
        """Enable (True) or disable (False) use of prototype for solution selection"""
        self.dr.set_prototype_selection(enabled)

    def _add_variable(self, var):
        if isinstance(var, Point):
            self._add_point(var)
        elif isinstance(var, Line):
            self._add_line(var)
        elif isinstance(var, Plane):
            self._add_plane(var)
        else:
            # assume point - depricated
            self._add_point(var)

    def _add_point(self, var):
        if var not in self._map:
            rigid = Rigid([var])
            self._map[var] = rigid
            self._map[rigid] = var
            self.dr.add(rigid)
            self._update_variable(var)
    
    def _add_line(self, var):
        diag_print("add line "+str(var),"GeometricSolver")
        # find coincident points
        points = list(self.problem.get_coincident_points(var))
        diag_print("on "+str(points),"GeometricSolver")
        if self.dimension == 2:
            if len(points) == 0:
                self._map_line_distance(var)
            elif len(points) >= 1:
                self._map_line_point_distance(var, points[0])
        elif self.dimension == 3:
            if len(points) == 0:
                self._map_line_3d_distance(var)
            elif len(points) >= 1:
                self._map_line_3d_point_distance(var, points[0])

    def _map_line_distance(self,line):
        # map a line (coincident with no points) to a distance cluster (on two new point variables)
        v = str(line)+"_vertex"
        n = str(line)+"_normal"
        dist = Rigid([v,n])
        # add add-hoc attributes to rigid, so we can distinguish vertex and normal! 
        dist.vertex = v
        dist.normal = n
        # add rigids for created points, needed for prototypes
        # NOTE: adding non-problem variables to mapping! 
        # TODO: clean up after removal of line
        vertex_rigid = Rigid([dist.vertex])
        self.dr.add(vertex_rigid)
        self._map[dist.vertex] = vertex_rigid 
        normal_rigid = Rigid([dist.normal])
        self.dr.add(normal_rigid)
        self._map[dist.normal] = normal_rigid 
        # add line to mapping
        self._map[line] = dist
        self._map[dist] = line
        self.dr.add(dist)
        diag_print("mapped "+str(line)+" to "+str(dist),"GeometricSolver")   
        # update configurations
        self._update_variable(line)
 
    def _map_line_point_distance(self,line, point):
        # map a line coincident with one point to a distance clusters (and one new point variable)
        v = list(self._map[point].vars)[0]
        n = str(line)+"_normal"
        dist = Rigid([v,n])
        # add add-hoc attributes to rigid, so we can distinguish vertex and normal! 
        dist.vertex = v
        dist.normal = n 
        # add rigids for created points, needed for prototypes
        # NOTE: adding non-problem variables to mapping! 
        # TODO: clean up after removal of line
        normal_rigid = Rigid([dist.normal])
        self.dr.add(normal_rigid)
        self._map[dist.normal] = normal_rigid 
        # add to mapping 
        self._map[line] = dist
        self._map[dist] = line
        self.dr.add(dist)
        diag_print("mapped "+str(line)+" to "+str(dist),"GeometricSolver")   
        self._update_variable(line)
 
    def _map_line_3d_distance(self,line):
        # map a line (coincident with no points) to a distance cluster (on two new point variables)
        v = str(line)+"_vertex"
        n1 = str(line)+"_normal1"
        n2 = str(line)+"_normal2"
        dist = Rigid([v,n1, n2])
        # add add-hoc attributes to rigid, so we can distinguish vertex and normal! 
        dist.vertex = v
        dist.normal1 = n1
        dist.normal2 = n2
        # add rigids for created points, needed for prototypes
        # NOTE: adding non-problem variables to mapping! 
        # TODO: clean up after removal of line
        vertex_rigid = Rigid([dist.vertex])
        self.dr.add(vertex_rigid)
        self._map[dist.vertex] = vertex_rigid 
        normal1_rigid = Rigid([dist.normal1])
        self.dr.add(normal1_rigid)
        self._map[dist.normal1] = normal1_rigid 
        normal2_rigid = Rigid([dist.normal2])
        self.dr.add(normal2_rigid)
        self._map[dist.normal2] = normal2_rigid 
        # add line to mapping
        self._map[line] = dist
        self._map[dist] = line
        self.dr.add(dist)
        diag_print("mapped "+str(line)+" to "+str(dist),"GeometricSolver")   
        # update configurations
        self._update_variable(line)

    def _map_line_3d_point_distance(self,line, point):
        # map a line coincident with one point to a distance clusters (and one new point variable)
        v = list(self._map[point].vars)[0]
        n1 = str(line)+"_normal1"
        n2 = str(line)+"_normal2"
        dist = Rigid([v,n1, n2])
        # add add-hoc attributes to rigid, so we can distinguish vertex and normal! 
        dist.vertex = v
        dist.normal1 = n1 
        dist.normal2 = n2 
        # add rigids for created points, needed for prototypes
        # NOTE: adding non-problem variables to mapping! 
        # TODO: clean up after removal of line
        normal1_rigid = Rigid([dist.normal1])
        self.dr.add(normal1_rigid)
        self._map[dist.normal1] = normal1_rigid 
        normal2_rigid = Rigid([dist.normal2])
        self.dr.add(normal2_rigid)
        self._map[dist.normal2] = normal2_rigid 
        # add to mapping 
        self._map[line] = dist
        self._map[dist] = line
        self.dr.add(dist)
        diag_print("mapped "+str(line)+" to "+str(dist),"GeometricSolver")   
        self._update_variable(line)

    def _rem_variable(self, var):
        diag_print("GeometricSolver._rem_variable","GeometricSolver")
        if var in self._map:
            self.dr.remove(self._map[var])      
            # Note: CLSolver automatically removes variables with no dependent clusters
            del self._map[var]

    def _add_constraint(self, con):
        if isinstance(con, AngleConstraint):
            # map to hedgdehog
            vars = list(con.variables());
            hog = Hedgehog(vars[1],[vars[0],vars[2]])
            self._map[con] = hog
            self._map[hog] = con
            self.dr.add(hog)
            # set configuration
            self._update_constraint(con)
        elif isinstance(con, DistanceConstraint):
            # map to rigid
            vars = list(con.variables());
            rig = Rigid([vars[0],vars[1]])
            self._map[con] = rig
            self._map[rig] = con
            self.dr.add(rig)
            # set configuration
            self._update_constraint(con)
        elif isinstance(con, RigidConstraint):
            # map to rigid
            vars = list(con.variables());
            rig = Rigid(vars)
            self._map[con] = rig
            self._map[rig] = con
            self.dr.add(rig)
            # set configuration
            self._update_constraint(con)
        elif isinstance(con, MateConstraint):
            # map to glueable cluster
            vars = list(con.variables());
            glue = Glueable(vars)
            self._map[con] = glue
            self._map[glue] = con
            self.dr.add(glue)
            # set configuration
            self._update_constraint(con)
        elif isinstance(con, FixConstraint):
            if self.fixcluster != None:
                self.dr.remove(self.fixcluster)
            self.fixvars.append(con.variables()[0])
            if len(self.fixvars) >= 1:
            #if len(self.fixvars) >= self.problem.dimension:
                self.fixcluster = Rigid(self.fixvars)
                self.dr.add(self.fixcluster)
                self.dr.set_root(self.fixcluster)
                self._update_fix()
        elif isinstance(con, SelectionConstraint):
            # add directly to clustersolver
            self.dr.add_selection_constraint(con)
        elif isinstance(con, CoincidenceConstraint):
            # re-map lines, planes, etc
            lines = filter(lambda var: isinstance(var,Line),con.variables())
            points = filter(lambda var: isinstance(var,Point),con.variables())
            if len(lines)==1 and len(points)==1:
                line = iter(lines).next()
                point = iter(points).next()
                # re-map line if needed
                #self._rem_variable(line)
                #self._add_line(line)
                # map coincience constraint of a point with a line
                line_rigid  = self._map[line]
                point_rigid = self._map[point]
                point_vertex = iter(point_rigid.vars).next()
                if point_vertex not in line_rigid.vars:
                    if self.dimension==2:
                        line_vertex = line_rigid.vertex
                        line_normal = line_rigid.normal  
                        angle_hog = Hedgehog(line_vertex,[line_normal, point_vertex])
                        self._map[con] = angle_hog 
                        self._map[angle_hog] = con
                        self.dr.add(angle_hog)
                        diag_print("mapped "+str(con)+" to "+str(angle_hog),"GeometricSolver")   
                        self._update_constraint(con)
                    elif self.dimension==3:
                        line_vertex = line_rigid.vertex
                        line_normal1 = line_rigid.normal1  
                        line_normal2 = line_rigid.normal2  
                        angle_hog = Hedgehog(line_vertex,[line_normal1, line_normal2, point_vertex])
                        self._map[con] = angle_hog 
                        self._map[angle_hog] = con
                        self.dr.add(angle_hog)
                        diag_print("mapped "+str(con)+" to "+str(angle_hog),"GeometricSolver")   
                        self._update_constraint(con)
                    #endif dimension
                #endif point_vertex
            #endif co(line,point)
        else:
            raise StandardError, "unknown constraint type"
            pass
         
    def _rem_constraint(self, con):
        diag_print("GeometricSolver._rem_constraint","GeometricSolver")
        if isinstance(con,FixConstraint):
            if self.fixcluster != None:
                self.dr.remove(self.fixcluster)
            var = con.variables()[0]
            if var in self.fixvars:
                self.fixvars.remove(var)
            #if len(self.fixvars) < self.problem.dimension:
            if len(self.fixvars) == 0:
                self.fixcluster = None
            else:
                self.fixcluster = Rigid(self.fixvars)
                self.dr.add(self.fixcluster)
                self.dr.set_root(self.fixcluster)
        elif isinstance(con, SelectionConstraint):
            # remove directly from clustersolver
            self.dr.rem_selection_constraint(con)
        elif con in self._map:
            self.dr.remove(self._map[con])
            del self._map[con]
  
    # update methods: set the value of the variables in the constraint graph

    def _update_constraint(self, con):
        if isinstance(con, AngleConstraint):
            # set configuration
            hog = self._map[con]
            vars = list(con.variables())
            v0 = vars[0]
            v1 = vars[1]
            v2 = vars[2]
            angle = con.get_parameter()
            p0 = vector.vector([1.0,0.0])
            p1 = vector.vector([0.0,0.0])
            p2 = vector.vector([math.cos(angle), math.sin(angle)])
            # pad vectors to right dimension
            if self.dimension == 3:
                p0.append(0.0)
                p1.append(0.0)
                p2.append(0.0)
            conf = Configuration({v0:p0,v1:p1,v2:p2})
            self.dr.set(hog, [conf])
            assert con.satisfied(conf.map)
        elif isinstance(con, DistanceConstraint):
            # set configuration
            rig = self._map[con]
            vars = list(con.variables())
            v0 = vars[0]
            v1 = vars[1]
            dist = con.get_parameter()
            #p0 = vector.vector([0.0,0.0])
            #p1 = vector.vector([dist,0.0])
            # use prototype to orient rigid - minimize difference solution and prototype
            p0 = self.problem.get_prototype(v0)
            v = self.problem.get_prototype(v1) - p0
            if vector.norm(v) != 0:
                v = v / vector.norm(v)
            else:
                v = vector.vector([0.0 for i in range(self.dimension)])
                v[0] = 1.0
            p1 = p0+v*dist
            conf = Configuration({v0:p0,v1:p1})
            self.dr.set(rig, [conf])
            assert con.satisfied(conf.map)
        elif isinstance(con, MateConstraint):
            # set configuration
            if self.dimension != 3:
                raise Exception, "MateConstraint only supported in 3D"
            glue = self._map[con]
            vars = list(con.variables())
            vo1 = vars[0]
            vx1 = vars[1]
            vy1 = vars[2]
            vo2 = vars[3]
            vx2 = vars[4]
            vy2 = vars[5]
            po1 = vector.vector([0.0,0.0,0.0])
            px1 = vector.vector([1.0,0.0,0.0])
            py1 = vector.vector([0.0,1.0,0.0])
            trans = con.get_parameter()
            po2 = transform_point(po1, trans)
            px2 = transform_point(px1, trans)
            py2 = transform_point(py1, trans)
            conf = Configuration({
                vo1:po1,
                vx1:px1,
                vy1:py1,
                vo2:po2,
                vx2:px2,
                vy2:py2,
            })
            self.dr.set(glue, [conf])
            assert con.satisfied(conf.map)
        elif isinstance(con, RigidConstraint):
            # set configuration
            rig = self._map[con]
            vars = list(con.variables())
            conf = con.get_parameter()
            self.dr.set(rig, [conf])
            assert con.satisfied(conf.map)
        elif isinstance(con, FixConstraint):
            self._update_fix()
        elif isinstance(con, CoincidenceConstraint):
            lines = filter(lambda var: isinstance(var,Line),con.variables())
            points = filter(lambda var: isinstance(var,Point),con.variables())
            if len(lines)==1 and len(points)==1:
                line = iter(lines).next()
                point = iter(points).next()
                if self.dimension == 2:
                    line_rigid = self._map[line]
                    point_rigid = self._map[point]
                    point_vertex = iter(point_rigid.vars).next()
                    print "point_vertex", point_vertex
                    line_vertex = line_rigid.vertex
                    line_normal = line_rigid.normal
                    angle_hog = self._map[con]
                    pv = vector.vector([1.0,0.0])
                    lv = vector.vector([0.0,0.0])
                    ln = vector.vector([0.0,1.0])
                    conf1 = Configuration({line_vertex:lv, line_normal:ln, point_vertex: 1.0*pv})
                    conf2 = Configuration({line_vertex:lv, line_normal:ln, point_vertex:-1.0*pv})
                    self.dr.set(angle_hog, [conf1,conf2])
                    diag_print("set "+str(angle_hog)+" to "+str(conf1),"GeometricSolver")
                    diag_print("set "+str(angle_hog)+" to "+str(conf2),"GeometricSolver")
                elif self.dimension == 3:
                    line_rigid = self._map[line]
                    point_rigid = self._map[point]
                    point_vertex = iter(point_rigid.vars).next()
                    print "point_vertex", point_vertex
                    line_vertex = line_rigid.vertex
                    line_normal1 = line_rigid.normal1
                    line_normal2 = line_rigid.normal2
                    angle_hog = self._map[con]
                    lv = vector.vector([0.0,0.0,0.0])
                    pv = vector.vector([1.0,0.0,0.0])
                    ln1 = vector.vector([0.0,1.0,0.0])
                    ln2 = vector.vector([0.0,0.0,1.0])
                    conf1 = Configuration({line_vertex:lv, line_normal1:ln1, line_normal2:ln2, point_vertex: 1.0*pv})
                    conf2 = Configuration({line_vertex:lv, line_normal1:ln1, line_normal2:ln2, point_vertex:-1.0*pv})
                    self.dr.set(angle_hog, [conf1,conf2])
                    diag_print("set "+str(angle_hog)+" to "+str(conf1),"GeometricSolver")
                    diag_print("set "+str(angle_hog)+" to "+str(conf2),"GeometricSolver")
               #endif dimension
        else:
            raise StandardError, "unknown constraint type"
    
    def _update_variable(self, var):
        if isinstance(var, Point):
            self._update_point(var)
        elif isinstance(var, Line):
            self._update_line(var)
        elif isinstance(var, Plane):
            self._update_plane(var)
        else:
            # assume point - depricated
            self._update_point(var)

    def _update_point(self, variable):
        cluster = self._map[variable]
        proto = self.problem.get_prototype(variable)
        conf = Configuration({variable:proto})
        self.dr.set(cluster, [conf])

    def _update_line(self, variable):
        cluster = self._map[variable]
        proto = self.problem.get_prototype(variable)
        if self.dimension == 2:
            line_vertex = cluster.vertex
            line_normal = cluster.normal
            # determine vertex and normal prototype coordinates
            p1 = proto[0:2]
            p2 = proto[2:4]
            v = p1
            n = perp_2d(p2-p1)
            # update prototypes of created point variables
            if line_vertex in self._map:
                vertex_rigid = self._map[line_vertex]
                conf = Configuration({line_vertex: v})
                self.dr.set(vertex_rigid, [conf])
                diag_print("set "+str(vertex_rigid)+" to "+str(conf),"GeometricSolver")
            if line_normal in self._map:
                normal_rigid = self._map[line_normal]
                conf = Configuration({line_normal: n})
                self.dr.set(normal_rigid, [conf])
                diag_print("set "+str(normal_rigid)+" to "+str(conf),"GeometricSolver")
            # update line configuration
            conf = Configuration({line_vertex:v, line_normal:n})
            self.dr.set(cluster, [conf])
            diag_print("set "+str(cluster)+" to "+str(conf),"GeometricSolver")
        elif self.dimension == 3:
            line_vertex = cluster.vertex
            line_normal1 = cluster.normal1
            line_normal2 = cluster.normal2
            # determine vertex and normal prototype coordinates
            p1 = proto[0:3]
            p2 = proto[3:6]
            v = p1
            n0 = p2
            n1,n2 = perp_3d(p2-p1)
            n1 = n1 + p1
            n2 = n2 + p1
            # update prototypes of created point variables
            if line_vertex in self._map:
                vertex_rigid = self._map[line_vertex]
                conf = Configuration({line_vertex: v})
                self.dr.set(vertex_rigid, [conf])
                diag_print("set "+str(vertex_rigid)+" to "+str(conf),"GeometricSolver")
            if line_normal1 in self._map:
                normal1_rigid = self._map[line_normal1]
                conf = Configuration({line_normal1: n1})
                self.dr.set(normal1_rigid, [conf])
                diag_print("set "+str(normal1_rigid)+" to "+str(conf),"GeometricSolver")
            if line_normal2 in self._map:
                normal2_rigid = self._map[line_normal2]
                conf = Configuration({line_normal2: n2})
                self.dr.set(normal2_rigid, [conf])
                diag_print("set "+str(normal2_rigid)+" to "+str(conf),"GeometricSolver")
            # update line configuration
            conf = Configuration({line_vertex:v, line_normal1:n1, line_normal2:n2})
            self.dr.set(cluster, [conf])
            diag_print("set "+str(cluster)+" to "+str(conf),"GeometricSolver")
        #endif dimension
    #fed _update_line

    def _update_fix(self):
        if self.fixcluster:
            vars = self.fixcluster.vars
            map = {}
            for var in vars:
                map[var] = self.problem.get_fix(var).get_parameter()
            conf = Configuration(map)
            self.dr.set(self.fixcluster, [conf])
        else:
            diag_print("no fixcluster to update","geometric")
            pass class DistanceConstraint(ParametricConstraint):
    """A constraint on the Euclidean distance between two points"""
    
    def __init__(self, a, b, dist):
        """Create a new DistanceConstraint instance
        
           keyword args:
            a    - a point variable name 
            b    - a point variable name
            dist - the distance parameter value
        """
        ParametricConstraint.__init__(self)
        self._variables = [a,b]
        self.set_parameter(dist)
    
    def satisfied(self, mapping):
        """return True iff mapping from variable names to points satisfies constraint""" 
        a = mapping[self._variables[0]]
        b = mapping[self._variables[1]]
        result = tol_eq(distance_2p(a,b), abs(self._value))
        return result

    def __str__(self):
        return "DistanceConstraint("\
            +str(self._variables[0])+","\
            +str(self._variables[1])+"="\
            +str(self._value)+")" class AngleConstraint(ParametricConstraint):
    """A constraint on the angle in point B of a triangle ABC"""
    
    def __init__(self, a, b, c, ang):
        """Create a new AngleConstraint instance. 
        
           keyword args:
            a    - a point variable name
            b    - a point variable name
            c    - a point variable name
            ang  - the angle parameter value
        """
        ParametricConstraint.__init__(self)
        self._variables = [a,b,c]
        self.set_parameter(ang)

    def satisfied(self, mapping):
        """return True iff mapping from variable names to points satisfies constraint""" 
        a = mapping[self._variables[0]]
        b = mapping[self._variables[1]]
        c = mapping[self._variables[2]]
        ang = angle_3p(a,b,c)
        if ang == None:
            # if the angle is indeterminate, its probably ok. 
            result = True
        else:
            # in 3d, ignore the sign of the angle
            if len(a) >= 3:
                cmp = abs(self._value)
            else:
                cmp = self._value
            result = tol_eq(ang, cmp)
        if result == False:
            diag_print("measured angle = "+str(ang)+", parameter value = "+str(cmp), "satisfied")
        return result

    def __str__(self):
        return "AngleConstraint("\
            +str(self._variables[0])+","\
            +str(self._variables[1])+","\
            +str(self._variables[2])+"="\
            +str(self._value)+")" class FixConstraint(ParametricConstraint):
    """A constraint to fix a point relative to the coordinate system"""

    def __init__(self, var, pos):
        """Create a new FixConstraint instance
        
           keyword args:
            var    - a point variable name 
            pos    - the position parameter
        """
        ParametricConstraint.__init__(self)
        self._variables = [var]
        self.set_parameter(vector.vector(pos))

    def satisfied(self, mapping):
        """return True iff mapping from variable names to points satisfies constraint""" 
        point = mapping[self._variables[0]]
        if len(point) != len(self._value):
            diag_print("warning: FixConstraint.satisfied: vectors of unequal length", "geometric.FixConstraint.satisfied")
            return False
        result = True;
        for i in range(len(self._value)):
            result &= tol_eq(point[i], self._value[i])
        return result

    def __str__(self):
        return "FixConstraint("\
            +str(self._variables[0])+"="\
            +str(self._value)+")" class vector(list):
        """
        A list based vector class
        """
        # no c'tor

        def __getslice__(self, i, j):
                try:
                        # use the list __getslice__ method and convert
                        # result to vector
                        return vector(super(vector, self).__getslice__(i,j))
                except:
                        raise TypeError, 'vector::FAILURE in __getslice__'
                
        def __add__(self, other):
                return vector(map(lambda x,y: x+y, self, other))

        def __neg__(self):
                return vector(map(lambda x: -x, self))
        
        def __sub__(self, other):
                return vector(map(lambda x,y: x-y, self, other))

        def __mul__(self, other):
            """
            Element by element multiplication
            """
            try:
                    return vector(map(lambda x,y: x*y, self,other))
            except:
                    # other is a const
                    return vector(map(lambda x: x*other, self))


        def __rmul__(self, other):
                return (self*other)


        def __div__(self, other):
            """
            Element by element division.
            """
            try:
                    return vector(map(lambda x,y: x/y, self, other))
            except:
                    return vector(map(lambda x: x/other, self))

        def __rdiv__(self, other):
            """
            The same as __div__
            """
            try:
                    return vector(map(lambda x,y: x/y, other, self))
            except:
                    # other is a const
                    return vector(map(lambda x: other/x, self))

        def size(self): return len(self)

        def conjugate(self):
            return vector(map(lambda x: x.conjugate(), self))

        def ReIm(self):
                """
                Return the real and imaginary parts
                """
                return [
                        vector(map(lambda x: x.real, self)),
                        vector(map(lambda x: x.imag, self)),
                        ]
        
        def AbsArg(self):
                """
                Return modulus and phase parts
                """
                return [
                        vector(map(lambda x: abs(x), self)),
                        vector(map(lambda x: math.atan2(x.imag,x.real), self)),
                        ]


        def out(self):
            """
            Prints out the vector.
            """
            print self

        def concatonated(self,other):
            """this vector concatonated with another"""
            return vector(list(self)+list(other)) def random_problem_2D(numpoints, radius=10.0, roundoff=0.0, angleratio=0.5):
    """Generate a random problem with given number of points, a roundoff
       value for the prototype points, a radius for the cloud of prototype points
       and a ratio of angle constraints over distance constraints"""
    group = {}
    problem = GeometricProblem(dimension=2)
    i = 0
    while i < numpoints:
        aname = 'p'+str(i)
        apoint = vector([
            _round(random.uniform(-radius,radius),roundoff),
            _round(random.uniform(-radius,radius),roundoff)
        ])
        unique = True
        for v in group:
            p = group[v]
            if tol_eq(apoint[0],p[0]) and tol_eq(apoint[1],p[1]):
                unique = False
                break
        if unique:
                problem.add_point(aname, apoint)
                group[aname] = apoint
                i = i + 1
    #next
    _constraint_group(problem, group, None, angleratio)
    return problem def diag_select(pattern):
    """Set regexp pattern to filter which messages are printed.""" 
    global diag_selector 
    diag_selector = re.compile(pattern) def diag_print(str, code=''):
    global diag_selector
    global diag_stream
    if diag_selector.match(code):
        diag_stream.write(str)
        diag_stream.write("\n") def distance_2p(p1, p2):
    """Returns the euclidean distance between two points
       arg keywords:
          p1 - a vector
          p2 - a vector
       returns: a number
    """
    return vector.norm(p2 - p1)
 def angle_3p(p1, p2, p3):
    """Returns the angle, in radians, rotating vector p2p1 to vector p2p3.
       arg keywords:
          p1 - a vector
          p2 - a vector
          p3 - a vector
       returns: a number
       In 2D, the angle is a signed angle, range [-pi,pi], corresponding
       to a clockwise rotation. If p1-p2-p3 is clockwise, then angle > 0.
       In 3D, the angle is unsigned, range [0,pi]
    """
    d21 = vector.norm(p2-p1)
    d23 = vector.norm(p3-p2)
    if tol_eq(d21,0) or tol_eq(d23,0):
        # degenerate, indeterminate angle
        return None
    v21 = (p1-p2) / d21
    v23 = (p3-p2) / d23
    t = vector.dot(v21,v23) # / (d21 * d23)
    if t > 1.0:             # check for floating point error
        t = 1.0
    elif t < -1.0:
        t = -1.0
    angle = math.acos(t)
    if len(p1) == 2:        # 2D case
        if is_counterclockwise(p1,p2,p3):
            angle = -angle
    return angle

import random
from test_generic import test
from geosolver.geometric import GeometricProblem, GeometricSolver, DistanceConstraint, AngleConstraint, FixConstraint
from geosolver.vector import vector 
from geosolver.randomproblem import random_problem_2D
from geosolver.diagnostic import diag_select, diag_print
from geosolver.intersections import distance_2p, angle_3p
#!/usr/bin/env python
"""This module provides some tests for the GeoSolver. 
These test are for 2D solving.
The tests are also simple examples of how to use of the GeomSolver API"""


# -------- 2D problems

def ddd_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))
    problem.add_constraint(DistanceConstraint('v1','v2',distance_2p(problem.get_point('v1'), problem.get_point('v2'))))
    problem.add_constraint(DistanceConstraint('v1','v3',distance_2p(problem.get_point('v1'), problem.get_point('v3'))))
    problem.add_constraint(DistanceConstraint('v2','v3',distance_2p(problem.get_point('v2'), problem.get_point('v3')))) 
    return problem

def dad_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))
    problem.add_constraint(DistanceConstraint('v1','v2',
        distance_2p(problem.get_point('v1'), problem.get_point('v2'))
    ))
    problem.add_constraint(DistanceConstraint('v2','v3',
        distance_2p(problem.get_point('v2'), problem.get_point('v3'))
    )) 
    problem.add_constraint(AngleConstraint('v1','v2','v3', 
       angle_3p(problem.get_point('v1'), problem.get_point('v2'), problem.get_point('v3'))
    ))
    return problem

def add_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))
    problem.add_constraint(DistanceConstraint('v1','v2',
        distance_2p(problem.get_point('v1'), problem.get_point('v2'))
    ))
    problem.add_constraint(DistanceConstraint('v2','v3',
        distance_2p(problem.get_point('v2'), problem.get_point('v3'))
    )) 
    problem.add_constraint(AngleConstraint('v3','v1','v2', 
       angle_3p(problem.get_point('v3'), problem.get_point('v1'), problem.get_point('v2'))
    ))
    return problem


def aad_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))
    problem.add_constraint(DistanceConstraint('v1','v2',distance_2p(problem.get_point('v1'), problem.get_point('v2'))))
    problem.add_constraint(AngleConstraint('v2', 'v3', 'v1', 
       angle_3p(problem.get_point('v2'), problem.get_point('v3'), problem.get_point('v1'))
    ))
    problem.add_constraint(AngleConstraint('v3', 'v1', 'v2',
       angle_3p(problem.get_point('v3'), problem.get_point('v1'), problem.get_point('v2'))
    ))
    return problem

def ada_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))
    problem.add_constraint(DistanceConstraint('v1','v2',distance_2p(problem.get_point('v1'), problem.get_point('v2'))))
    problem.add_constraint(AngleConstraint('v3', 'v1', 'v2', 
       angle_3p(problem.get_point('v3'), problem.get_point('v1'), problem.get_point('v2'))
    ))
    problem.add_constraint(AngleConstraint('v1', 'v2', 'v3',
       angle_3p(problem.get_point('v1'), problem.get_point('v2'), problem.get_point('v3'))
    ))
    return problem



def propagation_problem():
    problem = GeometricProblem(dimension=2)
    problem.add_point('v1', vector([random.random() for i in [1,2]]))
    problem.add_point('v2', vector([random.random() for i in [1,2]]))
    problem.add_point('v3', vector([random.random() for i in [1,2]]))

    problem.add_point('v4', vector([random.random() for i in [1,2]]))