Bits of Learning

Learning sometimes happens in big jumps, but mostly in little tiny steps. I share my baby steps of learning here, mostly on topics around programming, programming languages, software engineering, and computing in general. But occasionally, even on other disciplines of engineering or even science. I mostly learn through examples and doing. And this place is a logbook of my experiences in learning something. You may find several things interesting here: little cute snippets of (hopefully useful) code, a bit of backing theory, and a lot of gyan on how learning can be so much fun.

Thursday, May 04, 2017

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).

Wednesday, September 21, 2016

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()

Tuesday, September 13, 2016

APIs with Functions

This post introduces the concept of application programming interfaces (APIs) and what they are good for. As an example, we continue to use the Hisaab program.

As we have seen, functions are a very powerful abstraction mechanism to wrap around pieces of computation which can then be used and re-used as needed. This leads to a code which is more compact, readable, modifiable and flexible.

In Hisaab, the implementation we first created uses a particular representation of the graph data-structure called the edge list. For example, the following graph will be represented by the code right below that:



[
  ("A", 15, "F"),
  ("A", 20, "B"),
  ("B", 10, "D"),
  ("B", 10, "C"),
  ("D", 30, "C"),
  ("C", 30, "G"),
  ("G", 25, "E"),
  ("F", 100, "E"),
  ("E", 150, "D"),
  ("G", 90, "D")
]
The code of the algorithm is reproduced below:



 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
def reroute((n1, w, n2), n):
  if(n1 == n): return [(n2, -w, n1)]
  elif(n1 != n) and (n2 != n): return [(n1, w, n), (n2, -w, n)]
  else: return [(n1, w, n2)]

def starify_graph(g, n):
  new_g = []
  for e in g:
    new_edge = reroute(e, n)
    new_g.extend(new_edge)
  return new_g

# assumption: g is not a multigraph.
def merge((n1, w, n2), g):
  for i in range(len(g)):
    (m1, w1, m2) = g[i]
    if(m1 == n1 and m2 == n2):
      g[i] = (n1, w + w1, n2)
      return g
  g.append((n1, w, n2))
  return g

def graph_of_mgraph(mg):
  g = []
  for e in mg:
    g = merge(e, g)
  return g


The above code depends on the way we have implemented the graph data-structure. There are several examples in the above code which demonstrate the fact that the above code has been written with the knowledge about the internal implementation of the graph data-structure.

  • Line 1, 14: assumes that edges are triples.
  • Line 8, 10, 16, 18, 20, 24: graph is assumed to be an edge list.

The disadvantage with the above implementation of the algorithm is that it is strictly dependent on the way the graph data-structure is implemented. If the implementation of graph changes, the algorithm would have to be pretty much completely re-implemented.

Let's see how the above graph represented as an adjacency list would look like:

{
  'A': [('F', 15), ('B', 20)],
  'C': [('G', 30)],
  'B': [('D', 10), ('C', 10)],
  'E': [('D', 150)],
  'D': [('C', 30)],
  'G': [('E', 25), ('D', 90)],
  'F': [('E', 100)]
}

Below, we should the implementation of the algorithm when the graph is implemented as an adjacency list.


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
def starify_graph(g, n):
  new_g = {} # instantiating a dictionary
  for n1 in g: #initialising the empty graph
    new_g[n1] = []
  for n1 in g:
    # invert and add every node starting with n.
    if(n1 == n):
      for (n2, w) in g[n1]:
        add_edge(new_g[n2], (n1, -w))
    else:
      for (n2, w) in g[n1]:
        # add every node ending at n.
        if(n2 == n):
          add_edge(new_g[n1], (n2, w))
        # reroute all other nodes.
        else:
          add_edge(new_g[n1], (n, w))
          add_edge(new_g[n2], (n, -w))
  return new_g

As is evident, this implementation is significantly different from the previous implementation. Again, there are several evidences in this code to show that it strictly assumes an adjacency list implementation of the graph data-structure.
Line 2: Graph is a dictionary.
Line 4: The elements of the graph dictionary are lists.
Line 8, 9, 11: The elements of the graph dictionary are pairs.

Modularity

