Higher Order Functions and Tail Recursion

Calculating PI using Monte Carlo Simulation

Monte Carlo Simulation to Compute PI

Equations

Code

Follows a sample OCaml code to implement the above idea.

let one_point () = (Random.float (2.) -. 1., Random.float (2.) -. 1.)

let n_points n =
  let rec iter acc = function
      0 -> acc
    | n -> iter ((one_point ()) :: acc) (n - 1)
  in
  iter [] n

(* tail-recursive implementation of map to prevent stack-overflow. *)
let mymap f l =
  let rec iter f acc = function
      [] -> acc
    | h :: t -> iter f ((f h) :: acc) t
  in
  iter f [] l

(* tail-recursive implementation of filter to prevent stack-overflow. *)
let myfilter f l =
  let rec iter f acc = function
      [] -> acc
    | h :: t -> if (f h) then (iter f ((f h) :: acc) t)
                else (iter f acc t)
  in
  iter f [] l

let pi =
  let total = 10000000 in                                         (* number of times the dart is thrown *)
  let points = n_points total in                                  (* list of points all the darts hit *)
  let rs = mymap (fun (x, y) -> x *. x +. y *. y) points in       (* the distance from origin of all points *)
  let inner_points = myfilter (fun r -> r <= 1.) rs in            (* list of all points inside or on the unit circle *)
  let num_inner_points = List.length inner_points in              (* number of points inside or on the unit circle *)
  4. *. ((float_of_int num_inner_points) /. (float_of_int total)) (* value of PI. *)

OCaml Interpretor Interaction

