Saturday, November 28, 2015

C++ : Sums with constructors

C++ : Sums with constructors

I've been working recently on a type to model "sums with constructors" in C++ (ala OCaml). The implementation technique is "novel" in that it makes use of C++11's "unrestricted unions" feature. I learned it from the FTL library where the idea is credited to Björn Aili. FTL also shows how to provide a NEAT (for C++) syntax for pattern matching but, unless I just didn't get it, the FTL version doesn't admit recursive types "out-of-the-box". So, I extended Björn's work to admit recursive types by applying the recursive wrapper idea from Boost.Variant (Eric Friedman, Itay Maman). The resulting library, I call the "pretty good sum" library. It's C++14 but can be back-ported to C++11 (update : that's done and a lot of bug-fixes to). The code is online here if you want to play with it in your own programs.

There are a number of usage examples provided in the library tests/documentation. I'll provide a small one here - the ubiquitous option<> type (c.f. Boost.Optional and OCaml's builtin type α option).

In OCaml, the type definition is given by

type α option = Some of α | None
which is not recursive (see the other examples on github for that e.g. functional lists, abstract syntax trees) but I hope this example is still interesting in that it explores the type's monadic nature to implement so called "safe-arithmetic", that is, integer arithmetic that guards against overflow and division by zero (source : "Ensure that operations on signed integers do not result in overflow"). See this post for more on monads in C++.

The code in the example is fairly extensively commented so I hope you will excuse me this time if I don't provide my usual narrative (I've presented this program before in a Felix tutorial - there's a narrative there note to self : and some typos that I mean to get back to and fix).

Without further ado... A type for optional values in C++ using the "pretty good sum" type!

#include <pgs/pgs.hpp>

#include <gtest/gtest.h>

#include <iostream>
#include <cstdlib>
#include <climits>
#include <functional>

//type 'a t = Some of 'a | None

namespace {

using namespace pgs;

template <class T>
struct some_t { //Case 1
  T data;  
  template <class U>
  explicit some_t (U&& data) : data { std::forward<U> (data) }
  {}
};

struct none_t //Case 2
{};

//Options are a type that can either hold a value of type `none_t`
//(undefined) or `some_t<T>`
template<class T>
using option = sum_type<some_t<T>, none_t>;

//is_none : `true` if a `some_t<>`, `false` otherwise
template<class T>
bool is_none (option<T> const& o) {
  return o.template is<none_t> ();
}

//A trait that can "get at" the type `T` contained by an option
template <class>
struct option_value_type;
template <class T>
struct option_value_type<option<T>> { typedef T type; };
template <class T>
using option_value_type_t = typename option_value_type<T>::type;

//Factory function for case `none_t`
template <class T>
option<T> none () {
  return option<T>{constructor<none_t>{}};
}

//Factory function for case `some_t<>`
template <class T>
option<decay_t<T>> some (T&& val) {
  using t = decay_t<T>;
  return option<t>{constructor<some_t<t>>{}, std::forward<T> (val)};
}

//is_some : `false` if a `none_t`, `true` otherwise
template<class T>
inline bool is_some (option<T> const& o) {
  return o.template is<some_t<T>>();
}

//Attempt to get a `const` reference to the value contained by an
//option
template <class T>
T const& get (option<T> const & u) {
  return u.template match<T const&> (
   [](some_t<T> const& o) -> T const& { return o.data; },
   [](none_t const&) -> T const& { throw std::runtime_error {"get"}; }
  );
}

//Attempt to get a non-`const` reference to the value contained by an
//option
template <class T>
T& get (option<T>& u) {
  return u.template match<T&> (
   [](some_t<T>& o) -> T& { return o.data; },
   [](none_t&) -> T& { throw std::runtime_error {"get"}; }
   );
}

//`default x (Some v)` returns `v` and `default x None` returns `x`
template <class T>
T default_ (T x, option<T> const& u) {
  return u.template match<T> (
    [](some_t<T> const& o) -> T { return o.data; },
    [=](none_t const&) -> T { return x; }
  );
}

//`map_default f x (Some v)` returns `f v` and `map_default f x None`
//returns `x`
template<class F, class U, class T>
auto map_default (F f, U const& x, option<T> const& u) -> U {
  return u.template match <U> (
    [=](some_t<T> const& o) -> U { return f (o.data); },
    [=](none_t const&) -> U { return x; }
  );
}

//Option monad 'bind'
template<class T, class F>
auto operator * (option<T> const& o, F k) -> decltype (k (get (o))) {
  using result_t = decltype (k ( get (o)));
  using t = option_value_type_t<result_t>;
  return o.template match<result_t>  (
      [](none_t const&) -> result_t { return none<t>(); }, 
      [=](some_t<T> const& o) -> result_t { return k (o.data); }
  );
}

//Option monad 'unit'
template<class T>
option<decay_t<T>> unit (T&& a) {
  return some (std::forward<T> (a));
}

//map
template <class T, class F>
auto map (F f, option<T> const& m) -> option<decltype (f (get (m)))>{
  using t = decltype (f ( get (m)));
  return m.template match<option<t>> (
      [](none_t const&) -> option<t> { return none<t>(); }, 
      [=](some_t<T> const& o) -> option<t> { return some (f (o.data)); }
  );
}

}//namespace<anonymous>