Can we implement the above hisaab algorithm in a way so that it is independent of the internal details of the graph data-structure?

If we could do so, there would be the following advantages:
  • The hisaab algorithm would work fine irrespective of which implementation of graph we use.
  • The graph data-structure implementation can be modified subsequently without a risk of breaking the hisaab code.

The answer to the above question is Yes. By defining a set of functions which define a uniform way by which the rest of the code interacts with the graph code, we can ensure that there is separation between the view of the graph that the rest of the code sees, and its internal implementation. We do so, using a group of functions as follows:
  1. empty_graph: Returns a new empty graph with no nodes or edges
  2. make_graph: Returns a new graph having the edges which are provided in a list
  3. get_edges: Returns the list of edges in the graph
  4. add_edge: Adds an edge to a graph
The implementation of the above functions is specific to the specific implementation of the graph, but how they are used is not. Hence, the code which uses a graph only through the above functions is also independent of the specific implementation of the graph data-structure. The above set of functions can be termed as an application programming interface (API). This is because, this set defines the way one part of the program (called the client) interacts with another part of the program (called the server). Further functions like empty_graph and make_graph are called constructors (because they create new instances of the data-structure), get_edges are called extractors (because they reveal information about the data-structure without modifying it) and add_edge are called modifiers (because they modify the underlying data-structure). One way to look at the entire scenario would be as follows:

The client code (or user code) makes use of the facilities provided by the server code (or back-end code). In hisaab, the algorithm (implemented through the functions reroute and starify_graph) can be considered as the client code, while the graph data-structure can be viewed as the server code. As shown above, the API consisting of the four functions listed above form the interface through which the client talks to the server. Note that it is possible to freely change which of the two implementations the interface uses as the server without affecting the client. This gives us the flexibility of independently implementing the client without concerning ourselves about the internal details of the server.



 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
def reroute(e, n):
  n1 = get_start_node(e)
  w  = get_weight(e)
  n2 = get_end_node(e)
  if(get_start_node(e) == n):
    return [make_edge(n1 = n2, w = -w, n2 = n1)]
  elif(n1 != n) and (n2 != n):
    return [make_edge(n1 = n1, w = w, n2 = n), make_edge(n1 = n2, w = -w, n2 = n)]
  else:
    return [e]

def starify_graph(g, n):
  new_g = empty_graph()
  edges = get_edges(g)
  for e in edges:
    new_edges = reroute(e, n)
    for e in new_edges:
      add_edge(g = new_g, e = e)
  return new_g

The above code shows the hisaab algorithm re-implemented using the above API. The best way to get convinced this code is independent of any of the internal implementation details of graph data-structure is to look for a dependence. You will not find one!

The API allows us to implement hisaab as a re-usable algorithm independent of the implementation of the graph.

Large Scale Software Development and Software Design

In large scale software development scenario, software systems are developed by teams of multiple software developers. It would be too inefficient to develop each part of the system sequentially. It's important to ensure that all teams are kept busy throughout the project duration. All sub-teams in a project work in parallel to ensure an early completion of the project.

The way the above is achieved is by dividing the work into several modules which interact with each other to create the complete system. This process is called software design. The division is worked out in a way that there is comparable distribution of work between the various modules. More importantly, each module defines an interface (or API) through which the rest of the system interacts with that module. A good design will ensure that the interfaces are enough to allow powerful interaction between modules but doesn't expose unnecessary details of a module to the rest of the system. While the internal implementation of a module is likely to change rapidly over the execution of a project, the module interfaces must be decided early and frozen. This is important to ensure that all other parts of the system can proceed assuming this interface about the given module. Hence, every sub-team, working on an individual module, can proceed in parallel.

Prototypes and Improvisations

As mentioned above, API enable independent development of interacting modules. Suppose, one team is developing the client C, and another is developing the server S. A typical way this piece of software would be developed in practice is as follows: An API for S is defined first. As C is dependent on S, a quick-and-dirty implement of S -- say S1 -- is created and handed over to the client team so that C's development can proceed. While C's development proceeds, the server team continues working on a better-engineered version of S -- say S2 -- in parallel. S2 would follow the same API of S, however may score better than S1 in certain parameters like performance, security etc. By the time C's development gets over, S2 is also ready. The version of the product which gets shipped to the customer uses C as client S2 as the server.

