Wednesday, July 22, 2015

Interpreter of arithmetic expressions (Felix)

This is a tutorial implementing an interpreter of arithmetic expressions. No code generation tools are employed, the program is self-contained. Lexical analysis is performed using parser combinators, evaluation by the method of environments. An interactive program for testing the interpreter is also provided. Read more...

Thursday, June 25, 2015

Cumulative moving average

We have $nA_{N} = s_{1} + \cdots + s_{n}$ and $(n + 1)A_{n + 1} = s_{1} + \cdots + s_{n + 1}$ so $s_{n + 1} = (n + 1)A_{n + 1} - nA_{n}$ from which we observe $A_{n + 1} = \frac{s_{n + 1} + nA_{n}}{n + 1} = A_{n} + \frac{s_{n + 1} + nA_{n}}{n + 1} - A_{n} = A_{n} + \frac{s_{n + 1} - A_{n}}{n + 1}$.


let cumulative_moving_average (l : float list) : float =
  let f (s, i) x =
    let k = i + 1 in
    let s = (s +. (x -. s)/. (float_of_int k)) in
    (s, k) in
  fst @@ List.fold_left f (0., 0) l


fun cumulative_moving_average (l : list[double]) : double = {
  var f = fun (s : double, var i : int) (x : double) : double * int = {
    return (s + ((x - s) / (double)i), i);
  return (fold_left f (0.0, 0) l) .0 ;


#include <numeric>
#include <utility>

namespace detail {
    typedef std::pair<double, int> acc_t;
    inline acc_t cumulative_moving_average (acc_t const& p, double x) {
        double s = p.first;
        int k = p.second + 1;
        return std::make_pair (s + (x - s) / k, k); 
}//namespace detail

template <class ItT>
double cumulative_moving_average (ItT begin, ItT end) {
  return (
    std::accumulate (
      begin, end, std::make_pair (0.0, 0), 


template <class T>
auto cumulative_moving_average(T c) {
    return accumulate (
        std::begin (c), std::end (c), std::make_pair (0.0, 0), 
        [](auto p, auto x) {
          ++std::get<1> (p);
          return std::make_pair (
            std::get<0> (p) + (x - std::get<0> (p)) / 
                         std::get<1> (p), std::get<1> (p));
(credit : Juan Alday)

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

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



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


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
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) }; 
    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) };
    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=
  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$.

Sunday, May 24, 2015

Church Numerals

This is just a little fun. Jason Hickey in "Introduction to Objective Caml" poses some little end of chapter problems to define arithmetic operations for a type of unary (base-1) natural numbers. The type is

type num = Z | S of num
where Z represents the number zero and if i is a unary number, then S i is i + 1.

This formulation of Church numerals using a recursive type and pattern matching means in truth, the problems can be solved in less than 5 minutes or so. Of course, the real Church numerals are numbers encoded in functions

  • $c_{0} = \lambda s.\lambda z.\;z$
  • $c_{1} = \lambda s.\lambda z.\;s\;z$
  • $c_{2} = \lambda s.\lambda z.\;s\;(s\;z)$
  • $c_{3} = \lambda s.\lambda z.\;s\;(s\;(s\;z))$
  • $\cdots$
and their represenation including arithmetic operations can be formulated in OCaml too (and it's a good exercise but harder than what we are going to do here -- if you'd like to see more about that, have a look at this Cornell lecture).

Alright, without further ado, here we go then.

type num = Z | S of num

let scc (x : num) : num = S x
let prd : num -> num = function | S n -> n | _ -> Z

let rec add (x : num) (y : num) : num =
  match (x, y) with 
  | (Z, _) -> y
  | (_, Z) -> x
  | (S m, n) -> scc (add m n)

let rec sub (x : num) (y : num) : num =
  match (x, y) with
  | (Z, _) -> Z
  | (n, Z) -> n
  | (S m, n) -> sub m (prd n)

let rec mul (x : num) (y : num) : num = 
  match (x, y) with
  | (Z, _) -> Z
  | (_, Z) -> Z
  | (S Z, x) -> x
  | (x, S Z) -> x
  | (S m, n) -> add (mul m n) n

let rec to_int : num -> int = function | Z -> 0 | S n -> 1 + to_int n
let rec from_int (x : int)  = if x = 0 then Z else scc (from_int (x - 1))
For example, in the top-level we can write things like,
# to_int (mul (sub (from_int 23) (from_int 11)) (from_int 2));;
- : int = 24

The main thing I find fun about this little program though is how obvious its mapping to C++. Of course you need a discriminated union type my default choice being boost::variant<> (by the way, standardization of a variant type for C++ is very much under active discussion and development, see N4450 for example from April this year - it seems that support for building recursive types might not be explicitly provided though... That would be a shame in my opinion and if that's the case, I beg the relevant parties to reconsider!).

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

#include <stdexcept>
#include <iostream>

struct Z;
struct S;

typedef boost::variant<Z, boost::recursive_wrapper<S>> num;

struct Z {};
struct S { num i; };

int to_int (num const& i);

struct to_int_visitor 
  : boost::static_visitor<int> {
  int operator ()(Z const& n) const { return 0; }
  int operator ()(S const& n) const { return 1 + to_int (n.i); }

int to_int (num const& i) {
  return boost::apply_visitor (to_int_visitor (), i);

num from_int (int i) {
  if (i == 0){
    return Z {};
    return S {from_int (i - 1)};

num add (num l, num r);

struct add_visitor : boost::static_visitor<num> {
  num operator () (Z, S s) const { return s; }
  num operator () (S s, Z) const { return s; }
  num operator () (Z, Z) const { return Z {}; }
  num operator () (S s, S t) const { return S { add (s.i, t) }; }

num add (num l, num r) {
  return boost::apply_visitor (add_visitor (), l, r);

num succ (num x) { return S{x}; }

struct prd_visitor : boost::static_visitor<num>{
  num operator () (Z z) const { return z; }
  num operator () (S s) const { return s.i; }

num prd (num x) {
  return boost::apply_visitor(prd_visitor (), x);

num sub (num x, num y);

struct sub_visitor : boost::static_visitor<num> {
  num operator () (Z, Z) const { return Z {}; }
  num operator () (Z, S) const { return Z {}; }
  num operator () (S m, Z) const { return m; }
  num operator () (S m, S n) const { return sub (m.i, prd (n)); }

num sub (num x, num y) {
  return boost::apply_visitor (sub_visitor (), x, y);


int main () {

  num zero = Z {};
  num one = succ (zero);
  num two = succ (succ (zero));
  num three = succ (succ (succ (zero)));

  std::cout << to_int (add (two, three)) << std::endl;
  std::cout << to_int (sub (from_int (23), from_int (12))) << std::endl;

  return 0;
I didn't get around to implementing mul in the above. Consider it an "exercise for the reader"!

Thursday, March 19, 2015

Y Combinator

There is this blog post by Caltech computer scientist, Mike Vanier. The code in Mike's article uses the Scheme programming language. This one uses Python.

A $Y$ combinator is a higher-order function which, given argument $f$ (say,) satisfies the equation $Y f = f\;(Y f)$. They are said to compute a "fixed point" $f^{'}$ of their argument $f$ since $f^{'} = Y\;f = f\;(Y\;f) = f \;f^{'}$.

A $Y$ combinator takes a function that isn't recursive and returns a version of the function that is. That is, $Y$ is a function that takes $f$ and returns $f (f (f (\cdots)))$.

The existence of $Y$ combinators is amazing in that it tells us that it is possible to write recursive functions even in a programming language that says nothing about recursion!

The goal here is to derive a $Y$.

Start with the classic recursive definition of the factorial function.

  def fact (n) :
    if n == 0 : 
      return 1 
      return n * fact (n - 1)

We are trying to eliminate explicit recursion. To that end, factor out the recursive call and make it an application of a function argument.

def part_fact (this, n):
  if n == 0 :
    return 1
     return n * this (this, (n - 1))

fact = functools.partial(part_fact, part_fact)
That's sufficient to get rid of the explicit recursion but we mean to push on in search of a "library" solution to this problem, that is, some general result that can be re-used.

Next let's get this down to a function in one argument as is this way in the $\lambda$ calculus.

  def part_fact (this):
    return lambda n : 1 if n == 0 else n * (this (this)) (n - 1)

  fact = part_fact (part_fact)

We'd recover something tantalizingly close to the original factorial function if we factored out this (this) into a function of one argument.

def part_fact (this):
  f = this (this)
  return lambda n : 1 if n == 0 else n * f (n - 1)
fact = part_fact(part_fact)

This would be fine in a lazy language but Python is a strict language and exhibits infinite recursion because to evaluate part_fact (part_fact) requires evaluating f = part_fact (part_fact) and so on. The solution is to delay the evaluation until it's needed.

def part_fact (this):
  f = lambda y : (this (this)) y
  return lambda n : 1 if n == 0 else n * f (n - 1)
fact = part_fact(part_fact)

Refactor this into two parts.

def almost_fact (f):
  return lambda n : 1 if n == 0 else f (n - 1)

def part_fact (this):
  f = lambda y : (this (this))(y)
  return almost_fact (f)

fact = part_fact(part_fact)

Rephrase part_frac as a lambda and change the argument name to x.

def almost_fact (f):
  return lambda n : 1 if n == 0 else f (n - 1)

part_fract = lambda x : almost_fact (lambda y : (x (x))(y))

fact = part_fact (part_fact)

Eliminate 'part_fact'.

def almost_fact (f):
  return lambda n : 1 if n == 0 else f (n - 1)

fact = (lambda x : almost_fact (lambda y : (x (x))(y))) \
          (lambda x : almost_fact (lambda y : (x (x))(y)))

That's it! There's the $Y$ combinator. Generalize.

def Y (f):
 return (lambda x : f (lambda y : (x (x))(y))) \
           (lambda x : f (lambda y : (x (x))(y)))

def almost_fact (f):
  return lambda n : 1 if n == 0 else f (n - 1)

fact = Y (almost_fact)
That is, $Y = \lambda f.(\lambda x. f\;(\lambda y. (x\;x)\;y))\;(\lambda x. f\;(\lambda y. (x\;x)\;y))$. This $Y$ combinator is known s Curry's paradoxical combinator (in its "applicative order" form to account for the fact that Python is strict).

Try this on another function. Fibonacci numbers say.

def almost_fib (f) :
  return lambda n : 1 if n <= 2 else f (n - 1) + f (n - 2)

fib = Y (almost_fib)

print (str (fib (6))+"\n") #Prints '8'

Saturday, March 14, 2015

Labeled and optional arguments

I don't know why this is, but of all the features of the OCaml language, somehow remembering how to make use of labeled or optional arguments is difficult for me. I find when I need to use them, I have to go back and "relearn" the syntax every time. For that reason, I've written these summary notes for myself (based mostly on Jacques Garrigue's "Labels and Variants" chapter of the OCaml reference) and I hope they might be of help you too!


I've seen people write this function from time to time.
let flip (f : 'a -> 'b -> 'c) = fun y x -> f x y

flip applied to a function f in two arguments, produces a new function that when presented reverses their order in the application of f. For example, if sub is defined by let sub x y = x - y, then sub 3 2 will give the result 1 whereas (flip sub) 3 2 will produce the value -1.

The reason you'd find yourself wanting a function like flip is because you want to partially apply f but the argument you want to bind is inconveniently in the 'wrong' position. For example, the function that subtracts 2 from it's argument can be written (flip sub 2) (whereas sub 2 is the function that subtracts its argument from 2).

Labeled arguments do away with the need for functions like flip. You see, arguments in function types can be annotated.

val sub : (x:int) -> (y:int) -> int
In the above, x and y are labels. When you define a function with labeled arguments, you prefix the labels with a tilde ~.
let sub : (x:int -> y:int -> int) = fun ~x ~y -> x - y

Labels and variables are not the same things. A label is a name associated with a variable e.g. let u = 3 and v = 2 in sub ~x:u ~y:v (u, v are variables x, y labels) and this also works when the variables are replaced with literals as in sub ~x:3 ~y:2.

The expression ~x on it's own is shorthand for ~x:x which means that you can conveniently write let x = 3 and y = 2 in sub ~x ~y for example. This is "punning" - meaning if the variable name and the label are the same, the variable name is permitted to be elided in the application.

Labels enable making function application commutative (changing the order of the arguments does not change the result); one can write sub ~y ~x and that will yield the same result as sub ~x ~y. Where this is useful of course is partial application of a function on any any argument. For example, we can create functions from sub by binding either of x or y such as sub ~x:3 or sub ~y:2 without having to resort to ad-hoc trickery like flip sub.

Some technical details:

  • If more than one argument of a function has the same label (or no label), position prevails and they will not commute among themselves;
  • If an application is total, labels can be omitted (e.g. like in sub 3 2). Be wary though - when a function is polymorphic in its return type, it will never be considered as totally applied. Consider for example , let app ~f ~x = f x. Then app has type f:('a -> 'b') -> x:'a -> 'b. Now we attempt total application but omit the required labels app (fun x -> x) 1 and this types to f:('a -> ('b -> 'b) -> int -> 'c) -> x:'a -> 'c. You can see that the arguments have been "stacked up" in anticipation of presentation of f and x and the result although it types, is not at all what we were going for!
  • When a function is passed to a higher order function, labels must match in both types, labels may not be added or removed.

Optional arguments

In a function type involving an optional argument, the annotation ? appears.
val range : ?step:int -> int -> int -> int list
The ? appears again in the value syntax of a function definition involving an optional argument. Additionally, at that time, optional parameters can be given default values. When calling a function with a value for an optional parameter, well, you use the fact it is a labeled argument and provide a ~ as normal.
let rec range ?(step:int = 1) (m : int) (n : int) : int list = 
    if m >= n then [] 
    else m :: (range ~step (m + step) n)

Now, this bit is important. A function can not take optional arguments alone. There must be some non-optional ones there too and it is when a non-optional parameter that comes after an optional one in the function type is encountered is it determined if the optional parameter has been omitted or not. For example, if we reorder the argument list as in the below, we find we can't 'erase' the optional argument anymore.

let rec range (m : int) (n : int) ?(step : int = 1): int list = ...
Warning 16: this optional argument cannot be erased.
That said, optional parameters may commute with non-optional or unlabeled ones, as long as they are applied simultaneously. That is, going back to this definition,
let rec range ?(step:int = 1) (m : int) (n : int) : int list = 
    if m >= n then [] 
    else m :: (range ~step (m + step) n)
(range 0 ~step:2) 100 is valid whereas, (range 0) ~step:2 100 is not.

Optional arguments are implemented as option types. In the event no default is given, you can treat them as exactly that by matching on None or Some x to switch behaviors. Here's a schematic of what that looks like.

let f ?s y =
    match s with
      | None -> ...
      | Some x -> ...  

When you are just "passing through" an optional argument from one function call to another, you can prefix the applied argument with ? as in the following.

 let g ?s x = ...
 let h ?s x = g ?s x  (*[s] is passed through to [g]*)

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

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

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

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

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
      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 = 

(*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) () ->
       type t = s 
       let 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]