Showing posts with label Algorithms. Show all posts
Showing posts with label Algorithms. Show all posts

Sunday, August 11, 2019

Partitions of a set

Calculating the partitions of a set

Having "solved" a bunch of these divide & conquer problems, I'm the first to admit to having being lulled into a false sense of security. At first glance, the problem of this post seemed deceptively simple and consequently I struggled with it, sort of "hand-waving", not really engaging my brain and getting more and more frustrated how the dang thing wouldn't yield to my experience! I think the moral of the story is math doesn't care about your previous successes and so don't let your past practice trick you into laziness. Be guided by your experience but fully apply yourself to the problem at hand!

Suppose a set of two elements {2, 3}. There are only two ways it can be partitioned: (23), (3)(2). For meaning, you might think of these two partitions like this : in the first partition, there is a connection between the elements 2 and 3, in the second, 2 and 3 are isolated from each other.

Suppose a set of elements {1, 2, 3}. There are five partitions of this set : (123), (23)(1), (13)(2), (3)(21), (3)(2)(1) (I've carefully written them out this way to help with the elucidation). Maybe you want to break here and see if you can write an algorithm for calculating them before reading on?

Observe that we can get the partitions of {1, 2, 3} from knowledge of the partitions of {2, 3} by looking at each partition of {2, 3} in turn and considering the partitions that would result by inclusion of the element 1. So, for example, the partition (23) gives rise to the partitions (123) and (23)(1). Similarly, the partition (3)(2) gives rise to the partitions (13)(2), (3)(21) and (3)(2)(1). We might characterize this process as computing new partitions of {1, 2, 3} from a partition p of {2, 3} as "extending" p .

Suppose then we write a function extend x p to capture the above idea. Let's start with the signature of extend. What would it be? Taking (23)(1) as an exemplar, we see that a component of a partition can be represented as [a] and so a partition itself then as [[a]]. We know that extend takes an element and a partition and returns a list of (new) partitions so it must have signature extend :: a -> [[a]] -> [[[a]]] (yes, lists of lists of lists are somehow easy to get confused about).

Now for writing the body of extend. The base case is the easiest of course - extending the empty partition:

extend x [] = [[[x]]]
  
That is, a singleton list of partitions where that one partition has one component. The inductive case is the partition obtained by "pushing" x into the first component of p together with the extensions that leave the first component of p alone.
extend x (h : tl) = ((x : h) : tl) : map (h :) (extend x tl)

We can now phrase the function partition with signature partition :: [a] -> [[[a]]] like this:

partition [] = [[]]
partition (h : tl) = concatMap (extend h) (partition tl)
The base case says, the only partition of the empty set is the the empty partition.

Wrapping it all up, the algorithm in entirety is

partition :: [a] -> [[[a]]]
partition [] = [[]]
partition (h : tl) = concatMap (extend h) (partition tl)
  where
    extend :: a -> [[a]] -> [[[a]]]
    extend x [] = [[[x]]]
    extend x (h : tl) = ((x : h) : tl) : map (h :) (extend x tl)

Sunday, June 10, 2018

Bucket Sort

Bucket Sort

Bucket sort assumes input generated by a random process that distributes elements uniformly over the interval [0, 1).

The idea of bucket sort is to divide [0, 1) into n equal-sized subintervals or buckets, and then distribute the n input numbers into the buckets. To produce the output, sort the numbers in each bucket and then go through the buckets in order. Sorting a bucket can be done with insertion sort.

let rec insert x = function
  | [] -> [x]
  | h :: tl as ls ->
    if x < h then x :: ls else h :: insert x tl

let rec insertion_sort = function
  | [] | [_] as ls -> ls
  | h :: tl -> insert h (insertion_sort tl)

This code for bucket sort assumes the input is an n-element array a and that each element 0 ≤ a.(i) < 1. The code requires an auxillary array b.(0 .. n - 1) of lists (buckets).

let bucket_sort a =
  let n = Array.length a in
  let b = Array.make n [] in
  Array.iter
    (fun x ->
       let i =
         int_of_float (
           floor (float_of_int n *. x)
         ) in
        Array.set b i (x :: Array.get b i)
      ) a;
  Array.iteri
    (fun i l ->
       Array.set b i (insertion_sort l)
    ) b;
  Array.fold_left (fun acc bucket -> acc @ bucket) [] b
;;
bucket_sort [| 0.78; 0.17; 0.39; 0.26; 0.72; 0.94
             ; 0.21; 0.12; 0.23; 0.68|]
Bucket sort runs in linear time on the average.


References:
[1] "Introduction to Algorithms" Section 9.4:Bucket Sort -- Cormen et. al. (Second ed.) 2001.

Sunday, May 20, 2018

Dijkstra's algorithm

Shortest Path

This article assumes familiarity with Dijkstra's shortest path algorithm. For a refresher, see [1]. The code assumes open Core is in effect and is online here.

The first part of the program organizes our thoughts about what we are setting out to compute. The signature summarizes the notion (for our purposes) of a graph definition in modular form. A module implementing this signature defines a type vertex_t for vertices, a type t for graphs and type extern_t : a representation of a t for interaction between an implemening module and its "outside world".

module type Graph_sig = sig
  type vertex_t [@@deriving sexp]
  type t [@@deriving sexp]
  type extern_t

  type load_error = [ `Duplicate_vertex of vertex_t ] [@@deriving sexp]
  exception Load_error of load_error [@@deriving sexp]

  val of_adjacency : extern_t -> [ `Ok of t | `Load_error of load_error ]
  val to_adjacency : t -> extern_t

  module Dijkstra : sig
    type state

    type error = [
      | `Relax of vertex_t
    ] [@@deriving sexp]
    exception Error of error [@@deriving sexp]

    val dijkstra : vertex_t -> t -> [ `Ok of state | `Error of error ]
    val d : state -> (vertex_t * float) list
    val shortest_paths : state -> (vertex_t * vertex_t list) list
  end

end
A realization of Graph_sig provides "conversion" functions of_adjacency/to_adjacency between the types extern_t and t and nests a module Dijkstra. The signature of the sub-module Dijkstra requires concrete modules provide a type state and an implementation of Dijkstra's algorithm in terms of the function signature val dijkstra : vertex_t -> t -> [ `Ok of state | `Error of error ].

For reusability, the strategy for implementing graphs will be generic programming via functors over modules implementing s vertex type.

An implementation of the module type GRAPH defines a module type VERT which is required to provide a comparable type t. It further defines a module type S that is exactly module type Graph_sig above. Lastly, modules of type GRAPH provide a functor Make that maps any module of type VERT to new module of type S fixing extern_t to an adjacency list representation in terms of the native OCaml type 'a list and float to represent weights on edges.

module type GRAPH = sig
  module type VERT = sig
    type t[@@deriving sexp]
    include Comparable.S with type t := t
  end

  module type S = sig
    include Graph_sig
  end

  module Make : functor (V : VERT) ->
    S with type vertex_t = V.t
       and type extern_t = (V.t * (V.t * float) list) list
end
The two module types Graph_sig and GRAPH together provide the specification for the program. module Graph in the next section implements this specification.

Implementation of module Graph is in outline this.

module Graph : GRAPH = struct
  module type VERT = sig
    type t[@@deriving sexp]
    include Comparable.S with type t := t
  end

  module type S = sig
    include Graph_sig
  end

  module Make : functor (V : VERT) ->
    S with type vertex_t = V.t
       and type extern_t = (V.t * (V.t * float) list) list
    =

    functor (V : VERT) -> struct
       ...
    end
end
As per the requirements of GRAPH the module types VERT and S are provided as is the functor Make. It is the code that is ellided by the ... above in the definition of Make that is now the focus.

Modules produced by applications of Make satisfy S. This requires suitable definitions of types vertext_t, t and extern_t. The modules Map and Set are available due to modules of type VERT being comparable in their type t.

      module Map = V.Map
      module Set = V.Set

      type vertex_t = V.t [@@deriving sexp]
      type t = (vertex_t * float) list Map.t [@@deriving sexp]
      type extern_t = (vertex_t * (vertex_t * float) list) list
      type load_error = [ `Duplicate_vertex of vertex_t ] [@@deriving sexp]
      exception Load_error of load_error [@@deriving sexp]

While the external representation extern_t of graphs is chosen to be an adjacency list representation in terms of association lists, the internal representation t is a vertex map of adjacency lists providing logarithmic loookup complexity. The conversion functions between the two representations "come for free" via module Map.

      let to_adjacency g = Map.to_alist g

      let of_adjacency_exn l =  match Map.of_alist l with
        | `Ok t -> t
        | `Duplicate_key c -> raise (Load_error (`Duplicate_vertex c))

      let of_adjacency l =
        try
          `Ok (of_adjacency_exn l)
        with
        | Load_error err -> `Load_error err