Code:

hisaab1.py: Implements hisaab using a edge-list implementation
hisaab2.py: Implements hisaab using an adjacency-list implementation
hisaab3.py: Implements hisaab as a re-usable code by using a graph API

Friday, September 02, 2016

Hisaab -- A Program to Simplify Transactions among Friends

It's always fun when your technical knowledge comes in handy while to trying to solve some day-to-day problem. you get a kick in being able to describe a problem in its core mathematical essence, and in creating an algorithm to solve the problem. The size of the kick is independent of the size and complexity of the problem.

Follows a cute example of the above.

Problem Description

Often, a group of people may transact money between each other. A situation may arise when each person may have taken some money from someone else, thus owing the other person that much amount. To settle the accounts, there might be needed a very large number of transactions needed, requiring everyone to meet everyone else.

Solution Approach

good way to simplify the settlement process would be to make one of the people involved as the bank, or the star. Everyone makes a transaction with this person only, only the net amount she owes to anyone.

Mathematical Model






The original situation can be represented as a directed graph as shown above in the left side graph. Each node represents a person, and each directed edge indicates a transaction. For example, the edge with weight 15 between A and F means that A owes Rs. 15 to F.

The solution is shown in the right side graph above. Note that all edges now involve node D. Thus every node needs to make at most one transaction with node D, which is huge simplification of the process. The accounts will now be settled with only 5 transactions as opposed to 9 in the original graph. As an added benefit, nodes which need not transact at all also get identified. For example, node B gets isolated in the right side graph, meaning that it doesn't need to transact anything with node D.


Transforms

The solution is achieving through the following two simple graph transformations.

Rerouting



Suppose we wish to make node 1 the star node. This means that we want all transactions between any pair not involving node 1 should be routed through node 1, as shown above.

Multigraph to Graph




Implementation

Follows an implementation of the above idea in OCaml. We use a purely functional style (i.e. by using a completely side-effect free style of programming). If you remove the profuse amount in-line documentation, you will notice that the code is remarkably concise (just about 25 lines of code)!

Find an OCaml and Python implementation here.

  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
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
(*
The following program is useful to simplify the transactions amongst a group of people.
Problem Description: Often, a group of people may transact money between each other. A situation may arise when each person
may have taken some money from someone else, thus owing the other person that much amount. A
good way to simplify the settlement process would be to make one of the people involved as the
bank, or the star. Everyone makes a transaction with this person only, only the net amount she
owes to anyone.

User Instruction:
1. Open the file in a editor.
2. Edit the variable g to reflect what each person owes any other.
3. Save the file
4. Open the OCaml REPL
5. enter the command: #use "hisaab.ml";;

The answer will appear on the standard output.
*)

(*
  Description: Takes an edge e = (n1, w, n2) and a node n, and returns a list of edges, to
    apply starification transform on e.

  Signature: 'a * int * 'a -> 'a -> ('a * int * 'a) list

  Example:
  1)
# reroute (1, 10, 2) 1;;
- : (int * int * int) list = [(2, -10, 1)]
# reroute (2, 10, 1) 1;;
- : (int * int * int) list = [(2, 10, 1)]
# reroute (2, 10, 3) 1;;
- : (int * int * int) list = [(2, 10, 1); (3, -10, 1)]
*)
let reroute (n1, w, n2) n =
  if      n1 = n then [(n2, -w, n1)]
  else if n2 = n then [(n1, w, n2)]
  else                [(n1, w, n); (n2, -w, n)]

(*
  Description: Takes a list of lists and returns the corresponding list

  Signature: 'a list list -> 'a list

  Example:
# flatten_list [[1;2]; [3;4;5];[6]];;
- : int list = [1; 2; 3; 4; 5; 6]
*)
let rec flatten_list = function
    [] -> []
  | h :: t -> List.append h (flatten_list t)
 