TEST (pgs, option) {
  ASSERT_EQ (get(some (1)), 1);
  ASSERT_THROW (get (none<int>()), std::runtime_error);
  auto f = [](int i) { //avoid use of lambda in unevaluated context
    return some (i * i);   };
  ASSERT_EQ (get (some (3) * f), 9);
  auto g = [](int x) { return x * x; };
  ASSERT_EQ (get (map (g, some (3))), 9);
  ASSERT_TRUE (is_none (map (g, none<int>())));

  ASSERT_EQ (default_(1, none<int>()), 1);
  ASSERT_EQ (default_(1, some(3)), 3);
  auto h = [](int y) -> float{ return float (y * y); };
  ASSERT_EQ (map_default (h, 0.0, none<int>()), 0.0);
  ASSERT_EQ (map_default (h, 0.0, some (3)), 9.0);
}

namespace {

//safe "arithmetic"

std::function<option<int>(int)> add (int x) {
  return [=](int y) -> option<int> {
    if ((x > 0) && (y > INT_MAX - x) ||
        (x < 0) && (y < INT_MIN - x)) {
        return none<int>(); //overflow
      }
    return some (y + x);
  };
}

std::function<option<int>(int)> sub (int x) {
  return [=](int y) -> option<int> {
    if ((x > 0) && (y < (INT_MIN + x)) ||
        (x < 0) && (y > (INT_MAX + x))) {
      return none<int>(); //overflow
    }
    return some (y - x);
  };
}

std::function<option<int>(int)> mul (int x) {
  return [=](int y) -> option<int> {
    if (y > 0) { //y positive
      if (x > 0) {  //x positive
        if (y > (INT_MAX / x)) {
          return none<int>(); //overflow
        }
      }
      else { //y positive, x nonpositive
        if (x < (INT_MIN / y)) {
          return none<int>(); //overflow
        }
      }
    }
    else { //y is nonpositive
      if (x > 0) { // y is nonpositive, x is positive
        if (y < (INT_MIN / x)) {
          return none<int>();
        }
      }
      else { //y, x nonpositive 
        if ((y != 0) && (x < (INT_MAX / y))) {
          return none<int>(); //overflow
        }
      }
    }

    return some (y * x);
  };
}

std::function<option<int>(int)> div (int x) {
  return [=](int y) {
    if (x == 0) {
      return none<int>();//division by 0
    }

    if (y == INT_MIN && x == -1)
      return none<int>(); //overflow

    return some (y / x);
  };
}

}//namespace<\anonymous>

TEST(pgs, safe_arithmetic) {

  //2 * (INT_MAX/2) + 1 (won't overflow since `INT_MAX` is odd and
  //division will truncate)
  ASSERT_EQ (get (unit (INT_MAX) * div (2) * mul (2) * add (1)), INT_MAX);

  // //2 * (INT_MAX/2 + 1) (overflow)
  ASSERT_TRUE (is_none (unit (INT_MAX) * div (2) * add (1) * mul (2)));

  // //INT_MIN/(-1)
  ASSERT_TRUE (is_none (unit (INT_MIN) * div (-1)));
}

Saturday, October 31, 2015

C++ : Folds over variadic templates

C++ : Folds over variadic templates

Code like the following motivates the need to compute conjunctions (and disjunctions) of predicate packs.

template <class T, class... Ts>
struct recursive_union {

  // ...

