Posts Tagged: algorithms


6
Jan 14

Cyclical dependency detection in the database

We recently had a need to find cyclical dependencies in our database. This happens to be a rather straightforward graph algorithm issue. The database foreign key constraints form a directed graph. Finding a cycle in a directed graph is mostly detecting an already visited node in a DFS algorithm (back-edge). We mark nodes as visited and if the ancestor of a node in the tree is already visited, then a back-edge (cycle) exists.

In order to do this on our own, we’d have to read the metadata from the database for each table, construct a directed graph using the foreign keys and then run the algorithm discussed above, rather straightforward. Most of the complexity comes from the cross cutting concerns of database metadata munging. We can easily accomplish all of the above using sqlalchemy and its ability to perform a topological sort on the reflected tables. Topological sort fails in there is a cycle detected and the exception thrown includes the nodes that produce the back-edge. Using this simple trick, we allow sqlalchemy to detect the cycles for us.

You’ll need to install sqlalchemy (and your db driver), networkx and graphviz (for visualization).


31
Oct 11

Distributed locking made easy

There are various situations when one would use a simple mutex to ensure mutual exclusion to a shared resource. This is usually very simple to accomplish using your favorite language library, but that constrains you to a single process mutual exclusion. Even single machine mutual exclusion is rather straight forward, usually just locking a resource (i.e. file) and awaiting for the lock to be released. You can use that for an IPC mutex.

But what if one needs a distributed mutex to allow mutual exclusion amongst distributed clients? At that, the mutex has to offer various guarantees, as is with any shared state. We know shared state is hard to reason about and provokes a lot of bug-prone software, shared distributed state, is much harder, requiring distributed guaranteed consensus all while operating in a non-reliable network environment. There are various distributed consensus algorithms, Paxos being one of the more widely used ones.

But you can deploy your own distributed locking service, without having to implement your own. Apache Zookeeper, offers distributed synchronization and group services. Mutex/locking is just one of the things you can do with Zookeeper. Zookeeper is rather low level, so implementing a distributed lock, although trivial, requires some boilerplate.

Recently we needed a distributed lock service, to ensure only one person in our organization is performing a particular systems activity at any given point and time. We implemented it on top of a homegrown tool written in ruby. The example below is in ruby, though the api calls would translate to any language…