At this point the "scaffolding" for Dijkstra's algorithm, that part of GRAPH dealing with the representation of graphs is implemented.

The interpretation of Dijkstra's algorithm we adopt is functional : the idea is we loop over vertices relaxing their edges until all shortest paths are known. What we know on any recursive iteration of the loop is a current "state" (of the computation) and each iteration produces a new state. This next definition is the formal definition of type state.

      module Dijkstra = struct

        type state = {
          src    :                  vertex_t
        ; g      :                         t
        ; d      :               float Map.t
        ; pred   :            vertex_t Map.t
        ; s      :                     Set.t
        ; v_s    : (vertex_t * float) Heap.t
        }
The fields of this record are:
  • src : vertex_t, the source vertex;
  • g : t, G the graph;
  • d : float Map.t, d the shortest path weight estimates;
  • pre : vertex_t Map.t, π the predecessor relation;
  • s : Set.t, the set S of nodes for which the lower bound shortest path weight is known;
  • v_s : (vertex_t * float) Heap.t, V - {S}, , the set of nodes of g for which the lower bound of the shortest path weight is not yet known ordered on their estimates.

Function invocation init src g compuates an initial state for the graph g containing the source node src. In the initial state, d is everywhere except for src which is 0. Set S (i.e. s) and the predecessor relation π (i.e. pred) are empty and the set V - {S} (i.e. v_s) contains all nodes.

        let init src g =
          let init x = match V.equal src x with
            | true -> 0.0 | false -> Float.infinity in
          let d = List.fold (Map.keys g) ~init:Map.empty
              ~f:(fun acc x -> Map.set acc ~key:x ~data:(init x)) in
          {
            src
          ; g
          ; s = Set.empty
          ; d
          ; pred = Map.empty
          ; v_s = Heap.of_list (Map.to_alist d)
                ~cmp:(fun (_, e1) (_, e2) -> Float.compare e1 e2)
          }