(*
  Description: Takes a graph (could be multigraph) g and and a node n and returns
  and equivalent graph g' which is starred at n, and preserves the net input/output
  at each node of g.

  Signature: ('a * int * 'a) list -> 'a -> ('a * int * 'a) list
  Example:
# starify_graph [(1, 160, 2); (2, 440, 1); (1, 300, 3); (3, 160, 2)] 1;;
- : (int * int * int) list =
[(3, 160, 1); (2, -160, 1); (3, -300, 1); (2, 440, 1); (2, -160, 1)]
*)
let starify_graph g n =
  let list_of_lists = List.map (fun g -> reroute g n) g in
    flatten_list list_of_lists

(*
  Description: Takes an edge e = (n1, w, n2) and a graph g, and returns a 
  new graph g' such that g' is same as g except that if g contains an edge
  e' with same source and target nodes as e, then e and e' are merged by
  adding their weights.

  Signature: 'a * int * 'b -> ('a * int * 'b) list -> ('a * int * 'b) list

  Assumption: g is not a multi-graph.

  Example:
# let g2 = [(3, -140, 1); (2, 120, 1)];;
val g2 : (int * int * int) list = [(3, -140, 1); (2, 120, 1)]
# merge (3, 140, 1) g2;;
- : (int * int * int) list = [(3, 0, 1); (2, 120, 1)]
# merge (3, 140, 2) g2;; 
- : (int * int * int) list = [(3, -140, 1); (2, 120, 1); (3, 140, 2)]
# merge (3, 140, 2) [];;
- : (int * int * int) list = [(3, 140, 2)]
*)
let rec merge (n1, w, n2) = function
    [] -> [(n1, w, n2)]
  | (n1', w', n2') :: t ->
      if n1 = n1' && n2 = n2' then ((n1, w + w', n2) :: t)
      else (n1', w', n2') :: (merge (n1, w, n2) t)

(*
  Description: Takes a multigraph mg and returns its corresponding graph g by merging
    all edges between common pairs of nodes.

  Signature: ('a * int * 'b) list -> ('a * int * 'b) list

  Example:
# graph_of_mgraph [(1, 160, 2); (2, 300, 1); (2, 140, 1); (1, 300, 3); (3, 160, 2)];;
- : (int * int * int) list =
[(1, 160, 2); (2, 440, 1); (1, 300, 3); (3, 160, 2)]
*)  
let graph_of_mgraph mg =
  let rec iter g = function
    [] -> g
  | e :: t -> iter (merge e g) t
  in
  iter [] mg

(*
  Description: Turns the weight of every edge of a given graph to positive value by inverting the direction of the edge if its weight is negative.

  Signature: ('a * int * 'a) list -> ('a * int * 'a) list
  
  Example:
# abs_weight [("Taru", -140, "Shraddha"); ("Shilpi", 20, "Shraddha")];;
- : (string * int * string) list =
[("Shraddha", 140, "Taru"); ("Shilpi", 20, "Shraddha")]    
*)
let abs_weight g =
  List.map (
    fun (n1, w, n2) -> if w >= 0 then (n1, w, n2) else (n2, -w, n1)
  ) g

let remove_empty_edges g =
  List.filter (fun (_, w, _) -> w <> 0) g
(*
  Example graph (Provided by Shilpi)
  node 1 -> Shilpi
  node 2 -> Shraddha
  node 3 -> Taru

  (1, 160, 2) -> Shilpi owes Rs. 160 to Shraddha
*)
let g = [("Shilpi", 160, "Shraddha"); ("Shraddha", 300, "Shilpi"); ("Shraddha", 140, "Shilpi"); ("Shilpi", 300, "Taru"); ("Taru", 160, "Shraddha")]
let star_node = "Shraddha"
let _ = remove_empty_edges (
          abs_weight (
            graph_of_mgraph (
              starify_graph g star_node)))

let e1 = [
  ("r2", 10, "r3");
  ("r1", 20, "r2");
  ("r2", 10, "r4");
  ("r4", 30, "r3");
  ("r1", 15, "r6");
  ("r6", 100,"r5");
  ("r7", 25, "r5");
  ("r5", 150,"r4");
  ("r3", 30, "r7");
  ("r7", 90, "r4");
]

let star_node = "r4"
let _ = remove_empty_edges (
          abs_weight (
            graph_of_mgraph (
              starify_graph e1 star_node)))