# #use "pi.ml";;
val one_point : unit -> float * float = <fun>
val n_points : int -> (float * float) list = <fun>
val mymap : ('a -> 'b) -> 'a list -> 'b list = <fun>
val myfilter : ('a -> bool) -> 'a list -> bool list = <fun>
val pi : float = 3.1419352

Notes

A large number of points are generated within the bounding square. The probability distribution for the points over this area is uniform, which means that the probability of occurrence of a dot at all points in the given square is equal. With large values of N, the ratio of the number of points inside the unit circle to N starts approaching the ratio of the areas of the unit circle and the square. This idea is used to compute a close approximate to Pi.

As we use OCaml's library module Random to generate the points, the value of Pi computed in every simulation is non-deterministic and is different from the other by a small amount. The time needed to do the above computation was about 10 seconds on my laptop.

Most of the lines in the above code are really there to either enhance the readability or scalability of the program. As List.map and List.filter are non-tail recursive implementations, they lead to stack overflow for very large values of total (>= 1000000). The tail-recursive implementations of map and filter are provided only to prevent stack overflow for very large values of total (>= 1000000).

Higher Order Functions in Python

Pardon the hand-written notes. Will try and type out the whole thing in some future point.

Computing Unions and Differences for Sets of Time Durations

Timeline

This article presents an algorithm to compute set union and set minus in a specific context. Each set is a set of contiguous stretches on a real line.

Problem Description

Union

Consider a scenario where two people are teamed up to contribute to a particular task, say attending to an infant child. As long as any one of them is available to look after the child, it's fine. For other times, they probably would have to look for an alternative arrangement. The time any one person is available can be represented as a set of durations. The time for which at least one of the two persons is available is nothing but the union of the sets of durations during which each one of them is available.

Set Minus

Consider a scenario when two people are looking for a slot wherein both are available for a discussion. This can be thought of a difference between the the available time of one person and the busy time of the other.

Solution Approach

The solution approach makes use of a variable level which keeps track of the starts and ends of durations. After placing all the important time instants on a time line, we scan the time line from left to right.

While computing union, we use the following strategy to set the level during our scan of the timeline:

1. Each time any of the duration begins, we increment level by 1.
2. Each time any of the durations ends, we decrement level by 1.

As per the above, the stretches of the timeline with level = 0 indicate stretches when none of the two sets is present.

The rest of the timeline will have positive values of level indicating that at least one of the two sets is present there.

In another pass of the timeline, we now collect all the continuous stretches of the timeline where level is not zero. This gives us the set of durations corresponding to the union of the two sets.

While computing set minus, we employ the following strategy for maintaining level.

1. The starting point of a duration in the first set and end of a duration in the second set cause level to be incremented by 1.
2. The end point of a duration in the first set and start of a duration in the second set cause level to be decremented by 1.

As per the above, the stretches of the timeline with level = 1 indicates stretches in the difference between the first set and the second.

In another pass, we collect the set of contiguous portions of the timeline where level = 1. This is our answer of set minus.

Extensions

The union and setminus functions can be used to implement other useful set operations on the above kind of sets, viz. complement and intersection.

Implementation

Representation of the sets

We represent the sets of durations as lists of pairs. s1 and s2 in the above figure (Timeline) have been represented as follows:

  s1 = [(1, 5), (8, 10), (15, 20)]
  s2 = [(2, 4), (6, 12), (16, 21)]

Code Design

When we implement the above two algorithms as two completely separate functions union and setminus, we realise that they have a very common structure. This gives us the motivation to try and capture this commonality between the two functions in a piece of re-usable code. We do so by designing a higher-order function make_function, which is parameterised by another function which captures the specific elements of the two functions. f1 and f2 are the two functions capturing these specificities of union and set minus.

The make_function function prepares a function and returns it as a result. the two functions union and setminus are created by calling make_function with f1 and f2 as arguments respectively. Note that this could have been implemented alternatively by embedding a call to make_function in union and setminus. This is a perfectly good option when union and setminus get used only once. However, considering a scenario where union and setminus may get called several times in a larger program, the strategy used here will lead to a slightly improved speed, not to mention improved readability (arguable).

Code

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67

# utility function for the implementation of union def f1(s1, s2): START = 1 END = -1 l = [] for (s, e) in s1: l.append((s, START)) l.append((e, END)) for (s, e) in s2: l.append((s, START)) l.append((e, END)) return sorted(l, key = lambda (x, y): x) #utility function for implementation of setminus def f2(s1, s2): START1 = 1 END1 = -1 START2 = -1 END2 = 1 l = [] for (s, e) in s1: l.append((s, START1)) l.append((e, END1)) for (s, e) in s2: l.append((s, START2)) l.append((e, END2)) return sorted(l, key = lambda (x, y): x) # Higher order function to create union and setminus functions def make_function(f): # boilerplate function to be returned from make_function after being plugged with f. def g(s1, s2): level = 0 result = [] l = f(s1, s2) for (t, SE) in l: print (t, level) last = level level += SE if level == 1 and last == 0: start = t if level == 0 and last == 1: end = t result.append((start, end)) return result return g union = make_function(f1) setminus = make_function(f2) def t1(): s1 = [(1, 5), (8, 10), (15, 20)] s2 = [(2, 4), (6, 12), (16, 21)] print union(s1, s2) def t2(): s1 = [(1, 5), (8, 10), (15, 20)] s2 = [(2, 4), (6, 12), (16, 21)] print setminus(s1, s2) if __name__ == "__main__": t1() t2()

set.py

Counter in Python

Presenting here two implementations of a counter.

What does the counter do?

Given a base n, this counter will allow you to do base n counting. What's base-n counting? There are n digits: 0, 1, 2, ..., n - 1. At each place, the digit changes from i to i + 1 if i < n - 1, and to 0 if i = n - 1. This change happens when the change at that place is triggered. And when does the trigger happen. For the unit's place, the trigger happens whenever the function is called. For all other places, the trigger happens when the previous place's digit rolls over from n - 1 to 0. That's all. That's a bit technical description for the simple act of counting we learn in our KG. But if we want to implement the counter, it needs to be written in this algorithmic language.

A 3-bit binary (base-2) counter goes like:

000, 001, 010, 011, 100, 101, 110, 111, 000, ...

A 2-digit base-4 counter would count as:

00, 01, 02, 03, 10, 11, 12, 13, 20, 21, 22, 23, 30, 31, 32, 33, 00, ...

What do we usually do in day-to-day counting? decimal (base-10) counting.

How do we Use the Counter Program?

Below is a piece of code that could be used to implement a binary counter on a 3-bit number.

base = 2
inc = makeCounter_iter(base)
n = [0, 0, 0]
print n
for i in range(base ** len(n)):
    n = inc(n)
    print n

 The inc(n) gives n + 1 in the base n number system which is pre-decided while instantiating the counter:

base = 2

#!/usr/bin/python
def makeCounter_iter(base):
    def isZero(lst): return lst == [0] * len(lst)

    def inc(num):
        if(len(num) > 1):
            rem = inc(num[:-1])
            if(isZero(rem)):
                if(num[-1] == base - 1):    rem.extend([0]) 
                else:                rem.extend([num[-1] + 1])
            else:    rem.extend([num[-1]])
            return rem
        else:
            new = [0]
            if(not num[0] == base - 1):    new[0] = num[0] + 1
            return new
    return inc

Now follows the recursive version.

#!/usr/bin/python
def makeCounter_rec(base):
    def incDigit(num, pos):
        new = num[:]
        if(num[pos] == base - 1):
            new[pos] = 0
            if(pos < len(num) - 1):    return incDigit(new, pos + 1)
        else:
            new[pos] = new[pos] + 1
        return new

    def inc(num): return incDigit(num, 0)
    return inc

Any day, the recursive implementation trumps over the iterative version. The recursive implementation is smaller and more readable because it's very like the way the process of counting is formally defined in the beginning of the article.

Functional Programming Flavour

One thing to notice is that we use functional language aspect of python. The function makeCounter (makeCounter_iter and makeCounter_rec) is a higher order function that returns another function inc, the increment function. Also, the value of base is fixed in the creation context of inc inside the makeCounter function. The inc function uses the value of base at a place in the code which is typically outside the lexical scope where base is defined. This is known as the closure feature of functional programming. In a typical object-oriented scenario, this effect is achieved using the builder design pattern.