Relaxing an edge (u, v) with weight w (u, v) tests whether the shortest path to v so far can be improved by going through u and if so, updating d (v) and π (v) accordingly.

        type error = [
          | `Relax of vertex_t
        ] [@@deriving sexp]
        exception Error of error [@@deriving sexp]

        let relax state (u, v, w) =
          let {d; pred; v_s; _} = state in
          let dv = match Map.find d v with
            | Some dv -> dv
            | None -> raise (Error (`Relax v)) in
          let du = match Map.find d u with
            | Some du -> du
            | None -> raise (Error (`Relax u)) in
          if dv > du +. w then
            let dv = du +. w in
            (match Heap.find_elt v_s ~f:(fun (n, _) -> V.equal n v) with
            | Some tok -> ignore (Heap.update v_s tok (v, dv))
            | None -> raise (Error (`Relax v))
            );
            { state with
              d = Map.change d v
                  ~f:(function
                      | Some _ -> Some dv
                      | None -> raise (Error (`Relax v))
                    )
            ; pred = Map.set (Map.remove pred v) ~key:v ~data:u
            }
          else state
Here, relaxation can result in a linear heap update operation. A better implementation might seek to avoid that.

One iteration of the body of the loop of Dijkstra's algorithm consists of the node in V - {S} with the least shortest path weight estimate being moved to S and its edges relaxed.

        let dijkstra_exn src g =
          let rec loop ({s; v_s; _} as state) =
            match Heap.is_empty v_s with
            | true -> state
            | false ->
              let u = fst (Heap.pop_exn v_s) in
              loop (
                List.fold (Map.find_exn g u)
                  ~init:{ state with s = Set.add s u }
                  ~f:(fun state (v, w) -> relax state (u, v, w))
              )
          in loop (init src g)

        let dijkstra src g =
          try
            `Ok (dijkstra_exn src g)
          with
          | Error err -> `Error err

The shortest path estimates contained by a value of state is given by the projection d.

        let d state = Map.to_alist (state.d)

The shortest paths themselves are easily computed as,

   let path state n =
          let rec loop acc x =
            (match V.equal x state.src with
            | true -> x :: acc
            | false -> loop (x :: acc) (Map.find_exn state.pred x)
            ) in
          loop [] n

        let shortest_paths state =
          List.map (Map.keys state.g) ~f:(fun n -> (n, path state n))
      end
    end
which completes the implementation of Make.

The following program produces a concrete instance of the shortest path problem (with some evaluation output from the top-level).

module G : Graph.S with
  type vertex_t = char and type extern_t = (char * (char * float) list) list
  =
  Graph.Make (Char)

let g : G.t =
  match G.of_adjacency
          [ 's', ['u',  3.0; 'x', 5.0]
          ; 'u', ['x',  2.0; 'v', 6.0]
          ; 'x', ['v',  4.0; 'y', 6.0; 'u', 1.0]
          ; 'v', ['y',  2.0]
          ; 'y', ['v',  7.0]
          ]
  with
  | `Ok g -> g
  | `Load_error e -> failwiths "Graph load error : %s" e G.sexp_of_load_error
;;
let s = match (G.Dijkstra.dijkstra 's' g) with
  | `Ok s -> s
  | `Error e -> failwiths "Error : %s" e G.Dijkstra.sexp_of_error

;; G.Dijkstra.d s
- : (char * float) list =
[('s', 0.); ('u', 3.); ('v', 9.); ('x', 5.); ('y', 11.)]

;; G.Dijkstra.shortest_paths s
- : (char * char list) list =
[('s', ['s']); ('u', ['s'; 'u']); ('v', ['s'; 'u'; 'v']); ('x', ['s'; 'x']);
 ('y', ['s'; 'x'; 'y'])]


References:
[1] "Introduction to Algorithms" Section 24.3:Dijkstra's algorithm -- Cormen et. al. (Second ed.) 2001.

Saturday, August 12, 2017

Transpose

Transpose

If we are to represent a row of a matrix as a list of numbers, then a matrix can naturally be represented as a list of lists of numbers.

The transpose of a matrix $\mathbf{A}$ is a new matrix denoted $\mathbf{A^{T}}$. The traditional mathematical definition of $\mathbf{A^{T}}$ is expressed as saying the $i$ th row, $j$ th column element of $\mathbf{A^{T}}$ is the $j$ th row, $i$ th column element of $\mathbf{A}$:

$\left[\mathbf{A}\right]_{ij} = \left[\mathbf{A^{T}}\right]_{ji}$.

As definitions go, this isn't terribly helpful in explaining how to compute a transpose. A better equivalent definition for the functional programmer is : the matrix obtained by writing the columns of $\mathbf{A}$ as the rows of $\mathbf{A^{T}}$.

An elegant program for computing a transpose follows from a direct translation of that last definition.

let rec transpose (ls : 'a list list) : 'a list list =
  match ls with
  | [] | [] :: _ -> []
  | ls -> List.map (List.hd) ls :: transpose (List.map (List.tl) ls)

It is not at all hard to understand how the program works when you've seen an example:

transpose [[1; 2]; [3; 4;]; [5; 6]]
  = [1; 3; 5] :: transpose [[2]; [4;]; [6]]
  = [1; 3; 5] :: [2; 4; 6] :: transpose [[]; []; []]
  = [1; 3; 5] :: [2; 4; 6] :: []
  = [[1; 3; 5]; [2; 4; 6]]

Being as pretty as it is, one is inclined to leave things be but, as a practical matter it should be rephrased to be tail-recursive.

let rec transpose (ls : 'a list list) : 'a list list  =
  let rec transpose_rec acc = function
  | [] | [] :: _ -> List.rev acc
  | ls -> transpose_rec (List.map (List.hd) ls :: acc) (List.map (List.tl) ls)
  in transpose_rec [] ls


References:
"An Introduction to Functional Programming Systems Using Haskell" -- Davie A J T., 1992

Sunday, April 3, 2016

C++ : Streams

In this blog post, types and functions were presented in OCaml for modeling streams. This post takes the action to C++.

First, the type definition for a stream.

struct Nil {};
template <class T> class Cons;

template <class T>
using stream = sum_type <
    Nil
  , recursive_wrapper<Cons<T>>
>;
The definition is in terms of the sum_type<> type from the "pretty good sum" library talked about here.

The definition of Cons<>, will be in terms of "thunks" (suspensions). They're modeled as procedures that when evaluated, compute streams.

template <class T>
using stream_thunk = std::function<stream<T>()>;
To complete the abstraction, a function that given a suspension, "thaws" it.
template <class T> inline 
stream<T> force (stream_thunk<T> const& s) { 
  return s (); 
}

The above choices made, here is the definition for Cons<>.

template <class T>
class Cons {
public:
  using value_type = T;
  using reference = value_type&;
  using const_reference = value_type const&;
  using stream_type = stream<value_type>;

private:
  using stream_thunk_type = stream_thunk<value_type>;

public:
  template <class U, class V>
  Cons (U&& h, V&& t) : 
    h {std::forward<U> (h)}, t {std::forward<V> (t)}
  {}

  const_reference hd () const { return h; }
  stream_type tl () const { return force (t); }

private:
  value_type h;
  stream_thunk_type t;
};

Next, utility functions for working with streams.

The function hd () gets the head of a stream and tl () gets the stream that remains when the head is stripped off.

template <class T>
T const hd (stream<T> const& s) {
  return s.template match<T const&> (
      [](Cons<T> const& l) -> T const& { return l.hd (); }
    , [](otherwise) -> T const & { throw std::runtime_error { "hd" }; }
  );
}

template <class T>
stream<T> tl (stream<T> const& l) {
  return l.template match <stream<T>> (
    [] (Cons<T> const& s) -> stream <T> { return s.tl (); }
  , [] (otherwise) -> stream<T> { throw std::runtime_error{"tl"}; }
  );
}

The function take () returns the the first $n$ values of a stream.

template <class T, class D>
D take (unsigned int n, stream <T> const& s, D dst) {
  return (n == 0) ? dst :
    s.template match<D>(
       [&](Nil const& _) -> D { return  dst; },
       [&](Cons<T> const& l) -> D { 
         return take (n - 1, l.tl (), *dst++ = l.hd ()); }
    );
}

It's time to share a little "hack" I picked up for writing infinite lists.

  • To start, forget about streams;
  • Write your list using regular lists;
  • Ignore the fact that it won't terminate;
  • Rewrite in terms of Cons and convert the tail to a thunk.

For example, in OCaml the (non-terminating!) code

  let naturals = 
    let rec loop x = x :: loop (x + 1) in
  next 0
leads to this definition of the stream of natural numbers.
let naturals =
 let rec loop x = Cons (x, lazy (loop (x + 1))) in
loop 0

Putting the above to work, a generator for the stream of natural numbers can be written like this.

class natural_numbers_gen {
private:
  using int_stream = stream<int>;
    
private:
  int start;

private:
  int_stream from (int x) const {
    return int_stream{
      constructor<Cons<int>>{}, x, [=]() { return this->from (x + 1); }
    };
  }
  
public:
  explicit natural_numbers_gen (int start) : start (start) 
  {}

  explicit operator int_stream() const { return from (start); }
};
The first $10$ (say) natural numbers can then be harvested like this.
std::vector<int> s;
take (10, stream<int> (natural_numbers_gen{0}), std::back_inserter (s));

The last example, a generator of the Fibonacci sequence. Applying the hack, start with the following OCaml code.

  let fibonacci_numbers = 
    let rec fib a b = a :: fib b (a + b) in
    fib 0 1
The rewrite of this code into streams then leads to this definition.
let fibonnaci_sequence = 
  let rec fib a b = Cons (a, lazy (fib b (a + b))) in
fib 0 1
Finally, casting the above function into C++ yields the following.
class fibonacci_numbers_gen {
private:
  using int_stream = stream<int>;
    
private:
  int start;

private:
  int_stream loop (int a, int b) const {
    return int_stream{
      constructor<Cons<int>>{}, a, [=]() {return this->loop (b, a + b); }
    };
  }
    
public:
  explicit fibonacci_numbers_gen () 
  {}

  explicit operator int_stream() const { return loop (0, 1); }
  };

Saturday, April 2, 2016

Rotate

Rotate

This post is inspired by one of those classic "99 problems in Prolog".What we are looking for here are two functions that satisfy these signatures.

val rotate_left : int -> α list -> α list
val rotate_right : int -> α list -> α list 
rotate_left n rotates a list $n$ places to the left, rotate_right n rotates a list $n$ places to the right. Examples:
# rotate_left 3 ['a';'b';'c';'d';'e';'f';'g';'h'] ;;
- : char list = ['d'; 'e'; 'f'; 'g'; 'h'; 'a'; 'b'; 'c']

# rotate_left (-2) ['a';'b';'c';'d';'e';'f';'g';'h'] ;;
- : char list = ['g'; 'h'; 'a'; 'b'; 'c'; 'd'; 'e'; 'f']
Of course, rotate_left and rotate_right are inverse functions of each other so we expect, for any int $x$ and list $l$, rotate_right x @@ rotate_left x l $=$ rotate_left x @@ rotate_right x l $=$ l.

Well, there are a variety of solutions to this problem with differing degrees of verbosity, complexity and efficiency. My own attempt at a solution resulted in this.

let rec drop (k : int) (l : α list) : α list =
  match k, l with
  | i, _ when i <= 0 -> l
  | _, [] -> []
  | _, (_ :: xs) -> drop (k - 1) xs

let rec take (k : int) (l : α list) : α list =
  match k, l with
  | i, _ when i <= 0 -> []
  | _, [] -> []
  | _, (x :: xs)  -> x :: take (k - 1) xs

let split_at (n : int) (l : α list) : α list * α list = 
  (take n l), (drop n l)

let rec rotate_left (n : int) (l : α list) : α list =
  match n with
  | _ when n = 0 -> l
  | _ when n < 0 ->  rotate_right (-n) l
  | _ -> 
    let m : int = List.length l in
    let k : int = n mod m in
    let (l : α list), (r : α list) = split_at k l in 
    r @ l

and rotate_right (n : int) (l : α list) : α list =
  match n with
  | _ when n = 0 -> l
  | _ when n < 0 ->  rotate_left (-n) l
  | _ -> 
    let m : int = List.length l in
    let k : int = m - n mod m in
    let (l : α list), (r : α list) = split_at k l in 
    r @ l

So far so good, but then I was shown the following solution in Haskell.

rotateLeft n xs 
  | n >= 0     = take (length xs) $ drop n $ concat $ repeat xs
  | otherwise  = rotateLeft (length xs + n) xs

rotateRight n = rotateLeft (-n)
I found that pretty nifty! See, in the function rotateLeft, repeat xs creates an infinite list of lists, (each a copy of xs), "joins" that infinite list of lists into one infinite list, then the first $n$ elements are dropped from that the list and we take the next length xs which gets us the original list rotated left $n$ places.

I felt compelled to attempt to emulate the program above in OCaml.

The phrasing "works" in Haskell due to the feature of lazy evaluation. OCaml on the other hand is eagerly evaluated. Lazy evaluation is possible in OCaml however, you just need to be explicit about it. Here's a type for "lazy lists" aka "streams".

type α stream =  Nil | Cons of α * α stream Lazy.t
A value of type α Lazy.t is a deferred computation, called a suspension that has the result type α. The syntax lazy$(expr)$ makes a suspension of $expr$, without yet evaluating $expr$. "Forcing" the suspension (using Lazy.force) evaluates $expr$ and returns its result.

Next up, functions to get the head and tail of a stream.

let hd = function | Nil -> failwith "hd" | Cons (h, _) -> h
let tl = function | Nil -> failwith "tl" | Cons (_, t) -> Lazy.force t
Also useful, a function to lift an α list to an α stream.
let from_list (l : α list) : α stream =
  List.fold_right (fun x s -> Cons (x, lazy s)) l Nil

Those are the basic building blocks. Now we turn attention to implementing repeat x to create infinite lists of the repeated value $x$.

let rec repeat (x : α) : α stream = Cons (x, lazy (repeat x))

Now to implement concat (I prefer to call this function by its alternative name flatten).

The characteristic operation of flatten is the joining together of two lists. For eager lists, we can write a function join that appends two lists like this.

let rec join l m =
  match l with
  | [] -> m
  | h :: t -> h :: (join t m)
This generalizes naturally to streams.
let rec join (l : α stream) (m : α stream) =
  match l with
  | Nil -> m
  | Cons (h, t) -> Cons (h, lazy (join (Lazy.force t) m))
For eager lists, we can write flatten in terms of join.
let rec flatten : α list list -> α list = function
   | [] -> []
   | (h :: tl) -> join h (flatten tl)
Emboldened by our earlier success we might try to generalize it to streams like this.
let rec flatten (l : α stream stream) : α stream =
   match l with
   | Nil -> lazy Nil
   | Cons (l, r) ->  join l (flatten (Lazy.force r))
Sadly, no. This definition is going to result in stack overflow. There is an alternative phrasing of flatten we might try.
let rec flatten = function
  | [] -> []
  | [] :: t -> flatten t
  | (x :: xs) :: t -> x :: (flatten (xs :: t))
Happy to say, this one generalizes and gets around the eager evaluation problem that causes the unbounded recursion.
let rec flatten : α stream stream -> α stream = function
  | Nil -> Nil
  | Cons (Nil, t) -> flatten (Lazy.force t)
  | Cons (Cons (x, xs), t) ->
      Cons (x, lazy (flatten (Cons (Lazy.force xs, t))))

take and drop are straight forward generalizations of their eager counterparts.

let rec drop (n : int) (lst : α stream ) : α stream = 
  match (n, lst) with
  | (n, _) when n < 0 -> invalid_arg "negative index in drop"
  | (n, xs) when n = 0 -> xs
  | (_, Nil) -> Nil
  | (n, Cons (_, t)) -> drop (n - 1) (Lazy.force t)

let rec take (n : int) (lst : α stream) : α list = 
  match (n, lst) with
  | (n, _) when n < 0 -> invalid_arg "negative index in take"
  | (n, _) when n = 0 -> []
  | (_, Nil) -> []
  | (n, Cons (h, t)) -> h :: (take (n - 1) (Lazy.force t))

Which brings us to the lazy version of rotate expressed in about the same number of lines of code!

let rec rotate_left (k : int) (l : α list) : α list =
  let n = List.length l in
  if k >= 0 then
    l |> from_list |> repeat |> flatten |> drop k |> take n
  else rotate_left (n + k) l

let rotate_right (n : int) : α list -> α list = rotate_left (-n)

Saturday, May 30, 2015

Run length encoding data compression method

Functional programming and lists go together like Fred and Ginger. This little exercise is one of Werner Hett's "Ninety-Nine Prolog Problems". The idea is to implement the run length encoding data compression method.

Here's how we start. First we write a function that packs consecutive duplicates of a list into sublists e.g.

# B.pack ['a'; 'a'; 'a'; 'b'; 'c'; 'c'; 'd'] ;;
- : char list list = [['a'; 'a'; 'a']; ['b']; ['c'; 'c']; ['d'];]
Then, consecutive duplicates of elements are encoded as terms $(N, E)$ where $N$ is the number of duplicates of the element $E$ e.g.
# B.rle (B.pack ['a'; 'a'; 'a'; 'b'; 'c'; 'c'; 'd'; 'e'; 'e']) ;;
- : (int * char) list = [(3, 'a'); (1, 'b'); (2, 'c'); (1, 'd'); (2, 'e')]
We will of course require a function to decode compressed data e.g.
 # B.unrle(B.rle(
     B.pack ['a'; 'a'; 'a'; 'b'; 'c'; 'c'; 'd'; 'e'; 'e'])
    ) ;;
- : char list = ['a'; 'a'; 'a'; 'b'; 'c'; 'c'; 'd'; 'e'; 'e']
(Credit goes to Harvey Stein for the names rle and unrle by the way).

So that's it for the first iteration - here's some code that aims to implement these specifications.
module B = struct

let pack (x : α list) : α list list =
  let f (acc : α list list) (c : α) : α list list =
    match acc with
    | (((b :: _) as hd) :: tl) when c = b -> (c :: hd) :: tl
    |  _ -> [c] :: acc
  in List.fold_left f [] x

let rle (x : α list list) : (int * α) list =
  let f (acc : (int * α) list) (l : α list) : (int * α) list =
    (List.length l, List.hd l) :: acc
  in List.fold_left f [] x
  
let unrle (data : (int * α) list) =
  let repeat ((n : int), (c : α)) : α list =
    let rec aux acc i = if i = 0 then acc else aux (c :: acc) (i - 1) in
    aux [] n in
  let f (acc : α list) (elem : (int * α)) : α list =
    acc @ (repeat elem) in
  List.fold_left f [] data
  
end

Now, pack is just a device of course. We don't really need it so here's the next iteration that does away with it.

module E = struct

let rle (x : α list) : (int * α) list =
  let f (acc : (int * α) list) (c : α) : (int * α) list =
    match acc with
    | ((n, e) :: tl) when e = c -> (n + 1, c):: tl
    | _-> (1, c) :: acc
  in List.rev (List.fold_left f [] x)

let unrle (data : (int * α) list) =
  let repeat ((n : int), (c : α)) : α list =
    let rec aux acc i = if i = 0 then acc else aux (c :: acc) (i - 1) in
    aux [] n in
  let f (acc : α list) (elem : (int * α)) : α list =
    acc @ (repeat elem) in
  List.fold_left f [] data

end

Nifty!

Ok, the next idea is that when a singleton byte is encountered, we don't write the term $(1, E)$ instead, we just write $E$. Now OCaml doesn't admit heterogenous lists like Prolog appears to do so we need a sum type for the two possibilities. This then is the final version.

module F = struct

  type α t = | S of α | C of (int * α)

  let rle (bytes : α list) : α t list =
    let f (acc : α t list) (b : α) : α t list =
      match acc with
      | ((S e) :: tl) when e = b -> (C (2, e)) :: tl
      | ((C (n, e)) :: tl) when e = b -> (C (n + 1, b)) :: tl
      | _-> S b :: acc
    in List.rev (List.fold_left f [] bytes)

  let unrle (data : (α t) list) =
    let rec aux (acc : α list) (b : α) : (int -> α list) = 
      function | 0 -> acc | i -> aux (b :: acc) b (i - 1) in
    let f (acc : α list) (e : α t) : α list =
      acc @ (match e with | S b -> [b]| C (n, b) -> aux [] b n) in
    List.fold_left f [] data

end

Having worked out the details in OCaml, translation into C++ is reasonably straight-forward. One economy granted by this language is that we can do away with the data constructor S in this version.

#include <boost/variant.hpp>
#include <boost/variant/apply_visitor.hpp>
#include <boost/range.hpp>
#include <boost/range/numeric.hpp>

#include <list>

//Representation of the encoding
template <class A> struct C { std::pair <int, A> item; };
template <class A> using datum = boost::variant <A, C<A>>;
template <class A> using encoding = std::list<datum<A>>;

//Procedural function object that updates an encoding given a
//datum
template <class A>
struct update : boost::static_visitor<> {
  A c;
  encoding<A>& l;

  update (A c, encoding<A>& l) : c (c), l (l) 
  {}

  void operator ()(A e) const { 
    if (e == c) {
      l.back () = C<A>{ std::make_pair(2, c) }; 
      return;
    }
    l.push_back (c);
  }

  void operator ()(C<A> const& elem) const { 
    if (elem.item.second == c) {
      l.back () = C<A>{ std::make_pair (elem.item.first + 1, c) };
      return;
    }
    l.push_back (c);
  }
};

template <class R>
encoding<typename boost::range_value<R>::type> rle (R bytes) {
  typedef boost::range_value<R>::type A;

  auto f = [](encoding<A> acc, A b) -> encoding<A> {
    if (acc.size () == 0)
      acc.push_back (b);
    else   {
      boost::apply_visitor (update<A>(b, acc), acc.back ());
    }
    return acc;
  };

  return boost::accumulate (bytes, encoding<A> (), f);
}

I've left implementing unrle () as an exercise. Here's a little test though that confirms that we are getting savings in from the compression scheme as we hope for.

int main () {
  std::string buf=
    "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"
    "bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbb"
    "c"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "ddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd"
    "eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee"
    "z";
  std::list<char> data(buf.begin (), buf.end());
  encoding<char> compressed = rle (data);

  std::cout << sizeof (char) * (data.size ()) << std::endl;
  std::cout << sizeof (datum <char>) * (compressed.size ()) << std::endl;

  return 0;
}
On my machine, this program prints the values $484$ and $72$.

Friday, March 6, 2015

Heap sort

Given the existence of a priority queue data structure, the heap sort algorithm is trivially implemented by loading the unsorted sequence into a queue then successively pulling of the minimum element from the queue until the queue is exhausted.

There are many ways to implement a priority queue and so we seek an expression for a function for heap sorting that is polymorphic over those choices.

To begin, a module type for priority queues.

(**Priority queues over ordered types*)
module type PRIORITY_QUEUE = sig

    (**Output signature of the functor [Make]*)
    module type S = sig
        exception Empty

        type element (*Abstract type of elements of the queue*)
        type t (*Abstract type of a queue*)

        val empty : t (*The empty queue*)
        val is_empty : t -> bool (*Check if queue is empty*)
        val insert : t -> element -> t (*Insert item into queue*)
        val delete_min : t -> t (*Delete the minimum element*)
        val find_min : t -> element (*Return the minimum element*)
        val of_list : element list -> t
      end

    (**Input signature of the functor [Make]*)
    module type Ordered_type = sig
        type t
        val compare : t -> t -> int
      end

   (**Functor building an implementation of the priority queue structure
      given a totally ordered type*)
    module Make : 
       functor (Ord : Ordered_type) -> S with type element = Ord.t
  end

An implementation of this signature using "leftist heaps" is described for the interested in this Caltech lab but such details are omitted here.

module Priority_queue : PRIORITY_QUEUE = struct 
  module type S = sig .. end
  module type Ordered_type = sig .. end
  module Make (Elt : Ordered_type) : (S with type element = Elt.t) = struct .. end
end

What I really want to show you is this. We start with the following module type abbreviation.

type 'a queue_impl = (module Priority_queue.S with type element = 'a)
Then, the heap_sort function can be written such that it takes a module as a first class value and uses a locally abstract type to connect it with the element type of the list to be sorted.
let heap_sort (type a) (queue : a queue_impl) (l : a list) : a list =
  let module Queue = 
     (val queue : Priority_queue.S with type element = a) in
  let rec loop acc h =
    if Queue.is_empty h then acc
    else
      let p = Queue.find_min h in
      loop (p :: acc) (Queue.delete_min h) in
  List.rev (loop [] (Queue.of_list l))
There we have it. The objective has been achieved : we have written a heap sorting function that is polymorphic in the implementation of the priority queue with which it is implemented.

Usage (testing) proceeds as in this example.

(*Prepare an [Priority_queue.Ordered_type] module to pass as argument
  to [Priority_queue.Make]*)
module Int : Priority_queue.Ordered_type with type t = int = struct 
  type t = int let compare = Pervasives.compare 
end

(*Make a priority queue module*)
module Int_prioqueue : (Priority_queue.S with type element = int) = Priority_queue.Make (Int)

(*Make a first class value of the module by packing it*)
let queue = (module Int_prioqueue : Priority_queue.S with type element = int)

(*Now, pass the module to [heap_sort]*)
let sorted = heap_sort queue [-1; -2; 2] (*Produces the list [-2; -1; 2]*)

Addendum :

These ideas can be pushed a little further yielding a simpler syntax for the parametric heapsort algorithm.
(*Type abbreviations*)

type 'a order_impl = (module Priority_queue.Ordered_type with type t = 'a)
type 'a queue_impl = (module Priority_queue.S with type element = 'a)

(*Module factory functions*)

let mk_ord : 'a. unit -> 'a order_impl =
  fun (type s) () ->
    (module 
     struct 
       type t = s 
       let compare = Pervasives.compare 
     end : Priority_queue.Ordered_type with type t = s
    )

let mk_queue : 'a. unit -> 'a queue_impl =
  fun (type s) ord ->
    let module Ord = 
      (val mk_ord () : Priority_queue.Ordered_type with type t = s) in
    (module Priority_queue.Make (Ord) :
                               Priority_queue.S with type element = s)
For example, now we can write
# heap_sort (mk_queue ()) [-3; 1; 5] ;;
  - : int list = [-3; 1; 5]

Sunday, February 15, 2015

Fold left via fold right

The puzzle is to express fold_left entirely in terms of fold_right. For example, an attempted solution like

let fold_left f e s =
  List.rev (fold_right (fun a acc -> f acc a) e (List.rev s))
is inadmissible because it relies on List.rev and thus is not entirely in terms of fold_right.

Recall that given a function $f$, a seed $e$ and a list $[a_{1}; a_{2}; \cdots; a_{N}]$, fold_left computes $f\;(\cdots f\;(f\;e\;a_{1})\;a_{2}\;\cdots)\;a_{N}$ whereas fold_right computes $f\;a_{1}\;(f\;a_{2}\;(\cdots\;(f\;a_{N}\;e)\cdots))$. There's really no choice but to delay computation and the expression that solves this problem is this.

let fold_left f e s =
  List.fold_right (fun a acc -> fun x -> acc (f x a)) s (fun x -> x) e

For example, in the top-level

 # fold_left (fun acc x -> x * x :: acc) [] [1; 2; 3;] ;;
   - : int list = [9; 4; 1]

To see how this works, consider the right fold over the list $[a_{1}; a_{2}; a_{3}]$ (say)

  • On encountering $a_{3}$ we compute $f_{3} = \lambda x_{3} . i\; (f\;x_{3}\;a_{3})$;
  • On encountering $a_{2}$ we compute $f_{2} = \lambda x_{2} . f_{3}\;(f\;x_{2}\;a_{2})$;
  • On encountering $a_{1}$ we compute $f_{1} = \lambda x_{1} . f_{2}\;(f\;x_{1}\;a_{1})$;
but then, $f_{1}\;e = f_{2}\;(f\;e\;a_{1}) = f_{3}\;(f\;(f\;e\;a_{1})\;a_{2}) = f\;(f\;(f\;e\;a_{1})\;a_{2})\;a_{3}$ as desired.

Sunday, November 9, 2014

Complexity

Complexity

When you begin to program in OCaml, one of the first pieces of advice you encounter is to prefer the '::' operator over the '@' operator for list construction. The question is, does it matter? Even knowing that '@' exhibits complexity $O(N)$ as opposed to $O(1)$ for '::' should you care? I mean, how much of a difference can it make really?

The answer of course is yes. In practical terms, it matters. No, it really, really matters. Let's quantify it. Here's a version of the range function that uses '@'.

let range s e =
  let rec loop acc s e =
    if s >= e then acc
    else loop (acc @ [s]) (s + 1) e 
  in loop [] s e
As an aside, you'll note that it's been written to be tail-recursive (so as to avoid inducing stack overflow).

Timing this function for building lists of $2^{10} = 1,024$ elements to $2^{16} = 65,536$ elements produces the following table:

ntime (s)
100.011243
110.047308
120.197137
130.866350
144.130998
1522.769691
16142.506546

Ouch! To build a list of just $65,536$ elements, the program takes over 2 minutes!?!

Ok, here's the version of range that uses '::' to build the list (and necessarily does a List.rev on the result in order that the elements don't come out "back to front"):

let range s e =
  let rec loop acc s e =
    if s >= e then acc
    else loop (s :: acc) (s + 1) e 
  in List.rev (loop [] s e)

And the timings for this one?

ntime (s)
100.000037
110.000097
120.000143
130.000324
140.002065
150.001195
160.005341

That is, this one builds the list of $65,536$ elements in just $5$ milliseconds!

Convinced?

Sunday, September 21, 2014

Correlation Coefficient

Correlation coefficient

It's been a while since we looked at anything from the domain of statistics so here's another little bite-sized piece - a function to compute Pearson's "product moment correlation coefficient".

It's a measure of dependence between two data sets. We'll express it in terms of unbiased standard deviation which I didn't write out before so I'll include that function too.

let unbiased_standard_deviation t =
  (*http://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation

    In statistics and in particular statistical theory, unbiased
    estimation of a standard deviation is the calculation from a
    statistical sample of an estimated value of the standard deviation
    (a measure of statistical dispersion) of a population of values,
    in such a way that the expected value of the calculation equals
    the true value.

  *)
 let av = arithmetic_mean t in
 let squared_diffs =
   List.fold_left (fun acc xi -> ((xi -. av) *. (xi -. av)) :: acc) [] t
 in sqrt ((sum squared_diffs)/.((float_of_int (List.length t)) -. 1.0))

let correlation_coefficient x y =
  (*http://en.wikipedia.org/wiki/Correlation_and_dependence

    The most familiar measure of dependence between two quantities is
    the Pearson product-moment correlation coefficient, or "Pearson's
    correlation coefficient", commonly called simply "the correlation
    coefficient". It is obtained by dividing the covariance of the two
    variables by the product of their standard deviations.  
  *) 

  let x_b = arithmetic_mean x in
  let y_b = arithmetic_mean y in
  let s_x = unbiased_standard_deviation x in
  let s_y = unbiased_standard_deviation y in

  if s_x = 0. || s_y = 0. then 0.
  else
    let f acc x_i y_i =
      acc +. ((x_i -. x_b) *. (y_i -. y_b)) in
    let n = float_of_int (List.length x) in
    let s = List.fold_left2 f 0.0 x y  in
    s/.((n -. 1.) *. s_x *. s_y)

Saturday, August 9, 2014

Digits

Digits

Counting is a common problem in computer programming. I recently had a need for an algorithm to compute the number of base 10 digits of an integer (for example, there are 3 digits in the number 569). Here's what I worked out.
module type S=sig
    (**Count the number of base 10 digits of an integer*)
    val digits : int -> int 

    (**Prepend with zeroes as necessary to get the 'right'
       number of columns *)
    val digits_10 : int -> int -> string

end

module A : S = struct 
  let digits n =
    let rec loop acc n =
      let k = n / 10 in
      if k = 0 then acc else loop (acc + 1) k
    in loop 1 n

  let digits_10 i n =
    let l = digits n in
    let s = string_of_int (abs n) in
    let maybe_sign = if n < 0 then "-" else "" in
    if l >= i then (maybe_sign^s) 
    else
      let lz=string_of_int 0 in
      let rec aux acc j = 
        if j = 0 then acc 
        else 
          aux (lz^acc) (j - 1) in
      maybe_sign ^ aux s (i - l)
end

module Test_io = struct

  let test_0 () = 
    let print_line i = 
      Printf.printf "There are %d digits in %d\n" (A.digits i) i
    in
    let rec loop i n =
      if i = n then () 
      else
        let x  =
          int_of_float (10. ** (float_of_int i)) in
        print_line x; loop (i + 1) n
    in (*Range over 0 to to 10^19 *) loop 0 19

  let test_1 () =
    let print_line i = 
      Printf.printf "%s\n" i
    in
    let rec loop i n =
      if i = n then () 
      else
        let x = A.digits_10 3 i in
        print_line x; loop (i + 1) n
    in (*Range over 0 to 100, 3 columns*)loop 0 100

  let run () = test_0 (); test_1 ()

end

let () = Test_io.run ()
The program outputs
There are 1 digits in 1
There are 2 digits in 10
There are 3 digits in 100
There are 4 digits in 1000
There are 5 digits in 10000
There are 6 digits in 100000
There are 7 digits in 1000000
There are 8 digits in 10000000
There are 9 digits in 100000000
There are 10 digits in 1000000000
There are 11 digits in 10000000000
There are 12 digits in 100000000000
There are 13 digits in 1000000000000
There are 14 digits in 10000000000000
There are 15 digits in 100000000000000
There are 16 digits in 1000000000000000
There are 17 digits in 10000000000000000
There are 18 digits in 100000000000000000
There are 19 digits in 1000000000000000000
000
001
002
003
004
005
006
 .
 .
 .
097
098
099

Fun, fun, fun!

Friday, July 25, 2014

Cartesian product (redux)

Cartesian product (redux)

In this blog post we looked at programs to compute Cartesian Products. One algorithm (given here in OCaml) if you recall is
let prod l r =
  let g acc a =
    let f acc x =
      (a, x) :: acc
    in List.fold_left f acc r
  in List.fold_left g [] l |> List.rev

In that earlier post I translated the above program into C++. This time I'm doing the same straightforward translation but using the Boost.MPL library to compute such products at compile time. It looks like this:

#include <boost/mpl/pair.hpp>
#include <boost/mpl/vector.hpp>
#include <boost/mpl/push_front.hpp>
#include <boost/mpl/fold.hpp>
#include <boost/mpl/reverse.hpp>
#include <boost/mpl/placeholders.hpp>

using namespace boost::mpl::placeholders;

template <class L, class R>
struct product
{
  struct g {
    template <class AccT, class A>
    struct apply {
       typedef typename boost::mpl::fold <
         R, AccT
        , boost::mpl::push_front<_1, boost::mpl::pair<A, _2> > 
         >::type type;
    };
  };

  typedef typename boost::mpl::reverse <
    typename boost::mpl::fold <
                  L, boost::mpl::vector<>, g>::type>::type type;
};
The translation proceeds almost mechanically! Does it work though? You betcha! Here's a test we can convince ourselves with:
#include <boost/mpl/equal.hpp>
#include <boost/mpl/for_each.hpp>
#include <boost/mpl/int.hpp>

#include <iostream>

using namespace boost::mpl;

typedef vector<int_<1>, int_<2> > A;
typedef vector<int_<3>, int_<4>, int_<5> > B;
typedef product <A, B>::type result;

BOOST_MPL_ASSERT((
  equal<
    result
  , vector<
         pair<int_<1>, int_<3> >
       , pair<int_<1>, int_<4> >
       , pair<int_<1>, int_<5> >
       , pair<int_<2>, int_<3> >
       , pair<int_<2>, int_<4> >
       , pair<int_<2>, int_<5> >
      > 
  > ));

struct print_class_name {
    template <typename T>
    void operator()( T t ) const {
       std::cout << typeid(t).name() << " ";
    }
};

int main ()
{
  boost::mpl::for_each<result>(print_class_name());

  return 0;
}
The takeaway is, learning yourself some functional programming is a great way to improve your template meta-programming fu! (That of course should come as no surprise... )

Friday, July 18, 2014

Merge sort

Merge sort

Here's another divide and conquer algorithm. This one is for sorting a sequence.

Conceptually, a merge sort works like this (see http://en.wikipedia.org/wiki/Merge_sort):

  • Divide the unsorted list into n sub-lists, each containing one element (a list of one element is trivially sorted);
  • Repeatedly merge sublists to produce new sorted sub-lists until there is only one sub-list remaining : this will be the sorted list.

In the following OCaml, we are drawing on inspiration from the the Haskell standard prelude for the part of the algorithm concerned with dividing the unsorted list : functions take, drop and split.

(**{b Merge sort}, an {i O(n log n)} comparison based sorting
   algorithm *)
module type MERGESORT = sig

  (**[take n] applied to a list [xs], returns the prefix of [xs] of
     length [n] or [xs] itself if [n > List.length xs] e.g. [take 2
     [1; 2; 3]] {i = } [[1; 2]]*)
  val take : int -> 'a list -> 'a list

  (**[drop n] applied to a list [xs], returns the suffix of [xs]
     after the first [n] elements or, [[]] if [n > List.length xs]
     e.g. [drop 2 [1; 2; 3]] {i = } [[3]]*)
  val drop : int -> 'a list -> 'a list

  (**[split n xs] is equivalent to [take n xs, drop n xs] e.g.
     [split 2 [1; 2; 3]] {i = } [([1; 2], [3])]*)
  val split : int -> 'a list -> 'a list * 'a list

  (**[merge] given two {b sorted} sequences [xs] and [ys] returns a
     single sorted sequence of the elements of [xs] and [ys]
     e.g. [merge [1; 3] [2; 4]] {i = } [[1; 2; 3; 4]]*)
  val merge : 'a list -> 'a list -> 'a list

  (**[merge_sort] works by splitting a sequence into two parts,
     recursively sorting the two parts and merging the results into
     a single sorted sequence e.g. [merge_sort [1; 2; -1; 0; 3]] 
     {i = } [[-1; 0; 1; 2; 3]]*)
  val merge_sort : 'a list -> 'a list
end

module Merge_sort : MERGESORT = struct

  let rec take k l =
    match (k, l) with
    | n, _ when n <= 0 -> []
    | _, [] -> []
    | n, (x :: xs) -> x :: take (n - 1) xs

  let rec drop k l =
    match (k, l) with
    | n, xs when n <= 0 -> xs
    | _, [] -> []
    | n, (_ :: xs) -> drop (n - 1) xs

  let rec split k l = take k l, drop k l

  let rec merge l m =
    match (l, m) with
    | [], ys -> ys
    | xs, [] -> xs
    | ((x :: xs) as t), ((y :: ys) as s) -> 
      if x <= y then x :: (merge xs s) else y :: (merge t ys)
        
  let rec merge_sort l =
    let i = (List.length l) / 2 in
    if i = 0 then l
    else
      let u, v = split i l in
      let xs, ys = merge_sort u, merge_sort v in
      merge xs ys

end

We can test our little program in the top-level like this:

# let l = Merge_sort.merge_sort [-1; 2; 0; 4; 3; 5];;
val l : int list = [-1; 0; 2; 3; 4; 5]

Here are the functions in C++ phrased as range algorithms.

//"c:/users/fletch/desktop/setenv.cmd"
//cl /EHsc /Femerge.exe /I c:/project/boost_1_55_0 merge.cpp 

#include <boost/next_prior.hpp>
#include <boost/range.hpp>
#include <boost/range/algorithm/copy.hpp>

#include <vector>
#include <cstdlib>

namespace algo
{
  //take

  template <class S, class D>
  D take (std::size_t n, S src, D dst)
  {
    typedef boost::range_iterator<S>::type it_t;
    it_t curr = boost::begin (src), last = boost::end (src);
    if (n <= 0) return dst;
    if (boost::empty (src))  return dst;
    take (n-1, S (boost::next (curr), last), *dst++ = *curr);
    return dst;
  }

  //drop

  template <class S, class D>
  D drop (std::size_t n, S src, D dst)
  {
    typedef boost::range_iterator<S>::type it_t;
    it_t curr = boost::begin (src), last = boost::end (src);
    if (n <= 0) return boost::range::copy (src, dst);
    if (boost::empty (src))return dst;
    drop (n-1, S (boost::next (curr), last), dst);
    return dst;
  }

  //split

  template <class S, class D1, class D2>
  void split (int n, S src, D1 dst1, D2 dst2)
  { take (n, src, dst1); drop (n, src, dst2); }

  //merge

  template <class S1, class S2, class D>
  D merge (S1 src1, S2 src2, D dst)
  {
    typedef boost::range_iterator<S1>::type it1_t;
    typedef boost::range_iterator<S2>::type it2_t;
    typedef std::iterator_traits<it1_t>::value_type value1_t;
    typedef std::iterator_traits<it2_t>::value_type value2_t;
    if (boost::empty (src1)) return boost::copy (src2, dst);
    if (boost::empty (src2)) return boost::copy (src1, dst);
    it1_t curr1 = boost::begin (src1), last1 = boost::end (src1);
    it2_t curr2 = boost::begin (src2), last2 = boost::end (src2);
    value1_t x = *curr1;
    value2_t y = *curr2;
    if (x <= y)
      merge (S1 (boost::next (curr1), last1), src2, *dst++ = x);
    else
      merge (src1, S2 (boost::next (curr2), last2), *dst++ = y);
    return dst;
  }

  //merge_sort

  template <class S, class D>
  D merge_sort (S src, D dst)
  {
    typedef boost::range_iterator<S>::type it_t;
    typedef std::iterator_traits<it_t>::value_type value_t;
    std::size_t i = boost::size (src)/2;
    if (i == 0) return boost::range::copy (src, dst);
    std::vector<value_t> u, v;
    split (i, src, std::back_inserter (u), std::back_inserter (v));
    std::vector<value_t> xs, ys;
    merge_sort (u, std::back_inserter (xs));
    merge_sort (v, std::back_inserter (ys));
    merge (xs, ys, dst);
    return dst;
  }
  
}//namespace algo

The following usage example provides a lightweight test.

#include <boost/assign/list_of.hpp>

#include <utility>
#include <iterator>
#include <iostream>

int main ()
{
  using std::pair;
  using std::make_pair;
  using std::ostream_iterator;
  using boost::assign::list_of;

  int ord[] = {1, 2, 3, 4};

  auto src=make_pair(ord, ord + 4);
  auto dst=ostream_iterator<int>(std::cout, ", ");

  std::cout << "\ntake ():\n";

  algo::take (0u, src, dst); std::cout << std::endl;
  algo::take (1u, src, dst); std::cout << std::endl;
  algo::take (2u, src, dst); std::cout << std::endl;
  algo::take (3u, src, dst); std::cout << std::endl;
  algo::take (4u, src, dst); std::cout << std::endl;
  algo::take (5u, src, dst); std::cout << std::endl;

  std::cout << "\ndrop ():\n";

  algo::drop (5u, src, dst); std::cout << std::endl;
  algo::drop (4u, src, dst); std::cout << std::endl;
  algo::drop (3u, src, dst); std::cout << std::endl;
  algo::drop (2u, src, dst); std::cout << std::endl;
  algo::drop (1u, src, dst); std::cout << std::endl;
  algo::drop (0u, src, dst); std::cout << std::endl;

  std::cout << "\nsplit ():\n\n";

  algo::split (0u, src, dst, dst); std::cout << std::endl;
  algo::split (1u, src, dst, dst); std::cout << std::endl;
  algo::split (2u, src, dst, dst); std::cout << std::endl;
  algo::split (3u, src, dst, dst); std::cout << std::endl;
  algo::split (4u, src, dst, dst); std::cout << std::endl;
  algo::split (5u, src, dst, dst); std::cout << std::endl;

  std::cout << "\nmerge ():\n\n";

  std::vector <int> l=list_of(-1)(2), r=list_of(0)(1);
  algo::merge (l, r, dst); std::cout << std::endl;

  std::cout << "\nmerge_sort ():\n\n";

  int unord[] = {-1, 2, 0, 4, 3, 5};
  algo::merge_sort (make_pair (unord, unord + 6), dst);

  return 0;
}
The above program produces the following output.
take ():

1,
1, 2,
1, 2, 3,
1, 2, 3, 4,
1, 2, 3, 4,

drop ():


4,
3, 4,
2, 3, 4,
1, 2, 3, 4,

split ():

1, 2, 3, 4,
1, 2, 3, 4,
1, 2, 3, 4,
1, 2, 3, 4,
1, 2, 3, 4,
1, 2, 3, 4,

merge ():

-1, 0, 1, 2,

merge_sort ():

-1, 0, 2, 3, 4, 5,

Saturday, June 21, 2014

Permutations

Permutations

Combinations are selections that disregard order. The powerset P(S) of a set S is the set of all possible combinations of the n elements of S. There are 2n of these. Permutations are arrangements of sequences into orders. The number of permutations of a set is denoted n! For example, the set {x, y, z} has 6 permutations : {x, y, z}, {x, z, y}, {y, x, z}, {y, z, x}, {z, y, x} and {z, x, y}. Suppose f exists that for given k provides all permutations of S that start with sk. All permutations could then be found by f with k ranging over [0, N - 1]. But f trivially exists as it can be computed by prepending sk to the permutations of the smaller set S - { sk }. Here are two programs to compute the permutations of a set based on this idea.

Thursday, June 19, 2014

Powerset

Powerset

If S is a set, then P(S), the 'powerset' of S is the set of all subsets of S including the empty set and S itself. If S has cardinality N then the cardinality of P(S) is 2N (why?). That is, there are 2N subsets associated with S.

Here's a function to compute P(S) in OCaml. It's an instance of the 'divide and conquer' strategy of problem solving.

  let rec sets l =
    match l with
    | [] -> [[]]
    | x :: xs -> let l = sets xs in 
                   l @ (List.map (fun y -> x :: y) l)

This program translates to C++ naturally and with brevity thanks to C++ 11 lambdas.

#include <boost/utility.hpp>

#include <set>
#include <iterator>
#include <algorithm>

template <class I, class D>
D sets (I begin, I end, D dst)
{
  typedef typename std::iterator_traits<I>::value_type value_type;
  typedef std::set<value_type> set_t;

  if (begin == end)
    {
      *dst++ = set_t (); //the empty set
    }
  else
    {
      std::set<set_t> l;
      std::set<set_t>::iterator back=l.end ();
      sets (boost::next (begin), end, std::inserter (l, back));

      std::transform (l.begin (), l.end (), dst, 
        [](set_t const& s) -> set_t const& { return s; });
      std::transform (l.begin (), l.end (), dst, 
        [=](set_t s) -> set_t { s.insert (*begin); return s; });
    }

  return dst;
}

Thursday, May 15, 2014

Depth first search

Depth first search

Depth first search is an 'elementary graph algorithm'. This is a purely functional formulation.

(*depth-first-search

 "Introduction to Algorithms" - Cormen et. al., 1994

*)
module Char_map = Map.Make (Char)

type graph = (char list) Char_map.t

module type S = sig
  type state
  val string_of_state : state -> string
  val depth_first_search : graph -> state
end

module Dfs : S = struct

  type colors = White|Gray|Black

  type state = {
    d : int Char_map.t ; (*discovery time*)
    f : int Char_map.t ; (*finishing time*)
    pred : char Char_map.t ; (*predecessor*)
    color : colors Char_map.t ; (*vertex colors*)
  }

  let string_of_state {d; f; pred; color} =
    let open Printf in
    let bindings m fmt =
      let b = Char_map.bindings m in
      String.concat ", " (List.map (fun (x,y) -> sprintf fmt x y) b) in
    sprintf " d = {%s}\n f = {%s}\n pred = {%s}\n"
      (bindings d "'%c':'%d'") (bindings f "'%c':'%d'")
      (bindings pred "'%c':'%c'")

  let depth_first_search g =
    let node u (t, {d; f; pred; color}) =
      let rec dfs_visit t u {d; f; pred; color} =
        let edge (t, {d; f; pred; color}) v =
          if Char_map.find v color = White then
            dfs_visit t v {d; f; pred=(Char_map.add v u pred); color}
          else  (t, {d; f; pred; color})
        in
        let t, {d; f; pred; color} =
          let t = t + 1 in
          List.fold_left edge
            (t, {d=Char_map.add u t d; f;
                 pred; color=Char_map.add u Gray color})
            (Char_map.find u g)
        in
        let t = t + 1 in
        t , {d; f=(Char_map.add u t f); pred; color=Char_map.add u Black color}
      in
      if Char_map.find u color = White then dfs_visit t u {d; f; pred; color}
      else (t, {d; f; pred; color})
    in
    let v = List.fold_left (fun acc (x, _) -> x::acc) [] (Char_map.bindings g) in
    let initial_state= 
       {d=Char_map.empty;
        f=Char_map.empty;
        pred=Char_map.empty;
        color=List.fold_right (fun x->Char_map.add x White) v Char_map.empty}
    in
    snd (List.fold_right node v (0, initial_state))

end

(* Test *)

let () =
  let g =
       List.fold_right
          (fun (x, y) -> Char_map.add x y)
          ['u', ['v'; 'x'] ;
           'v',      ['y'] ;
           'w', ['z'; 'y'] ;
           'x',      ['v'] ;
           'y',      ['x'] ;
           'z',      ['z'] ;
          ]
          Char_map.empty
  in
  let s = Dfs.depth_first_search g in
  Printf.printf "%s\n" (Dfs.string_of_state s)

Saturday, March 29, 2014

Power means

Power means

Before we get going, let me mention that my friends over at What does the quant say? have started a blog! This weeks post is about? You guessed it, statistics! You can check out "Getting our hands dirty with beta functions" by clicking here.

This post picks up where we left off in "Statistics", looking at different definitions for the mean of a sequence of numbers.

The geometric mean together with the arithmetic and harmonic means of the last post together are called the Pythagorean means.

let geometric_mean t =
  (*http://en.wikipedia.org/wiki/Geometric_mean

    The geometric mean is defined as the nth root (where n is the
    count of numbers) of the product of the numbers.

  *)
  let n = List.length t in
  let prod = List.fold_left (fun acc xi -> acc *. xi) 1.0 t in
  prod ** (1.0/.(float_of_int n))
Note the presence of the power function powf (operator ( ** ))! How to compute values of this function was the focus of this post.

Here is yet another common mean, the quadratic mean.

let quadratic_mean t =
  (*http://en.wikipedia.org/wiki/Root_mean_square

  In mathematics, the root mean square (abbreviated RMS or rms),
  also known as the quadratic mean, is a statistical measure of the
  magnitude of a varying quantity. It is especially useful when
  variates are positive and negative, e.g., sinusoids. RMS is used
  in various fields, including electrical engineering.  It can be
  calculated for a series of discrete values or for a continuously
  varying function. Its name comes from its definition as the square
  root of the mean of the squares of the values.

  *)
  let squares = List.fold_left 
    (fun acc xi -> acc @ [xi *. xi]) [] t in
  sqrt (arithmetic_mean squares)

The generalized mean or "power mean" includes all of the means we have considered to date as special cases.

let power_mean p t =
  (*http://en.wikipedia.org/wiki/Generalized_mean

    In mathematics, generalized means are a family of functions for
    aggregating sets of numbers, that include as special cases the
    arithmetic, geometric, and harmonic means

  *)
  let powers = List.fold_left 
    (fun acc xi -> acc @ [( ** ) xi p]) [] t in
    (arithmetic_mean powers)**(1.0/.p)
Note: Recovering the geometric mean from this definition requires looking at the limit of the expression as p -> 0 and application of L'Hopital's rule (see http://en.wikipedia.org/wiki/Generalized_mean for the details).

Friday, March 21, 2014

Powers

Powers

bn is familiar when n is a whole integer (positive or negative). Here's an algorithm to compute it.
let rec pow (i, u) = 
  if i < 0 then 1.0 / (pow (-i, u))
  else if i = 0 then 1.0 else u * pow ((i-1), u)
It is also familiar to us when c/d = n that is, for n a rational number (then the solution is found by taking the dth root of bc). Less so though I think when x = n is some arbitrary real. That is, generally speaking how is bx computed?

Logarithms and exponentiation. Since b = e log b then bx = (elog b)x = e x (log b).
let powf (b, x) = exp (x * log b)
This only works when b is real and positive. The case negative b leads into the theory of complex powers.

“Sometimes,' said Pooh, 'the smallest things take up the most room in your heart.”

Sunday, November 3, 2013

Matrix multiplication (functional approach in Python)

Matrix multiplication

The last entry looked at Gaussian elimination modeling matrices and vectors as tuples and using some idioms of functional programming. The linear algebra theme is continued here with programs to compute matrix/vector and matrix/matrix multiplications.

First some utilities to warm up on. Functions to project the rows or columns of a matrix.
def proj_r (A, i): return A[i]
def proj_c (A, j): return tuple (map (lambda t: (t[j],), A))
The function that scales a row vector by a constant is here again along with its symmetric partner for scaling a column vector.
def scale_r (u, lam): 
  return tuple (map (lambda ui: lam*ui, u))

def scale_c (c, k):
  def f (e):
    (xi,) = e
    return k*xi,
  return tuple (map (f, c))
For example, if c=((1,),(2,)), then scale_c (c, 2) is ((2,), (4,)).

The utility flatten reinterprets a sequence of column vectors as a matrix of row vectors.
def flatten (t):
  def row(i):
    def f (acc, y): 
        return acc + y[i]
    return functools.reduce (f, t, ())
  return tuple (map (row, range (0, len (t[0]))))
Here is an example. If t1=(1,),(2,) and t2=(3,),(4,) then s=(t1, t2) is a column vector sequence and flatten (s) is the matrix ((1, 3), (2, 4)).

Matrix vector multiplication comes next.
def mat_vec_mul (A, x):
  def f (i):
      (xi,) = x[i]
      return scale_c (proj_c (A, i), xi)
  M = tuple (map (f, range (0, len (x))))
  return tuple (map (lambda r : (sum (r), ), flatten (M)))
This computation utilizes the interpretation of the product as the sum of the column vectors of A weighted by the components of x.

Matrix multiplication follows immediately by seeing it as the flattening of a sequence of matrix/vector multiplications (that is of A on the individual columns of B).
def mat_mat_mul (A, B):
  f = lambda j : mat_vec_mul (A, proj_c (B, j))
  return flatten (tuple (map (f, range (0, len (B[0])))))
The transcript below shows the function at work in the Python top-level.
>>> A = ((2, 3), (4, 0) )
>>> B = ((1,  2, 0), (5, -1, 0) )
>>> print (str (mat_mat_mul (A, B)))
((17, 1, 0), (4, 8, 0))