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 I think. 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 <gtest/gtest.h>

#include "sum_type.hpp"

#include <iostream>
#include <cstdlib>

//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>;

//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<T> some (T&& val) {
  return option<T>{constructor<some_t<T>>{}, std::forward<T> (val)};
}

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

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

//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.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.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 const& default_ (T const& x, option<T> const& u) {
  return u.match<T const&> (
    [](some_t<T> const& o) -> T const& { return o.data; },
    [=](none_t const&) -> T const& { 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.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.match<result_t>(
      [](none_t const&) { return none<t>(); }, 
      [=](some_t<T> const& o) { return k (o.data); }
  );
}

//Option monad 'unit'
template<class T>
option<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.match<option<t>>(
      [](none_t const&) { return none<t>(); }, 
      [=](some_t<T> const& o) { 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){ 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"

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

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

auto mul (int x) {
  return [=](int y) {

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

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