  //'U' is not 'T' but 'T' is a recursive wrapper and 'U' is the type
  //contained in 'T'
  template <class U, class... Args,
  std::enable_if_t<
     and_<
      is_recursive_wrapper<T>
    , std::is_same<U, unwrap_recursive_wrapper_t<T>>>::value, int> = 0
  >
  explicit recursive_union (constructor<U>, Args&&... args)
    noexcept (std::is_nothrow_constructible<U, Args...>::value)
  : v (std::forward<Args>(args)...)
  {}

  // ...

};
I was much helped by a suitable implementation of and_<> credited to Jonathan Wakely. A more general approach to is to use fold.
#include <type_traits>

namespace pgs {

template<class F, class Acc, class... Ts>
struct fold_left : Acc {
};

template <class F, class Acc, class T, class... Ts>
struct fold_left<F, Acc, T, Ts...> : 
    fold_left <F, typename F::template apply<Acc, T>::type, Ts...> {
};

//or

struct or_helper {
  template <class Acc, class T>
  struct apply : std::integral_constant<bool, Acc::value || T::value> {
  };
};

template <class... Ts>
struct or_ : fold_left <or_helper, std::false_type, Ts...> {
};

//and

struct and_helper {
  template <class Acc, class T>
  struct apply : std::integral_constant<bool, Acc::value && T::value> {
  };
};

template <class... Ts>
struct and_ : fold_left <and_helper, std::true_type, Ts...> {
};

}//namespace pgs

Friday, October 9, 2015

Expression algebras (Python)

This is a Python version of the C++ program presented in this earlier blog entry on expression algebras.

import functools
_isconst = \
  lambda x :functools.reduce \
    (lambda acc, c : acc and isinstance (c, _const), x, True)

class _float :
  def __neg__ (self) :
    return _const (-self.f) if _isconst ([self]) else _neg (self)
  def __add__ (self, x) : 
      return _const (self.f + x.f) if _isconst ([self, x]) else \
      x if _isconst ([self]) and self.f == 0 else               \
      self if _isconst ([x]) and x.f == 0 else _add (self, x)
  def __sub__ (self, x) : 
      return _const (self.f - x.f) if _isconst ([self, x]) else \
      const (-x.f) if _isconst ([self]) and self.f == 0 else    \
      self if _isconst ([x]) and x.f == 0.0 else _sub (self, x)
  def __mul__ (self, x) : 
      return _const (self.f * x.f) if _isconst ([self, x]) else \
      x if _isconst ([self]) and self.f == 1 else               \
      self if _isconst ([x]) and x.f == 1 else _mul (self, x)
  def __div__ (self, x) : 
      return _const (self.f / x.f) if _isconst ([self, x]) else \
      self if _isconst([x]) and x.f == 1 else _div (self, x)

class _neg (_float):
  def __init__ (self, f) : self.f = f
  def __str__ (self) : return "-" + "(" + str(self.f) + ")"
class _fix (_float):
  def __init__ (self, d, f) : self.d = d; self.f = f
  def __str__ (self) : return "fix(" + str(self.d) +", " + str(self.f) + ")"
class _add (_float) :
  def __init__ (self, lhs, rhs) : self.lhs = lhs; self.rhs = rhs
  def __str__ (self) : return str(self.lhs)+ " + " + str(self.rhs)
class _sub (_float):
  def __init__ (self, lhs, rhs) : self.lhs = lhs; self.rhs = rhs
  def __str__ (self) : return str(self.lhs)+ " - " + str(self.rhs)
class _mul (_float):
  def __init__ (self, lhs, rhs) : self.lhs = lhs; self.rhs = rhs
  def __str__ (self) : return str (self.lhs)+ " * " + str (self.rhs)
class _div (_float):
  def __init__ (self, lhs, rhs) : self.lhs = lhs; self.rhs = rhs
  def __str__ (self) : return str (self.lhs)+ " / " + str (self.rhs)
class _const (_float):
  def __init__ (self, f) : self.f = f;
  def __str__ (self) : return str (self.f)
class _obs (_float):
  def __init__ (self, tag) : self.tag = tag
  def __str__ (self) : return "observation \"" + str(self.tag) + "\""
class _max (_float):
  def __init__ (self, lhs, rhs) : 
      self.lhs = lhs; self.rhs = rhs
  def __str__ (self) : 
    return "max(" + str (self.lhs) + ", " + str (self.rhs) + ")"