The usage of this Lock class is such:

    Lock.new("localhost:2181").with_lock(
      "/myapp",
      5, ## Timeout in seconds
      lambda { ## Timeout callback
        abort("Couldn't acquire lock.  Timeout.") },
      lambda { ## Do what ever you want here  }
    )

The details of the algorithm are outlined here.

Of course, before you use it, you must install Zookeeper and create the root path /myapp in order to be able to use it.

Also, please note, I have removed the access control part from the example. In order to use this in production, I strongly encourage you read this.


10
Feb 11

Merge sort in lisp

I’m currently taking (online) MIT’s algorithms course. I needed a good brush-up on algorithms, starting with the basic theory. I’m going to jump in to help a buddy working on a graph database in lisp Vivace Graph, which should be an open source competitor to AllegroGraph. Allegro has probably the best graph db out there right now. There are some competitors in the java world, but they are either expensive or don’t provide enough features, or don’t provide APIs outside of java. Allegro has an extensive list of client libs. I’ll write up more on the graph db in some future, for now, I got back into lisp and implemented the simplest alg (merge sort) in it.

The reason I’m posting is to get some comments. I’d love to hear your thoughts on how a lisp nOOb can improve this 🙂


21
Jun 10

Lazy cheap flight calculations with priority queues

There is an interesting problem of utilizing priority queues to figure out the best price combination in a set of flight legs. The problem is as follows:

We need to calculate the cheapest combination of flight legs (connections) for a flight to a particular destination. We’re given a price ordered N set of flight legs and we need to find the winning combination. Each combination would be evaluate for eligibility and would either pass or fail, so the cheapest combination doesn’t necessarily reflect the cheapest possibly combination of prices from the legs. A black box predicate function is consulted to ensure the combination is eligible. This reflects various airline rules, like overlapping times, specials that are only available to certain people, routes, or connections.

Solution: A naive approach for say a two leg flight is to say construct a (n x m) ordered matrix and evaluate each priced ordered combination through the black box predicate routing until one passes. The problem with this approach is that we unendingly construct a full matrix when in many cases one of the first combinations is enough to present the cheapest “valid” price. The key to reducing this is to construct a lazy data structure which will prioritize the cheapest flights and can then be iterated to find one that’s valid. We do so at runtime while constructing matrix combinations. The solutions is generalized, so the same can be used for n leg flights.

The algorithm goes something like this…

Construct the first set of combinations which can reflect the cheapest flight. The first cheapest combinations is always n1 + m1. If that doesn’t pass, the next possible set of cheapest combinations is either n2 + m1 or n1 + m2. We then continue to n1 + ma and na + m1, where a is incremented until the end of the route leg set for either leg.

The worst case running time is quadratic O(n2), but because of the lazy data structure, the algorithm runs in rather constant time, depending on how lucky we are that the first few combinations will yield a “rule valid” price combination.

This problem idea came from reading The Algorithm Design Manual by Steven S. Skiena. I recommend this book for anyone wishing to delve into the world of more advanced algorithm design.

Here is the solution in python. You’ve probably noticed I’ve been using a lot of python. Besides the fact that I like the language, python is an incredibly good language for conveying algorithmic ideas in a concise but very readable way.

The only two functions that matter, are cheapest_price and _pick_combo, the rest are just auxiliary functions used to support an OO structure and running a sample.

  import heapq, random, time

  class Route(object):
      """docstring for TicketFinder"""
      def __init__(self):
          self.heap = []
          self.unique = dict()
          self.legs = []
          self.max_leg_len = 0
          self._counter = 0
          self._loop_counter = 0

      def add_leg(self, leg):
          leg.sort()
          self.legs.append(leg)
          leg_len = len(leg)
          if leg_len > self.max_leg_len:
              self.max_leg_len = leg_len

      def cheapest_price(self, pred_func=lambda x: True):
          for i in range(0, self.max_leg_len):
              combo = self._pick_combo(i, pred_func)
              if combo: return combo

      def print_stats(self):
          print("""Legs: %s
  Combos examined: %s
  Loops: %s
  """ % (len(self.legs), self._counter, self._loop_counter))

      def _pick_combo(self, curr_idx, pred_func):
          num_legs = len(self.legs)
          price_combo = [ leg[curr_idx] for leg in self.legs if not curr_idx >= len(leg) ]
          self._add_combo(price_combo)
          cheapest_price = self._eval_price_combo(pred_func)
          if cheapest_price: return cheapest_price
          for idx in range(1, self.max_leg_len-curr_idx):
              for j in range(0, num_legs):
                  if len(self.legs[j]) &lt= (curr_idx+idx): continue
                  combo = []
                  for k in range(0, num_legs):
                      self._loop_counter += 1
                      if j == k:
                          combo.append(self.legs[k][curr_idx+idx])
                      elif curr_idx &lt len(self.legs[k]):
                          combo.append(self.legs[k][curr_idx])
                  self._add_combo(combo)

              cheapest_price = self._eval_price_combo(pred_func)
              if cheapest_price: return cheapest_price

      def _add_combo(self, combo):
          self._counter += 1
          if len(combo) == len(self.legs) and not self.unique.has_key(str(combo)):
              heapq.heappush(self.heap, combo)
              self.unique[str(combo)] = True

      def _eval_price_combo(self, pred_func):
          for i in range(0, len(self.heap)):
              least_combo = heapq.heappop(self.heap)
              if pred_func(least_combo):
                  print("Winning combo: %s" % [ "%.2f" % l for l in least_combo ])
                  return sum(least_combo)
          return None


  ############### Samples below ##################

  def sample_run(num_legs, pred_func):
      print(("#" * 30) + " Sample Run " + ("#" * 30))
      route = Route()
      for i in range(0, num_legs):
          route.add_leg( [ random.uniform(100, 500) for i in range(0, 100) ] )

      start = time.clock()
      price = route.cheapest_price(pred_func)
      calc_time = time.clock() - start

      if price:
          print("Cheapest price: %.2f" % price)
      else:
          print("No valid route found")
      route.print_stats()
      print(("#" * 72) + "\n")

  if __name__ == '__main__':
      sample_run(2, lambda x: True)
      def pred(x):
          for price in x:
              if price &lt 150: return False
          return True
      sample_run(3, pred)

I haven’t thoroughly tested this for correctness besides numerous runs and some basic validation so let me know if you see anything apparently wrong here.

Running the above yields

    ############################## Sample Run ##############################
    Winning combo: ['103.62', '106.40']
    Cheapest price: 210.03
    Legs: 2
    Combos examined: 1
    Loops: 0

    ########################################################################

    ############################## Sample Run ##############################
    Winning combo: ['150.74', '150.25', '173.95']
    Cheapest price: 474.95
    Legs: 3
    Combos examined: 2852
    Loops: 8523

    ########################################################################

For the first sample run, we use a predicate function which yields True, so we never examine anything other than the first combo n1 + m1. For the second sample, I add a predicate function which only accepts any price combination where all legs are above $150. (Of course this is not anything resembling airline rules, just good enough to simulate some sample cases, where the first n combinations are rejected). In the second sample run, we utilized 3 legs and examined 2852 combinations before coming up with the winning leg combination for the route. Each price within the combination is the smallest possible price above $150 for each leg.


27
May 10

Random points in polygon generation algorithm

I needed to generate a set of random points within a polygon, including convex and concave. The need arouse in a geospatial domain where polygons are rather small (on a geo-scale) and wouldn’t span more than say 10 miles, though the benefit of employing more complex algorithms to deal with spheroid properties are negligible. Plane geometry provided enough to meet this requirement. Point-in-Polygon tests are rather simple and are used to test whether a point exists in a polygon. The test is performed using a Ray casting algorithm which test the intersections of a ray across the x-axis starting from the point in question.

Another concept is the Minimum Bounding Rectangle (Bounding Box), which is the minimal rectangle needed to enclose a geographical object (i.e. polygon).

So, one can generate random points within a polygon by…

  1. Generating a bounding box
  2. Generating a point within the bounding box. This is a simple algorithm.
  3. Using Point-in-Polygon to test whether this point exists within the polygon.

Because of the random sampling nature and false positives from step 2, which must be tested in step 3, the above must be performed in a loop until the Point-in-Polygon test passes.

This works quite well for generating test data, as there are no tight bounds on the performance characteristics of random generation. One could also use the above algorithm in production as long as the ration of polygon to bounding box is rather large, which is usually the case for convex polygons. The ratio might be too small convex polygons, though causing a more than acceptable number of false positives in step #2.

I’ve implemented this in the geo-utils python package and made available on github. Feel free to use and provide any feedback.

To utilize the geo-utils to generate random points within a polygon, you would do the following:

  from vtown import geo
  from vtown.geo.polygon import Polygon


  polygon = Polygon(  geo.LatLon(42.39321,-82.92114),
                      geo.LatLon(42.39194,-82.91669),
                      geo.LatLon(42.39147,-82.91796),
                      geo.LatLon(42.39090,-82.91974),
                      geo.LatLon(42.39321,-82.92114))

  point = polygon.random_point()

The above polygon is generated using lat/lon coordinates, but you can generate them using simple x/y coordinates with geo.Point(x,y)

Here are some code snippets from the implementation. I only pasted the relevant parts. For boilerplate and relevant data structures, see the geo-utils package.

class BoundingBox(object):

    def __init__(self, *points):
        """docstring for __init__"""
        xmin = ymin = float('inf')
        xmax = ymax = float('-inf')
        for p in points:
            if p.x < xmin: xmin = p.x
            if p.y < ymin: ymin = p.y
            if p.x > xmax: xmax = p.x
            if p.y > ymax: ymax = p.y
        self.interval_x = Interval(xmin, xmax)
        self.interval_y = Interval(ymin, ymax)

    def random_point(self):
        x = self.interval_x.random_point()
        y = self.interval_y.random_point()
        return Point(x, y)

class Polygon:
  ## __init__ omitted here...

  def contains(self, point):
        seg_counter = private.SegmentCounter(point)
        for i in range(1, len(self.points)):
            line = Line(*self.points[i-1:i+1])
            if seg_counter.process_segment(line):
                return True
        return seg_counter.crossings % 2 == 1

  def random_point(self):
        bb = BoundingBox(*self.points)
        while True:
            print("GENERATING RANDOM POINT...")
            p = bb.random_point()
            if self.contains(p):
                return p

class SegmentCounter(object):

    def __init__(self, point):
        self.point = point
        self.crossings = 0

    def process_segment(self, line):
        p, p1, p2 = self.point, line.point1, line.point2
        if p1.x < p.x and p2.x < p.x:
            return False

        if (p.x == p2.x and p.y == p2.y):
            return True

        if p1.y == p.y and p2.y == p.y:
            minx = p1.x
            maxx = p2.x
            if minx > maxx:
                minx = p2.x
                maxx = p1.x
            if p.x >= minx and p.x <= maxx:
                return True
            return False


        if ((p1.y > p.y) and (p2.y <= p.y)) \
                or ((p2.y > p.y) and (p1.y <= p.y)):
            x1 = p1.x - p.x
            y1 = p1.y - p.y
            x2 = p2.x - p.x
            y2 = p2.y - p.y

            det = numpy.linalg.det([[x1, y1], [x2, y2]])
            if det == 0.0:
                return True
            if y2 < y1:
                det = -det

            if det > 0.0:
                self.crossings += 1