class _min (_float):
  def __init__ (self, lhs, rhs) : 
      self.lhs = lhs; self.rhs = rhs
  def __str__ (self): 
    return "min(" + str (self.lhs) + ", " + str (self.rhs) + ")"

def visit (f, acc, xpr):
  if isinstance (xpr, _const) : return f._const (acc, xpr)
  if isinstance (xpr, _neg) : return f._neg (acc, xpr)
  if isinstance (xpr, _fix) : return f._fix (acc, xpr)
  if isinstance (xpr, _obs) : return f._obs (acc, xpr)
  if isinstance (xpr, _add) : return f._add (acc, xpr)
  if isinstance (xpr, _sub) : return f._sub (acc, xpr)
  if isinstance (xpr, _mul) : return f._mul (acc, xpr)
  if isinstance (xpr, _div) : return f._div (acc, xpr)
  if isinstance (xpr, _max) : return f._max (acc, xpr)
  if isinstance (xpr, _min) : return f._min (acc, xpr)

  raise RuntimeError ("Expression match failure")

const = lambda c : _const (c)
observation = lambda s : _obs (s)
max_ = lambda a, b : _max (a, b)
min_ = lambda a, b : _min (a, b)

def fix (d, x):

  class __fix_visitor:
    def __init__ (self, d) : 
      self.d = d
    def _const (self, _, xpr) : 
      return xpr
    def _obs (self, _, xpr) : 
      return _fix (self.d, xpr)
    def _fix (self, _, xpr) : return xpr
    def _neg (self, _, xpr) : 
      return _neg (visit (self, _, xpr.f))
    def _add (self, _, xpr) : 
      return _add (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _sub (self, _, xpr) : 
      return _sub (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _mul (self, _, xpr) : 
      return _mul (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _div (self, _, xpr) : 
      return _div (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _max (self, _, xpr) : 
      return _max (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _min (self, _, xpr) : 
      return _min (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))

    return visit (__fix_visitor (d), None, x)

def simplify (fs, x):

  class _apply_fixings_visitor :
    def __init__(self, fs) : self.fs = fs
    def _const (self, _, xpr) : return xpr
    def _obs (self, _, xpr) : return xpr
    def _fix (self, _, xpr) : 
      fs = [f for f in self.fs if f[0] == xpr.f.tag and f[1] == xpr.d]
      return xpr if len (fs) == 0 else _const (fs[0][2])
    def _neg (self, _, xpr) : 
      return _neg (visit (self, _, xpr.f))
    def _add (self, _, xpr) : 
      return _add (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _sub (self, _, xpr) : 
      return _sub (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _mul (self, _, xpr) : 
      return _mul (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _div (self, _, xpr) : 
      return _div (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _max (self, _, xpr) : 
      return _max (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))
    def _min (self, _, xpr) : 
      return _min (visit (self, _, xpr.lhs), visit (self, _, xpr.rhs))

  class _simplify_visitor:
    def _const (self, _, xpr) : 
      return xpr
    def _fix (self, _, xpr) : 
      return xpr
    def _obs (self, _, xpr) : 
      return xpr
    def _neg (self, _, xpr) : 
      f = visit (self, _, xpr.f)
      return xpr if not _isconst ([f]) else -f
    def _add (self, _, xpr) : 
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (l.f + r.f)
    def _sub (self, _, xpr) :
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (l.f - r.f)
    def _mul (self, _, xpr) :
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (l.f * r.f)
    def _div (self, _, xpr) :
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (l.f / r.f)
    def _max (self, _, xpr) :
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (max (l.f, r.f))
    def _min (self, _, xpr) :
      l = visit (self, _, xpr.lhs); r = visit (self, _, xpr.rhs)
      return xpr if not _isconst([l, r]) else const (min (l.f, r.f))

  return visit ( \
    _simplify_visitor (), None, visit (_apply_fixings_visitor (fs), None, x))

Sunday, October 4, 2015

List comprehensions in C++ via the list monad

Monads

As explained in Monads for functional programming by Philip Wadler, a monad is a triple $(t, unit, *)$. $t$ is a parametric type, $unit$ and $*$ are operations:

  val unit : α -> α t
  val ( * ) : α t -> (α -> β t) -> β t
We can read expressions like

$m * \lambda\;a.n$

as, "perform computation $m$, bind $a$ to the resulting value, and then perform computation $n$". Referring to the signatures of $*$ and $unit$, in terms of types we see $m$ has the type α t, $\lambda\;a.n$ has type α -> β t and the whole expression has type β t.

In order for $(t, unit, *)$ to be a monad the operations $unit$ and $*$ need satisfy three laws :

  • Left unit. Compute the value $a$, bind $b$ to the result, and compute $n$. The result is the same as $n$ with value $a$ substituted for variable $b$.

    $unit\;a * \lambda\;b.n = n[a/b]$.

  • Right unit. Compute $m$, bind the result to $a$, and return $a$. The result is the same as $m$.

    $m * \lambda\;a.unit\;a = m$.

  • Associative. Compute $m$, bind the result to $a$, compute $n$, bind the result to $b$, compute $o$. The order of parentheses doesn't matter.

    $m * (\lambda\;a.n * \lambda\;b.o) = (m * \lambda\;a.n) * \lambda\;b.o$.

The list monad

Lists can be viewed as monads.That is, there exist operations $unit$ and $*$ that we may define for lists such that the three monad laws from the preceding section hold.

#include <list>
#include <iterator>
#include <type_traits>
#include <algorithm>
#include <iostream>

/*
  The list monad
*/

//The unit list containing 'a'
/*
  let unit : 'a -> 'a t = fun a -> [a]
*/
template <class A> 
std::list<A> unit (A const& a) { return std::list<A> (1u, a); }

//The 'bind' operator
/*
  let rec ( * ) : 'a t -> ('a -> 'b t) -> 'b t =
    fun l -> fun k ->
      match l with | [] -> [] | (h :: tl) -> k h @ tl * k
*/
template <class A, class F>
typename std::result_of<F(A)>::type 
operator * (std::list<A> a, F k) {
  typedef typename std::result_of<F(A)>::type result_t;

  if (a.empty ())
    return result_t ();

  result_t res = k (a.front ());
  a.pop_front ();
  res.splice (res.end (), a * k);

  return res;
}
The invocation $unit\;a$ forms the unit list containing $a$. The expression, $m * k$ applies $k$ to each element of the list $m$ and appends together the resulting lists.

There are well known derived forms. For example, $join\;z$ is the expression $z * \lambda\;m. m$. In the list monad, it results in a function that concatenates a list of lists.

//'join' concatenates a list of lists
/*
    let join : 'a t t z = z * fun m -> m
*/
template <class A>
std::list <A> join (std::list<std::list<A>> const& z) {
  return z * [](auto m) { return m; };
}
The function $map$ is defined by the expression $map\;f\;m = m * \lambda\;a.unit\;(f\;a)$.
//'map' is the equivalent of 'std::transform'
/*
    let map : ('a -> b') -> 'a t -> 'b t =
      fun f -> fun m -> m * fun a -> unit (f a)
*/
template <class A, class F>
std::list<A> map (F f, std::list<A> const& m) {
  return m * [=](auto a) { return unit (f (a)); };
}

List comprehensions

List comprehensions are neatly expressed as monad operations. Here are some examples.
int main () {

  //l = [1, 2, 3]
  std::list<int> l = {1, 2, 3};
  
  //m = [1, 4, 9]
  auto m = l * [](int x) { return unit (float (x * x)); };

  //n = l x m = [(1, 1), (1, 4), (1, 9), (2, 1), (2, 4), (2, 9), ...]
  auto n = l * ([&m](int x){ return m * ([=](float y){ return unit (std::make_pair (x, y)); });});

  return 0;
}

Sunday, September 20, 2015

Expression algebras

Expression algebras

Some phenomena can be conveniently described as expressions in observations of one or more processes. To do this requires developing an algebra of real numbers with "placeholders" where a placeholder represents a process.

For example, consider writing software for some kind of "smart watch" fitness device that includes a heart rate monitor. We might define heart rate performance over an interval of time by something like

$r/r_{t} - 1.0$ .

If $r$ "fixes" at some time $t +\Delta t$ to a value greater than $r_{t}$, performance will be positive. If it fixes to the same value as $r_{t}$ then performance will be $0$ and if $r_{t + \Delta t}$ is observed to be less than $r_{t}$, degradation has occurred.

In our smart watch application we may in some contexts be interested only in the positive part of heart-rate performance, i.e.

$(r/r_{t} - 1.0)^{+} = max (0, r/r_{t} - 1.0)$ .

We can create further derived expressions that scale and "cap" performance

$100 * min (0.2, max (0, r/r_{t} - 1.0))$ .

Of course, for the above formula to result in physicals amounts we must specify concretely when the heart-rate is to be observed e.g.

$fix (t + \Delta t, 100 * min (0.2, max (0, r/r_{t} - 1.0))) = $

$100 * min (0.2, max (0, r_{t + \Delta t}/r_{t} - 1.0))$ .

Let's look at one way to express all this in C++. The core idea is to model each primitive expression case as a C++ type. The recursive expression type itself is a discriminated union.

#include <string>
#include <iostream>
#include <functional>
#include <algorithm>

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

struct _neg;   
struct _fix;   
struct _add;   
struct _sub; 
struct _mul;   
struct _div;   
struct _max;   
struct _min;
struct _const; 
struct _obs; 

using expression = 
  boost::variant<
      _const
    , _obs
    , boost::recursive_wrapper<_neg>
    , boost::recursive_wrapper<_fix>
    , boost::recursive_wrapper<_add>
    , boost::recursive_wrapper<_sub>
    , boost::recursive_wrapper<_mul>
    , boost::recursive_wrapper<_div>
    , boost::recursive_wrapper<_min>
    , boost::recursive_wrapper<_max>  
>;

Now that we've got the definition of expression out of the way, we can "go back" and fill in the skipped details of the expression classes.

struct _const { double f; };
struct _obs { std::string tag; };
struct _neg { expression f; };
struct _fix { std::string d; expression f; };
struct _add { expression lhs; expression rhs; };
struct _sub { expression lhs; expression rhs; };
struct _mul { expression lhs; expression rhs; };
struct _div { expression lhs; expression rhs; };
struct _max { expression lhs; expression rhs; };
struct _min { expression lhs; expression rhs; };
I'm using a simple std::string to represent a fixing time $t$. Of course in the real world that would probably be some kind of datetime.

The stream insertion operators aren't pretty but they are regular.

std::ostream& operator << (std::ostream& os, _neg const& e) {
  return os << "-(" << e.f << ")";
}
std::ostream& operator << (std::ostream& os, _fix const& e) {
  return os << "fix(" << e.d << ", " << e.f << ")";
}
std::ostream& operator << (std::ostream& os, _add const& e) {
  return os << e.lhs << " + " << e.rhs;
}
std::ostream& operator << (std::ostream& os, _sub const& e) {
  return os << e.lhs << " - " << e.rhs;
}
std::ostream& operator << (std::ostream& os, _mul const& e) {
  return os << e.lhs << " * " << e.rhs;
}
std::ostream& operator << (std::ostream& os, _div const& e) {
  return os << e.lhs << " / " << e.rhs;
}
std::ostream& operator << (std::ostream& os, _max const& e) {
  return os << "max (" << e.lhs << ", " << e.rhs << ")";
}
std::ostream& operator << (std::ostream& os, _min const& e) { 
  return os << "min (" << e.lhs << ", " << e.rhs << ")";
}
std::ostream& operator << (std::ostream& os, _obs const& e) {
  return os << "observation(\"" << e.tag << "\")";
}
std::ostream& operator << (std::ostream& os, _const const& e) {
  return os << e.f;
}

Now things start to get interesting. Definitions of the elementary arithmetic operators. This bit is worth reading closely to make sure you understand it before moving on.

namespace detail {

  struct neg_visitor :
    boost::static_visitor {
    expression operator()(_const lhs) const {
      return _const {-lhs.f};
    }

    template  
    expression operator ()(T const& x) const { 
      _neg res; res.f = x; return res; 
    }
  };

  template <class Op>
  struct binop_visitor 
    : boost::static_visitor<expression>{
    typedef std::function<double(double, double)> op_t;
    op_t op;
    binop_visitor (op_t const& op) : op(op)
    {}
    expression operator()(_const lhs, _const rhs) const { 
      return _const {op (lhs.f, rhs.f)}; 
    }
    template <class T, class U> 
    expression operator ()(T const& x, U const& y) const { 
      Op res; res.lhs = x; res.rhs = y; return res;
    }
  };

}//namespace detail

expression operator - (expression const& x) {
  return boost::apply_visitor (detail::neg_visitor (), x);
}

expression operator + (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_add> (
       [](double x, double y) ->double { return x + y; })
     , lhs, rhs);
}
expression operator - (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_sub> (
       [](double x, double y) ->double { return x - y; })
     , lhs, rhs);
}
expression operator * (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_mul> (
       [](double x, double y) -> double { return x * y; })
     , lhs, rhs);
}
expression operator / (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_div> (
       [](double x, double y) -> double { return x / y; })
     , lhs, rhs);
}
expression max (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_max> (
        [](double x, double y) -> double { return (std::max) (x, y); })
     , lhs, rhs);
}
expression min (expression const& lhs, expression const& rhs) {
  return boost::apply_visitor (detail::binop_visitor<_min> (
        [](double x, double y) -> double { return (std::min) (x, y); })
     , lhs, rhs);

}

Here are the basic expression factory functions modulo fix.

expression const_ (double f) { return _const{f}; }
expression observation (std::string const& tag) { _obs m; m.tag = tag; return m; }

Here then is fix.

expression fix (std::string const& d, expression const& x);

namespace detail {

  struct fix_visitor :
    boost::static_visitor<expression> {
    std::string d;
    fix_visitor (std::string const& d) : d(d)
    {}
    expression operator ()(_const const& x) const
    { return x; }
    expression operator ()(_obs const& x) const
    {  _fix res;
     res.d = d; res.f = x; return res; }
    expression operator ()(_neg const& x) const
    { return x; }
    expression operator ()(_fix const& x) const
    { return x; }
    expression operator ()(_add const& x) const
    { _add res; res.lhs = fix (d, x.lhs);
      res.rhs = fix (d, x.rhs); return res; }
    expression operator ()(_sub const& x) const
    { _sub res; res.lhs = fix (d, x.lhs);
       res.rhs = fix (d, x.rhs); return res; }
    expression operator ()(_mul const& x) const
    { _mul res; res.lhs = fix (d, x.lhs);
      res.rhs = fix (d, x.rhs); return res; }
    expression operator ()(_div const& x) const
    { _div res; res.lhs = fix (d, x.lhs);
      res.rhs = fix (d, x.rhs); return res; }
    expression operator ()(_max const& x) const
    { _max res; res.lhs = fix (d, x.lhs);
      res.rhs = fix (d, x.rhs); return res; }
    expression operator ()(_min const& x) const
    { _min res; res.lhs = fix (d, x.lhs);
      res.rhs = fix (d, x.rhs); return res; }
  };

}//namespace detail

expression fix (std::string const& d, expression const& x) {
  return boost::apply_visitor (detail::fix_visitor (d), x);
}
You see what's going on there right? In order to fix an expression you have to recurse down into the expression replacing unfixed observations with fix expressions.

Ok, let's give it a go!

int main () {

  expression rate = observation ("heart rate");
  expression scale =  const_ (100);
  expression score = scale * fix ("t + 5", 
      min (const_ (0.2), max (rate / fix ("t", rate) - const_ (1.0), const_ (0.0))))
  ;

  std::cout << score << std::endl;

  return 0;
}
That program produces the output
100 * min (0.2, max(fix(t + 5, observation("heart rate")) / fix(t, observation("heart rate")) - 1, 0.0))
(note that the fixing to $t + 5$ has moved into the interior of the original expression as we'd hope).

So that's it - a basic expression algebra with placeholders. What's missing at this point is a function to simplify an expression in the context of a set of fixings. If we assume fixings are described as tuples of (1) the observation tag, (2) the time-stamp and (3) the observed value, the signature for that function might read something like,

expression simplify (
   std::vector<std::tuple<std::string, std::string, double>> const& fixings
  , expression const& xpr);
I'll leave the implementation of this as an exercise for the reader.

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}$.

OCaml:

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

Felix:

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

C++03:

#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), 
        detail::cumulative_moving_average)
  ).first;
}

C++11:

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

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 {};
  }
  else{
    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);
}

//Tests

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 
    else: 
      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
   else:
     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!

Labels

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

Saturday, February 7, 2015

Recursive lists in C++

Earlier this week, I had a need for a recursive list, that is, a list defined in terms of itself. I think, "back in the day" implementing a data structure of that sort would have been a snap for the everyday C programmer. Today, in this modern C++ world I found myself struggling a little and came to think that maybe the old ways are fading :)

For motivation, here's a couple of examples of the sort of thing I'm talking about.

(1) The list [0; 1; 0; 1; 0; 1; ...] is a list with a cycle in it. In OCaml you'd write that as let rec l = 0 :: 1 :: l.

(2) An interpreter using the technique of environments and closures can require an environment ((string * value) list) to contain a closure where the closure contains the environment. In OCaml you'd write let rec vars = (tag, V_closure (vars, xpr)) :: !env; env := vars.

Of course with pointers, it's not hard to implement recursive structures in C++. The trouble is having to concern yourself with their memory management due to the absence of garbage collection.

Alright, here is what I came up with. The code is pretty short.

#include <boost/variant.hpp>

#include <memory>
#include <stdexcept>

template <class T> struct node;
template <class T> using node_ptr=typename node<T>::node_ptr;
template <class T> using node_weak_ptr=typename node<T>::weak_ptr;
template <class T> using node_shared_ptr=typename node<T>::shared_ptr;
template <class T> struct ptr_t;
template <class T> using list=ptr_t<node<T>>;
template <class T> using list_ref=node_weak_ptr<T>;

template <class T> list<T> nil (); 
template <class T> bool empty (list<T> l);
template <class T> list<T> cons (T val, list<T> l);
template <class T> T& hd (list<T> l);
template <class T> list<T>& tl (list<T> l);
template <class T> list_ref<T> ref (list<T> src);
template <class T> bool is_ref (list<T> src);

The idea behind the implementation is generalize a pointer to node as a union with two variants, a shared pointer or a weak pointer.

template <class T> struct ptr_t :
  boost::variant <std::shared_ptr<T>, std::weak_ptr<T>> {
  typedef boost::variant <std::shared_ptr<T>, std::weak_ptr<T>> base;
  ptr_t () {}
  ptr_t (std::weak_ptr<T> p) : base (p) {}
  ptr_t (std::shared_ptr<T> p) : base (p) {}
};

template <class T>
struct node {
  typedef ptr_t<node> node_ptr;
  typedef std::weak_ptr<node> weak_ptr;
  typedef std::shared_ptr<node> shared_ptr;

  T data;
  node_ptr next;
};

This little bit of implementation detail comes up a couple of times so it's handy to factor it out.

namespace {
  //'get' at the raw pointer in the union of a smart/weak pointer
  template <class T>
  T* get (ptr_t<T> l) {
    if (std::shared_ptr<T>* p=
        boost::get<std::shared_ptr<T>>(&l)) {
      return p->get ();
    }
    return boost::get<std::weak_ptr<T>>(l).lock ().get ();
  }
}//namespace<anonymous>

The rest of the implementation is basically a set of "one-liners".

template <class T> list<T> nil (){ 
  return node_shared_ptr<T> (); 
}

template <class T> bool empty (list<T> l) { 
  return (get (l)) == nullptr; 
}

template <class T> list<T> cons (T val, list<T> l) {
  return node_shared_ptr<T> (new node<T>{val, l});
}

template <class T> T& hd (list<T> l) {   
  if (empty (l))
    throw std::runtime_error ("hd");
  return get (l) -> data;
}

template <class T> list<T>& tl (list<T> l) { 
  if (empty (l))
    throw std::runtime_error ("tl");
  return get (l) -> next; 
}

template <class T> bool is_ref (list<T> src) {
  return boost::get<list_ref<T>>(&src)!=nullptr;
}

template <class T> node_weak_ptr<T> ref (list<T> src)  {
  return node_weak_ptr<T>(boost::get<node_shared_ptr<T>>(src));
}

OK, well, that's about it. Let's see, regarding usage, (1) could be expressed like this

list<int> l = cons (0, cons (1, nil<int> ())); tl (tl (l)) = ref (l);
or, if we assume the existence of a 'last' function with an obvious definition, could be tidied up to read
list<int> l = cons (0, cons (1, nil<int> ())); tl (last (l)) = ref (l);
and (2) can be stated like this
typedef std::pair<std::string, value_t> p_t;
list<p_t> vars = node_shared_ptr<p_t>(new node<p_t>);
hd (vars) = std::make_pair (tag, V_closure {ref (vars), xpr});
tl (vars) = *env; 
*env = vars;