Saturday, December 7, 2013

Pipelining with the |> operator in OCaml

Pipelining with |> in OCaml

Reading through Real World OCaml (which I'm really enjoying and can heartily recommend by the way), I came across this clever definition for pipe notation.
let ( |> ) x f = f x
This little operator allows for (maybe misusing a term I think may have been coined by my friend MS,) "reverse application". To get a sense of the impact on readability that it may have, I applied it to the program of this earlier blog entry. At the heart of the program is the following little fragment.
match (filesin !root) with
| Some names ->
  let n=String.length !suffix in
  let pred e =
    let i = (String.length e - n) in
    (i >= 0) && (String.sub e i n) = !suffix
  in process_files (List.filter pred (Array.to_list names) )
Now, rewritten to use pipelining by virtue of |> we get this.
match (filesin !root) with
| Some names ->
  let n=String.length !suffix in
  let pred e =
    let i = (String.length e - n) in
    (i >= 0) && (String.sub e i n) = !suffix
  in names |> Array.to_list |> List.filter pred |> process_files 
 
I might be easily amused but... wow!

Coincidentally, reading from the OCaml PRO blog I learned of this related operator (also shown independently to me by my buddy -- thanks Stefano!).
let ( @@ ) f x = f x
This one can be used to good effect reducing the syntactic noise of lots of parentheses as per their example.
List.iter print_int @@ List.map (fun x -> x + 1 ) @@ [1; 2; 3]
(I assume it's obvious what this does and also, I added spaces in accordance with this advice).

Of course, if I'd had my wits about me, I'd have realized that I've known this operator all this time as $ from the Felix programming language... No excuse man, I mean, it's right there in the introductory tutorial!
println$ "Hello " + Env::getenv "USER";
Sorry 'bout being so slow to catch on Felix!

Now, I must admit I prefer $ to @@ but careful examination of the syntax reference seems to me that it won't be appropriate for OCaml (in that $ is left associative). As a final note, |> and @@ are "builtin" operators in OCaml compilers today.

Sunday, November 3, 2013

Matrix multiplication (functional approach in Python)

Matrix multiplication

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

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

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

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

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

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

Saturday, November 2, 2013

Gaussian elimination (functional approach in Python)

Gaussian elimination

The goal here is to implement simple Gaussian elimination in Python, in a functional style just using tuples.

We view (a, b, c) a row vector and interpret ((a,),(b,),(c,)) as a column vector.

These first two elementary operations (scaling a row by a scalar and subtracting one row from another) come easily.
def scale_row (u, lam): 
    return tuple (
      map (lambda ui: lam*ui, u))

def subtract_rows (u, v) :
  return tuple (
    map (lambda i : u[i] - v[i], range (0, len (u))))
These alone enable us to proceed directly to forward elimination of an augmented coefficient matrix.
def elimination_phase (t):
  def elimination_step (tp, k): 
    pr = tp[k]
    def g (e):
      (i, r) = e
      if i <= k : 
        return r
      else:
        lam = r[k]/pr[k]
        return subtract_rows (r, scale_row (pr, lam))
    return tuple (
             map (g, tuple (
                      zip (range (0, len (tp)), tp))))
  return functools.reduce (
              elimination_step, range (0, len (t)), t)
Vector dot products are a straight forward fold.
def dot (u, v):
  def f (acc, i):
    (a, (b,)) = (u[i], v[i])
    return acc + a*b
  return functools.reduce (f, range (0, len (u)), 0.0)
With this, the back substitution phase can be written like this.
def back_substitution_phase (t):
  n = len (t)
  def back_substitution_step (x, k):
    bk = (t[k])[n]
    akk = (t[k])[k]
    xk = bk/akk if k == n-1 \
      else (bk - dot ((t[k][k+1:])[0:n-k-1], x))/akk
    return ((xk,),) + x
  return functools.reduce (
    back_substitution_step, range (len (t)-1, -1, -1), ())
Can we test it? Yes we can!
#Solve Ax=b with 
#
#        6   -4     1
#  A =  -4    6    -1
#        1   -4     6
#
#and b = (-14, 36, 6)^T.

A = (
    ( 6.,  -4.,  1.,  -14.)
  , (-4.,   6., -4.,   36.)
  , ( 1.,  -4.,  6.,    6.))

print (str (back_substitution_phase (elimination_phase (A))))

#Should give (10, 22, 14)^T.

Term matching

Term matching

Extending one substitution with another can only work if the two substitutions are consistent with each other.
exception Term_match_exc

let extend_subst s1 s2 =
  let f acc (x, t)= 
    try 
      let u = List.assoc x acc in
      if t = u then acc
      else raise Term_match_exc
    with Not_found -> (x, t)::acc
in List.fold_left f s1 s2
A pattern match of a term u by a term t is a substitution that transforms t into u. If t is a variable s then the substitution (s, u) is an obvious pattern match for example. More generally though, if t is of the form f (t1, ..., tn), then
  • u must be of the form f (u1, ..., un);
  • for each i there must exist a match of ui with ti;
and all the matches must be mutually compatible for a solution to exist.
let term_match (t1, t2) =
  let rec term_match' subst = function
    | (Var v, t) -> extend_subst [v, t] subst
    | (t, Var v) -> raise Term_match_exc
    | (Term (f, l), Term (g, r)) ->
      if f = g then 
        List.fold_left term_match' subst (List.combine l r)
      else raise Term_match_exc
  in term_match' [] (t1, t2)
For example, a pattern match of h(y(z, k(x, s))) and h(y(l(m, n,o()), k(t(), u))) is {(x,t()),(z, l(m,n,o()), (s, u)}.

Saturday, October 19, 2013

Substitutions on terms with variables

Substitutions on terms with variables

This picks up on the earlier terms with variables blog entry. Remember the term_trav function? I'll recall it here but this time sharpen it up a bit with named parameters.
let rec term_trav ~f ~g ~x ~v = function
  | Term (a, tl) -> 
    let l = List.map (term_trav ~f ~g ~x ~v) tl in
    let res = List.fold_right g l x in
    f (a, res)
  | Var b -> v b

Application of substitutions

Substitutions are modeled as values of type ('b, ('a, 'b) term) list. Application of a substitution is about the finding of some of the variables of a term and replacing them with different terms to find the term's image under the substitution.
let apply_substitution subst t = 
  term_trav 
    ~f:(fun (f, l) -> Term (f, l)) 
    ~g:(fun x acc -> x::acc) 
    ~x:[] 
    ~v:(fun s -> try List.assoc s subst with _ -> Var s) 
    t

Of course, the definition of apply_substitution is equivalent to this one
let rec apply_substituion' subst = function
  | Term (f, ns) -> Term (f, List.map (apply_subst subst) ns)
  | Var x as v -> try List.assoc x subst with _ -> v
and I guess is more efficient too (in that a fold_right has been eliminated). 

For example, this program
let s = ["x", term_of_string "g(z)"] in
print_term (apply_substitution s (term_of_string "f (x, y, x)"))
shows that f (x, y, x) -> f (g(z), y, g(z)) under the substitution {(x, g(z))}.  

Composition of substitutions

If sig1 and sig2 are two substitutions, we can compute their composite sig1 o sig2 like this
let compose_subst sig1 sig2 =
  (List.map (fun (v, t) -> (v, apply_substitution sig1 t)) sig2)
  @(let vs = List.map fst sig2 in 
    List.filter (fun (x, t) -> not (List.mem x vs)) sig1)

That's quite subtle! Informally this is saying : first apply to the terms in sig2, the substitutions of sig1. Then, filter out those substitutions in sig1 that are in variables modified by sig2. The result is the concatenation of those lists. 

Here's an example. If sig1={(x, g(x,y))} and sig2={(y, h(x,z)), (x, k(x))} then we'd expect sig1 o sig2 = {(y, h(g(x,y), z)), (x, k(g(x, y)))}. The following program can be used to test our expectations.

 let sig1 = ["x", term_of_string "g (x, y)"] in
let sig2 = ["y", term_of_string "h (x, z)"; "x", term_of_string "k(x)"] in
let subst=compose_subst sig1 sig2 in
print_string "y -> " ; print_term (List.assoc "y" subst) ; print_newline ();
print_string "x -> " ; print_term (List.assoc "x" subst) ; print_newline ()

Saturday, October 12, 2013

Tuple matching

Tuple matching

Here is an algorithm that can be used to match tuples. For example, on matching (a, b, (c, d)) against (1, 2, (3, 4)) we'd hope to get the set (a = 1, b = 2, c = 3, d = 4).
module type S = sig

  type expr =
  | E_var of string
  | E_const of int
  | E_tuple of expr list

  val const : int -> expr
  val var : string -> expr
  val tuple : expr list -> expr

  val tuple_match : (string * expr) list -> expr -> expr -> (string * expr) list
  val tuple_of_expr : expr -> expr

end

module Expr : S = 
struct

  type expr =
  | E_var of string
  | E_const of int
  | E_tuple of expr list

  let var s = E_var s
  let const i = E_const i
  let tuple t = E_tuple t

  let as_tuple = function
    | E_tuple l as t -> t
    | _ ->  failwith "Tuple expected"

  let rec tuple_match acc x y =
    match x with
    | E_var (s) -> (s, y)::acc
    | E_tuple ((h::t)) ->
      let (E_tuple l) = as_tuple y in
      let acc = tuple_match acc (tuple t) (tuple (List.tl l)) in
      tuple_match acc h (List.hd l)
    | E_tuple [] -> acc
    | _ as unk -> failwith "Match failure"

end
For example,
(* (a, b, t) |- (1, 2, (3, 4)) *)
let x = Expr.tuple 
  [Expr.var "a"; Expr.var "b"; Expr.var "t" ] in
let y = Expr.tuple [Expr.const 1; Expr.const 2; 
   Expr.tuple [Expr.const 3; Expr. const 4] ] in 
Expr.tuple_match [] x y
;;
in the top-level prints
- : (string * Expr.expr) list =
[("a", Expr.E_const 1); ("b", Expr.E_const 2);
 ("t", Expr.E_tuple [Expr.E_const 3; Expr.E_const 4])]
whereas
(* (a, b, (c, d)) |- (1, 2, (3, 4))*)
let x = Expr.tuple 
  [Expr.var "a"; Expr.var "b"; 
   Expr.tuple [Expr.var "c"; Expr.var "d"]] in
let y = Expr.tuple 
  [Expr.const 1; Expr.const 2; 
   Expr.tuple [Expr.const 3; Expr. const 4]] in 
Expr.tuple_match [] x y
;;
results in this
- : (string * Expr.expr) list =
[("a", Expr.E_const 1); ("b", Expr.E_const 2); ("c", Expr.E_const 3);
 ("d", Expr.E_const 4)]

Saturday, September 28, 2013

Terms with Variables

Cosineau and Mauny present a type suitable for modeling terms with variables convenient for applications involving variables and substitutions. It is defined like this:
(*term_types.ml*)
type (α, β) term = 
| Term of α * (α, β) term list
| Var of β

Analyzing syntax

To proceed, let us restrict our attention to the type (string, string) term. The concrete syntax we use for such expressions is parenthesized.
  • a (t1, ..., tn) is the syntax of Term (a, [t1; ...; tn])
We disambiguate variables from constants by representing constants as functions applied to an empty list of arguments
  • a (b (), c) is the syntax of Term ("a", [Term ("b", []); Var "c"])
We start towards a function to analyze syntax. First the lexer.
(*lexer.mll*)
{
  open Parser  (*The type token is defined in parser.mli.*)
}

let alpha=['a'-'z' 'A'-'Z']

rule token = parse
    [' ' '\t' '\n']         (* skip blanks *) { token lexbuf }
  | '('                                             { LPAREN }
  | ')'                                             { RPAREN }
  | ','                                              { COMMA }
  | ((alpha)(alpha)*) as s                         { STRING s}
  | eof                                                { EOI }
Now the parser.
/*parser.mly*/

%token  STRING
%token LPAREN RPAREN COMMA EOI

%start main
%type <(string, string) Term_types.term> main
%%
main:
 | exp EOI { $1 }
 ;
exp:
 | tag exp_more           
     { 
       match $2 with
       | None -> Term_types.Var $1
       | Some l -> Term_types.Term ($1, l)
     }
 ;
exp_more:
 | /*nothing*/ { None }
 | LPAREN exp_list RPAREN { Some $2 }
;
exp_list:
 | /*nothing*/ { [] }
 | exp exp_list_more { $1::$2 }
 ;
exp_list_more:
 | /*nothing*/ {[]}
 | COMMA exp exp_list_more { $2::$3 }
 ;
tag:
 | STRING { $1 }
 ;
With these pieces in place we can write the analyzer term_of_string like so
(*term.ml*)

let term_of_string s =
  let parse_buf lexbuf =
    try Parser.main Lexer.token lexbuf
    with 
    | Parsing.Parse_error ->
      begin
        let curr = lexbuf.Lexing.lex_curr_p in
        let line = curr.Lexing.pos_lnum in
        let cnum = curr.Lexing.pos_cnum - curr.Lexing.pos_bol in
        let tok = Lexing.lexeme lexbuf in
        raise 
          (Failure
             (Printf.sprintf 
                "file \"\", line %d, character %d\n\Error : Syntax error \"%s\"" line cnum tok))
      end
  in parse_buf (Lexing.from_string s)
Reverse parsing is done with string_of_term:
let rec string_of_term (t:(string, string)term) : string = 
  let string_of_list (f:α -> string) (l:α list) : string = 
    "(" ^ String.concat "," (List.map f l) ^ ")"  in
  match t with
  | Term (s, tl) -> Printf.sprintf "%s" (s^(string_of_list string_of_term tl))
  | Var s -> s
A simple function for printing a term then goes like this:
let print_term (t:(string, string) term) : unit = Printf.printf "%s" (string_of_term t)
For example, this program let tt = term_of_string "a(b(), c)" in print_term tt prints "a(b(), c)" recovering the input as we'd hope.

Traversal

The function term_trav is a workhorse for doing computations over terms. It generalizes all functions defined by recursion on the structure of a term.
(* [term_trav] is a generalization of all functions defined by
   recusion on the structure of an [(α, β) term] - a homomorphism.

   The definition of [term_trav] below gives the following signature.

   {[
   val term_trav : (α * β -> γ) ->
   (γ -> β -> β) -> β -> (δ -> γ) -> (α, δ) term -> γ
   ]}

   Proof:

   Assume [term_trav] to operate on arguments of type [(α, δ)
   term]. It "returns" either via [f] or [v], and so [f], [v] and
   [term_trav] must share the same return type, [γ] say.

   Since the intermediate list results from recursive calls to
   [term_trav] then it must have type [γ list]. This means [g] must
   be of type [γ -> β -> β] for some type [β] which fixes [x] to
   [β] whilst completing [f] as [α * β -> γ].

   Lastly [v] takes arguments of type [δ] meaning it types to [δ ->
   γ].

*)
let rec term_trav f g x v = function
  | Term (a, tl) -> 
    let l = List.map (term_trav f g x v) tl in
    let res = List.fold_right g l x in
    f (a, res)
  | Var b -> v b
For a trivial example, we can express a function to compute the "size" of a term this way:
let term_size (t:(α, β) term) : int = 
  (term_trav (fun (_, s) -> s + 1) (fun x y -> x + y) 0 (fun _ -> 1))
    t
I appreciate there's an efficiency gain to be had by eliminating the intermediate list in term_trav but the resulting function definition is more challenging to reason about (well, at least for me).

Saturday, September 21, 2013

Real world OCaml

Real world OCaml

For many tasks, I find OCaml as useful as Python and like "a better C".

This program matches files in a directory with a provided suffix and invokes a program on them.
(*e.g.
  Match files in directory with given suffix and invoke a program on
  them.

$ ./build/bin/filesin.exe -d \
    c:/project/examples -p \
    c:/project/someprog.exe
*)

(*Globals*)
let version,root,prog,suffix=ref false,ref "",ref "",ref ""

let (read_args:unit -> unit) = fun () ->
  let specification =
    [("-v", Arg.Set version, "Print the version number") ;
     ("-d", Arg.String (fun s -> root := s), "Directory") ;
     ("-s", Arg.String (fun s -> suffix := s), "Suffix") ;
     ("-p", Arg.String (fun s -> prog := s), "Path to the prog") ;
    ]
  in Arg.parse specification
  (fun s -> Printf.printf "Warning : Ignoring unrecognized argument \"%s\"\n" s)
  ("Usage : "^(Sys.argv.(0))^" options")

let (filesin:string -> string array option) = fun s ->
  if not (Sys.file_exists s)
       || not (Sys.is_directory s) then
    None
  else (Some (Sys.readdir s))

let (read_and_close:Unix.file_descr -> string)  = fun fdr ->
  let buf = Buffer.create 1024 in
  let bytes_to_read = ref true in
  while !bytes_to_read do
    let s  = String.create 1024 in
    let num_bytes_read = Unix.read fdr s 0 124 in
    Buffer.add_string buf (String.sub s 0 num_bytes_read) ;
    if num_bytes_read < 1024 then bytes_to_read := false ;
    ()
  done ;
  Unix.close fdr ;
  Buffer.contents buf

let (process_file:string -> unit) = fun f ->
  let p = Filename.concat (!root) f in
  let cmd = Printf.sprintf "%s -l -e %s" (!prog) p in
  (*I tried to use anonymous pipes here but that seems broken on
    Windows when [Unix.system] is involved*)
  let tf=Filename.temp_file "tmp" "" in
  let fdw = Unix.openfile tf [Unix.O_WRONLY; Unix.O_TRUNC; Unix.O_CREAT] 0o666 in
  let tmp = Unix.dup Unix.stdout in
  let stat = (
    Unix.dup2 fdw Unix.stdout ;
    Printf.printf "\n*** - %s:\n" p; flush stdout;
    Unix.system cmd  (*Invoke the intpreter on f*)
  )
  in
  Unix.dup2 tmp Unix.stdout ; Unix.close tmp; Unix.close fdw;
  match stat with
  | Unix.WEXITED 127 ->
    Printf.printf "The command \"%s\" could not be executed" cmd
  | _ ->
    let fdr = Unix.openfile tf [Unix.O_RDONLY] 0o666 in
    Printf.printf "%s" (read_and_close fdr) ; flush stdout (*; Sys.remove tf*)

let (process_files:string list -> unit) = fun fs ->
  List.iter (fun f -> process_file f) fs

let (main:unit) =
  read_args () ; Printf.printf "%s" !prog ;
  if !version then print_string "0.0.0\n"
  else
    if !suffix = "" then
    raise (Failure
      "Missing argument : suffix") ;
   if not (Sys.file_exists (!prog)) then
      raise (Failure (Printf.sprintf
        "Bad argument : The file '%s' does not exist" !prog)) ;
   try
     match (filesin !root) with
     | Some names ->
       let n=String.length !suffix in
       let pred e =
         let i = (String.length e - n) in
         (i >= 0) && (String.sub e i n) = !suffix
       in process_files (List.filter pred (Array.to_list names) )
     | None -> ()
   with
   | _  as e -> Printf.printf "Exception : %s\n" (Printexc.to_string e)

Friday, September 13, 2013

Sieve of Eratosthenes

Sieve of Eratosthenes

Find all the prime numbers less than or equal to k.

Python

import functools

k=1000

def sieve (k):

  def f (p, lst):
    s = list (filter (lambda i : i%p != 0 or i/p == 1, lst))
    r = list (filter (lambda i : i > p, s))
    if len (r) == 0:
        return s
    else:
        return f (r[0], s)

  return f (2, range (2, k + 1))

print (str(sieve (k)))

Sunday, September 8, 2013

Sub-arrays of maximal sum

Sub-arrays of maximal sum

Another puzzle problem from this blog. "You are given an array of n integers, each of which may be positive, negative or zero. Give an algorithm to identify the start and end index, i and j, of the interval whose elements form the maximal sum of all possible intervals. Assume j >=i. e.g. {1 3 -8 2 -1 10 -2 1} -> i=3, j=5 - sum=11.

Python

'''Find sub-arrays of maximal sum

   e.g. {1 3 -8 2 -1 10 -2 1} -> i=3, j=5 - sum=11

   http://basicalgos.blogspot.com/2012/05/52-find-sub-array-of-max-sum.html
'''

import functools

#Compute all the intervals of the sequence t per the definition
#above. e.g. for [1, 2, 3] compute [[1], [1, 2], [1, 2, 3], [2],
#[2, 3], [3]]. 
#Decorate the intervals with their start and end indices and the
#interval sum.
def intervals (t):
  def f (acc, i):
    def row (i, t):
      s = t[i:]
      def f (acc, j):
         indices=(i, i+j-1)
         interval=s[0:j]
         return acc + [(indices, interval, sum (interval, 0))]
      return functools.reduce (f, list (range (1, len (s) + 1)), [])
    return acc + row (i, t)
  return functools.reduce (f, list (range (0, len (t))), [])
  
#x, the input sequence.
x=[1, 3, -8, 2, -1, 10, -2, 1]

#Print the intervals of x with maximal sum.
res = intervals (x)
def f (m, t): 
  (_,_,n)=t ; return max (m, n)
maxsum = functools.reduce (f, res, 0)
for t in res:
  (indices, interval, total)=t
  if (total == maxsum):
    print (str(t)) #An interval of maximal sum.

C++

#include <vector>
#include <numeric>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <functional>

typedef std::pair<std::size_t, std::size_t> indicies_t;
inline std::size_t& fst (indicies_t& i) { return i.first; }
inline std::size_t& snd (indicies_t& i) { return i.second; }
inline std::size_t fst (indicies_t const& i) { return i.first; }
inline std::size_t snd (indicies_t const& i) { return i.second; }
void print_indicies (indicies_t const& p);

typedef std::vector<int> interval_t;
inline std::size_t len (interval_t const& i) { return i.size(); }
interval_t take (int k, interval_t const& lst);
interval_t drop (int k, interval_t const& lst);
inline int sum (interval_t const t)
{ 
  int tot=std::accumulate (t.begin(), t.end(), 0, std::plus<int>());

  return tot;
}
interval_t range (int begin, int end);
void print_interval (interval_t const& s);

struct interval_info_t
{
  indicies_t indicies; //(0, 2)
  interval_t interval; //e.g. [1, 2, 3]
  int total; //6
};
void print_interval_info (interval_info_t const& i);

typedef std::vector<interval_info_t> interval_info_set_t;

struct inner_f
{
  interval_t s;
  std::size_t i;

  inner_f (interval_t const& s, std::size_t i) 
    : s (s), i (i) 
  {}

  interval_info_set_t operator() 
  (interval_info_set_t acc, std::size_t j)
  {
    indicies_t indicies;
    fst (indicies) = i;
    snd (indicies) = i + j - 1;
    interval_t interval = take (j, s);
    int total = sum (interval);

    interval_info_t interval_info;
    interval_info.indicies = indicies;
    interval_info.interval = interval;
    interval_info.total = total;

    acc.push_back (interval_info);
  
    return acc;
 }
};

interval_info_set_t row (std::size_t i, interval_t const& t)
{
  interval_t s = drop (i, t);
  interval_t rg = range (1, len (s) + 1);
  
  return std::accumulate (
    rg.begin(), rg.end(), interval_info_set_t(), inner_f (s, i));
}

struct outer_f
{
  interval_t t;
  outer_f (interval_t const t) : t (t) {}

  interval_info_set_t operator() 
    (interval_info_set_t acc, std::size_t i) const
  {
    interval_info_set_t res=row (i, t);
    std::copy (res.begin(), res.end(), std::back_inserter (acc));
  
    return acc;
  }
};

interval_info_set_t intervals (interval_t const& t)
{
  interval_t rg = range (0, len (t));
  
  return std::accumulate (
    rg.begin(), rg.end(), interval_info_set_t (), outer_f (t));
}

inline int find_max (int m, interval_info_t const& t)
{
  return (std::max)(m, t.total);
}

//--driver

int main ()
{
  int data[] = {1, 3, -8, 2, -1, 10, -2, 1};
  interval_t input(data, data + sizeof(data)/sizeof(int));
  interval_info_set_t ivals = intervals (input);
  int maxsum = std::accumulate (ivals.begin(), ivals.end(), 0, find_max);

  for (std::size_t i = 0; i < ivals.size(); ++i)
    {
      if (ivals[i].total==maxsum)
 print_interval_info (ivals[i]); //An interval of maximal sum.
    }
  
  return 0;
}

//-- tools

interval_t take (int k, interval_t const& lst)
{
  if (k <= 0) return interval_t ();
  if (lst.size() == std::size_t (0)) return interval_t ();

  interval_t t = take (k-1, interval_t (++(lst.begin()), lst.end()));
  interval_t res (t.size()+1);
  res[0] = lst[0];
  for (std::size_t i = 0; i < len (t); ++i)
    res [i+1] = t[i];

  return res;
}

interval_t drop (int k, interval_t const& lst)
{
  if (k <= 0) return lst;
  if (lst.size () == 0) return interval_t ();

  return drop (k - 1, interval_t (++(lst.begin()), lst.end()));
}

interval_t range (int begin, int end)
{
  interval_t buf;
  for(int i = begin; i < end; ++i) buf.push_back (i);

  return buf;
}

void print_indicies (indicies_t const& p)
{
  std::cout << "(" << fst (p) << ", " << snd (p) << ")";
}

void print_interval (interval_t const& s)
{
  std::cout << "["; 
  std::ostream_iterator<int> out_it (std::cout, ", ");
  std::copy (s.begin(), s.end(), out_it);
  std::cout << "],"; 
}

void print_interval_info (interval_info_t const& i)
{
  std::cout << "(";
  print_indicies (i.indicies);
  print_interval (i.interval);
  std::cout << i.total;
  std::cout << ")\n";
}

Saturday, September 7, 2013

Permutations of a list

Permutations of a list

Find all the permutations of a list. For example, given the list [1,2, 3] compute all of [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1],[3, 1, 2], [3, 2, 1]].

The reasoning here goes like this. Suppose I know how to compute all the permutations of a list that start with the element at position k. Then the list of all permutations can be obtained by calling that function with k = 0, ..., N-1 where N is the length of the list. How do you compute the permutations of a list that start with the element at position k? Well, you compute all the permutations of the reduced list that doesn't contain the kth element and prepend them all with the kth element.

Python

def take(k, lst):
   '''The list slice, lst[0:k].
  
      >>> print (str(take(0, [1, 2, 3])))
      []
      >>> print (str(take(1, [1, 2, 3])))
      [1]
      >>> print (str(take(2, [1, 2, 3])))
      [1, 2]
      >>> print (str(take(3, [1, 2, 3])))
      [1, 2, 3]
      >>> print (str(take(4, [1, 2, 3])))
      [1, 2, 3]
  
   '''
   if(k <= 0):
       return []
   if(len(lst) == 0):
       return []
   return [lst[0]]+take(k-1, lst[1:])
  
def drop(k, lst):
    '''The list slice, lst[k:].
  
       >>> print (str(drop(0, [1, 2, 3])))
       [1, 2, 3]
       >>> print (str(drop(1, [1, 2, 3])))
       [2, 3]
       >>> print (str(drop(2, [1, 2, 3])))
       [3]
       >>> print (str(drop(3, [1, 2, 3])))
       []
       >>> print (str(drop(4, [1, 2, 3])))
       []
  
    '''
    if(k <= 0):
        return lst
    if(len(lst)==0):
        return []
    return drop(k-1, lst[1:])
  
def find_all_permutations(lst):
  '''Find all permutations of a list.
  
     >>> print (str(find_all_permutations([1, 2, 3])))
     [[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
  
  '''
  def find_k_permutations(k, lst):
    '''Find all permutations of the given list that start with the
       element at index k.
  
    '''
    x = lst[k]
    listp=take(k, lst)+drop(k+1, lst)
    ps=find_all_permutations(listp)
    s = list (map (lambda perm : [x] + perm, ps))
 
    return s
  
  def loop(k, acc, lst):
      if(k == len(lst)):
          return acc
      else:
          return loop(k+1, acc + find_k_permutations(k, lst), lst)
  
  return ([[]], loop(0, [], lst))[len(lst)!=0]
  
if __name__ == "__main__":
    import doctest
    doctest.testmod()

C++

#include <iostream>
#include <iterator>
#include <algorithm>
#include <vector>

typedef std::vector<int> perm_t;
typedef std::vector<int>::iterator perm_iterator_t;
typedef std::vector<int>::const_iterator perm_const_iterator_t;

typedef std::vector<perm_t> perm_set_t;
typedef std::vector<perm_t>::iterator perm_set_iterator_t;
typedef std::vector<perm_t>::const_iterator perm_set_const_iterator_t;

//-- print functions

void print_perm (perm_t const& s)
{
  std::cout << "["; 
  std::ostream_iterator<int> out_it (std::cout, ", ");
  std::copy (s.begin(), s.end(), out_it);
  std::cout << "],"; 
}

void print_perm_set (perm_set_t const& ts)
{
  std::cout << "[";
  perm_set_const_iterator_t begin=ts.begin (), end = ts.end ();
  perm_set_const_iterator_t t=begin;
  while (t != end)
    {
      print_perm (*t);
      ++t;
    }
  std::cout << "]";
}

//-- take, drop

perm_t take (int k, perm_t const& lst)
{
  if (k <= 0)
    {
      return perm_t ();
    }
  if (lst.size() == std::size_t (0))
    {
      return perm_t ();
    }

  perm_t t = take (k-1, perm_t (++(lst.begin()), lst.end()));
  perm_t res (t.size()+1);
  res[0] = lst[0];
  for (std::size_t i = 0; i < t.size (); ++i)
    res [i+1] = t[i];

  return res;
}

perm_t drop (int k, perm_t const& lst)
{
  if (k <= 0)
    {
      return lst;
    }

  if (lst.size () == 0)
    {
      return perm_t ();
    }

  return drop (k - 1, perm_t (++(lst.begin()), lst.end()));
}

//-- helpers

perm_set_t find_all_permutations (perm_t const& lst);

//Prepend x to a perm.

struct prepend
{
  int x;

  prepend (int x) : x (x)
  {}

  perm_t operator ()(perm_t const& p) const
  {
    perm_t r (p.size() + 1);
    r[0] = x;
    for (std::size_t i = 0; i < p.size(); ++i)
      r[i + 1] = p[i];

    return r;
  }
};

//Find all permutations of the given list that
//start with the element at index k.

perm_set_t find_k_permutations (int k, perm_t const& lst)
{
  int x = lst[k];
  perm_t left = take (k, lst);
  perm_t right = drop (k + 1, lst);
  perm_t listp;
  std::copy (left.begin(), left.end(), std::back_inserter (listp));
  std::copy (right.begin(), right.end(), std::back_inserter (listp));
  perm_set_t ps = find_all_permutations (listp);
  perm_set_t s;
  std::transform (ps.begin(), ps.end(), std::back_inserter(s), prepend (x));

  return s;
}

perm_set_t loop (int k, perm_set_t acc, perm_t lst)
{
  if (k == lst.size())
    return acc;
  else
    {
      perm_set_t t = find_k_permutations (k, lst);
      std::copy (t.begin (), t.end (), std::back_inserter (acc));
      
      return loop (k+1, acc, lst);
    }
}

perm_set_t find_all_permutations (perm_t const& lst)
{
  if (lst.size() == 0)
    {
      perm_set_t res;
      res.push_back (perm_t ());//The empty sequence.

      return res;
    }

  return loop (0, perm_set_t (), lst);
}

//-- test

int main()
{
  int data[] = {1, 2, 3};
  perm_t input (data, data + sizeof(data)/sizeof(int));

  print_perm_set (find_all_permutations (input));

  return 0;
}

Longest increasing sub-sequence

Longest increasing sub-sequence


According to this blog, "The longest increasing sub-sequence problem is to find a sub-sequence of a given sequence in which the sub-sequence elements are in sorted order, lowest to highest, and in which the sub-sequence is as long as possible". As an example, for 0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15 they give one longest increasing sub-sequence as 0, 2, 6, 9, 13, 15.

Python

'''Find a longest increasing subsequence
http://basicalgos.blogspot.com/2012/03/37-longest-increasing-sequence-lis.html

'''
import functools

def ics (s):  #s is s sequence
  """Find the set of increasing subsequences of s.

  """
  if len (s) == 1 :
      return [s] #s itself is an ics of length 1

  #This is the key point.
  ts = ics (s[1:]) #ts are the ics values of s - it's first element

  #Extend ts with any new subsequences that can be formed by
  #prepending s[0] to t for t in ts.
  def f (acc, t) :
    if s[0] < t[0]:
       nt = [s[0]]+t
       return acc + [nt]
    else:
        return acc
  res = functools.reduce (f, ts, ts) + [[s[0]]]

  return res

#s, the input sequence.

s=[0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15]

#Print all longest increasing subsequences of s.

res = ics (s) #Find every ics of s
lens=list (map (lambda x : len (x), res)) #Now their lengths
lismax= \
  functools.reduce (\
    lambda acc, y: max (acc, y), lens, 0) #The max length
for f in res:
    if len(f) == lismax: print (str(f)) #An ics of max length

C++

/*
Find a longest increasing subsequence
http://basicalgos.blogspot.com/2012/03/37-longest-increasing-sequence-lis.html
 */
#include <iostream>
#include <vector>
#include <iterator>
#include <algorithm>
#include <numeric>

typedef std::vector<int> ics_t;
typedef ics_t::const_iterator ics_const_iterator_t;
typedef ics_t::iterator ics_iterator_t;

typedef std::vector<ics_t> ics_set_t;
typedef ics_set_t::const_iterator ics_set_const_iterator_t;
typedef ics_set_t::iterator ics_set_iterator_t;

//-- print functions

void print_ics (ics_t const& s)
{
  std::cout << "["; 
  std::ostream_iterator<int> out_it (std::cout, ", ");
  std::copy (s.begin(), s.end(), out_it);
  std::cout << "],"; 
}

void print_ics_set (ics_set_t const& ts)
{
  std::cout << "[";
  ics_set_const_iterator_t begin=ts.begin (), end = ts.end ();
  ics_set_const_iterator_t t=begin;
  while (t != end)
    {
      print_ics (*t);
      ++t;
    }
  std::cout << "]";
}

//-- Calculate the length of an ics

std::size_t len_ics (ics_t const& s)
{
  return s.size ();
}

//-- Calculate the max ics length in a set

std::size_t max_ics_len (ics_set_t const& s)
{
  std::vector<std::size_t> lens;
  std::transform (
    s.begin (), s.end (), std::back_inserter (lens), len_ics);

  return std::accumulate (
        lens.begin (), lens.end (), 0, std::max<std::size_t>);
}

//  Extend ts with any new subsequences that can be formed by
//  prepending s[0] to t for t in ts.

struct f // f : ics_set_t * ics_t -> ics_set_t
{
  int s0; 

  f (int s0) : s0 (s0) 
  {}
  ics_set_t operator()(ics_set_t acc, ics_t const& t) const
  {
    if (s0 < t[0])
      {
        ics_t nt (t.size() + 1);
        nt[0] = s0;
        for (std::size_t i=0; i < t.size(); ++i)
          nt[i+1] = t[i];
        acc.push_back (nt);

        return acc;
      }
    else
      {
        return acc;
      }
  }
};

//Find the set of increasing subsequences of s.
ics_set_t ics (ics_t s)
{
  if (s.size () == 1)
    {
      ics_set_t res;
      res.push_back (s); //s itself is an ics of length 1

      return res;
    }

  //This is the key point.
  ics_set_t ts = ics (ics_t(++(s.begin()), s.end())); 
  //ts are the ics values of s - it's first element

  //Extend ts with any new subsequences that can be formed by
  //prepending s[0] to t for t in ts.
  ics_set_t t=std::accumulate (ts.begin (), ts.end(), ts, f(s[0]));
  t.push_back (ics_t(1, s[0]));

  return t;
}

int main ()
{
  int data[] = 
  {
    0, 8, 4, 12, 2, 10, 6, 14, 1,9, 5, 13, 3, 11, 7, 15
  };
  ics_t s (data, data + sizeof(data)/sizeof(int));

  ics_set_t res = ics (s);
  std::size_t n=max_ics_len (res);
  for (std::size_t i = 0; i < res.size (); ++i)
    {
      ics_t const& x=res[i];
      if (x.size () == n)
      {
        print_ics (x);
        std::cout << std::endl;
      }
    }
  
  return 0;
}

Saturday, July 13, 2013

Static OCaml embedding on Win32

A procedure for OCaml embedding in a DLL at least is documented on the web in this thread. As far as I know though, there isn't any explicit information anywhere on how to embed OCaml statically on Win32. This tip for how to do it was taught to me by Alain Frisch. The build commands are given as CMake rules. First we compile our OCaml sources and have ocamlopt.opt output an object file of the result.
add_custom_command( 
  OUTPUT foo_ocaml.obj
  DEPENDS foo.mli foo.ml
  COMMAND ${OCAMLOPT_OPT} -o foo_ocaml.obj
          foo.mli
          foo.ml
  )
We package the output object into a library archive with this:
add_library (foo STATIC foo_ocaml.obj)
set_target_properties (foo PROPERTIES LINKER_LANGUAGE CXX)
The C wrapper is compiled with the C compiler into an archive with a file flexdll_stubs_c.c:
add_library (foo_api
  STATIC
  foo.c
  flexdll_stubs_c.c
)
flexdll_stubs_c.c provides "what it says on the tin":
void flexdll_dump_exports(void* u)
{ 
  (void)u; 
}

void *flexdll_dlopen(char const* file, int mode)   
{ 
  (void)file; (void)mode; return (void*)0; 
}

void flexdll_dlclose(void* u)                                                      
{ 
  (void)u; 
}

void* flexdll_dlsym(void* u, const char * n)             
{ 
  (void)u; (void)n; return (void*)0; 
}

char* flexdll_dlerror()                                        
{ 
  static char* flexdll_error_buffer = "flexdll is not availble";

  return flexdll_error_buffer; 
}
Executable targets link with the two archives above and the OCaml runtime via libasmrun.lib.

Saturday, June 22, 2013

Maybe

There are different approaches to the issue of not having a value to return. One idiom to deal with this in C++ is the use of boost::optional<T> or std::pair<bool, T>.
class boost::optional<T> //Discriminated-union wrapper for values.

Maybe is a polymorphic sum type with two constructors : Nothing or Just a.
Here's how Maybe is defined in Haskell.

  {- The Maybe type encapsulates an optional value. A value of type
  Maybe a either contains a value of type a (represented as Just a), or
  it is empty (represented as Nothing). Using Maybe is a good way to
  deal with errors or exceptional cases without resorting to drastic
  measures such as error.

  The Maybe type is also a monad.
  It is a simple kind of error monad, where all errors are
  represented by Nothing. -}

  data Maybe a = Nothing | Just a

  {- The maybe function takes a default value, a function, and a Maybe
  value. If the Maybe value is Nothing, the function returns the default
  value. Otherwise, it applies the function to the value inside the Just
  and returns the result. -}

  maybe :: b -> (a -> b) -> Maybe a -> b
  maybe n _ Nothing  = n
  maybe _ f (Just x) = f x

I haven't tried to compile the following OCaml yet but I think it should be roughly OK.
 type 'a option = None | Some of 'a  ;;

  let maybe n f a =
    match a with
      | None -> n
      | Some x -> f x
      ;;

Here's another variant on the Maybe monad this time in Felix. It is applied to the problem of "safe arithmetic" i.e. the usual integer arithmetic but with guards against under/overflow and division by zero.

  union success[T] =
    | Success of T
    | Failure of string
    ;

  fun str[T] (x:success[T]) =>
    match x with
      | Success ?t => "Success " + str(t)
      | Failure ?s => "Failure " + s
    endmatch
    ;

  typedef fun Fallible (t:TYPE) : TYPE => success[t] ;

  instance Monad[Fallible]
  {
    fun bind[a, b] (x:Fallible a, f: a -> Fallible b) =>
      match x with
        | Success ?a => f a
        | Failure[a] ?s => Failure[b] s
      endmatch
      ;

    fun ret[a](x:a):Fallible a => Success x ;
  }

  //Safe arithmetic.

  const INT_MAX:int requires Cxx_headers::cstdlib ;
  const INT_MIN:int requires Cxx_headers::cstdlib ;

  fun madd (x:int) (y:int) : success[int] =>
    if x > 0 and y > (INT_MAX - x) then
        Failure[int] "overflow"
    else
      Success (y + x)
    endif
    ;

  fun msub (x:int) (y:int) : success[int] =>
    if x > 0 and y < (INT_MIN + x) then
      Failure[int] "underflow"
    else
      Success (y - x)
    endif
    ;

  fun mmul (x:int) (y:int) : success[int] =>
    if x != 0 and y > (INT_MAX / x) then
      Failure[int] "overflow"
    else
      Success (y * x)
    endif
    ;

  fun mdiv (x:int) (y:int) : success[int] =>
      if (x == 0) then
          Failure[int] "attempted division by zero"
      else
        Success (y / x)
      endif
      ;

  //--
  //
  //Test.

  open Monad[Fallible] ;

  //Evalue some simple expressions.

  val zero = ret 0 ;
  val zero_over_one = bind ((Success 0), (mdiv 1)) ;
  val undefined = bind ((Success 1),(mdiv 0)) ;
  val two = bind((ret 1), (madd 1)) ;
  val two_by_one_plus_one = bind (two , (mmul 2)) ;

  println$ "zero = " + str zero ;
  println$ "1 / 0 = " + str undefined ;
  println$ "0 / 1 = " + str zero_over_one ;
  println$ "1 + 1 = " + str two ;
  println$ "2 * (1 + 1) = " + str (bind (bind((ret 1), (madd 1)) , (mmul 2))) ;
  println$ "INT_MAX - 1 = " + str (bind ((ret INT_MAX), (msub 1))) ;
  println$ "INT_MAX + 1 = " + str (bind ((ret INT_MAX), (madd 1))) ;
  println$ "INT_MIN - 1 = " + str (bind ((ret INT_MIN), (msub 1))) ;
  println$ "INT_MIN + 1 = " + str (bind ((ret INT_MIN), (madd 1))) ;

  println$ "--" ;

  //We do it again, this time using the "traditional" rshift-assign
  //syntax.

  syntax monad //Override the right shift assignment operator.
  {
    x[ssetunion_pri] := x[ssetunion_pri] ">>=" x[>ssetunion_pri] =># "`(ast_apply ,_sr (bind (,_1 ,_3)))";
  }
  open syntax monad;

  println$ "zero = " + str (ret 0) ;
  println$ "1 / 0 = " + str (ret 1 >>= mdiv 0) ;
  println$ "0 / 1 = " + str (ret 0 >>= mdiv 1) ;
  println$ "1 + 1 = " + str (ret 1 >>= madd 1) ;
  println$ "2 * (1 + 1) = " + str (ret 1 >>= madd 1 >>= mmul 2) ;
  println$ "INT_MAX = " + str (INT_MAX) ;
  println$ "INT_MAX - 1 = " + str (ret INT_MAX >>= msub 1) ;
  println$ "INT_MAX + 1 = " + str (ret INT_MAX >>= madd 1) ;
  println$ "INT_MIN = " + str (INT_MIN) ;
  println$ "INT_MIN - 1 = " + str (ret INT_MIN >>= msub 1) ;
  println$ "INT_MIN + 1 = " + str (ret INT_MIN >>= madd 1) ;
  println$ "2 * (INT_MAX/2) = " + str (ret INT_MAX >>= mdiv 2 >>= mmul 2 >>= madd 1) ; //The last one since we know INT_MAX is odd and that division will truncate.
  println$ "2 * (INT_MAX/2 + 1) = " + str (ret INT_MAX >>= mdiv 2 >>= madd 1 >>= mmul 2) ;

  //--
That last block using the <<= syntax produces (in part) the following output (the last two print statments have been truncated away -- the very last one produces an expected overflow).

Tuesday, June 4, 2013

Powerset

This is an algorithm to compute a powerset. For example, given the set [1, 2, 3] the program should compute [[1, 2, 3], [1, 2], [1, 3], [2,3], [1], [2], [3], []]. The program is written in the Felix programming language.
fun find_seq[T] (k:size) (l:list[T]) : list[list[T]] =
  {
    return
      if k > len l
      then
        list[list[T]] () //There are no subsets of length k.
      elif k == 0uz
      then
        list[list[T]] (list[T] ()) //The empty set (the one subset of length zero).
      else
        match l with
          | Cons(?x, ?xs) => 
             join 
              (map (fun (l:list[T]):list[T]=>
                           (join (list[T] x) l)) (find_seq (k - 1) xs))
              (find_seq k xs)
        endmatch
      endif
    ;
  }
  
  fun power_set[T] (lst:list[T]):list[list[T]] =
  {
    fun loop[T] (k:size) (lst:list[T]) (acc:list[list[T]]):list[list[T]] =
    {
      return
        if k == size(-1) then acc
        else loop (k - 1) lst (join acc (find_seq k lst))
        endif
      ;  
    }
  
   return loop (len lst) lst (list[list[T]] ());
  }
  
  println$ str (power_set (list (1, 2, 3)));

Sunday, June 2, 2013

Taylor Polynomials

I find this little calculation a pleasing, succinct reminder of some the (many) things I like about functional programming. It concerns a formulation of the trigonometric sin function as the first few terms of its Taylor polynomial about x = 0:
The "definition" above can be minimally and naturally "coded up" from "first principles" (I mean, using nothing more than elementary arithmetic operations) as simply as this!:

 let rec power  = 
  fun u -> fun i ->
    if i < 0 then 1.0 / (power u (-i))
    else if i = 0 then 1.0 else u * (power u (i-1))

let rec fac = fun n -> if n = 0 then 1 else n * fac (n-1)
let pow = fun i -> fun u -> power u (2 * i + 1)
let coef = fun i -> (power (-1) i)/(fac (2 * i + 1))

let rec taylor = 
  fun i -> fun n -> fun u ->
    if i = n then (coef i) * (pow i u)
    else ((coef i) * (pow i u)) + (taylor (i + 1) n u)
  
let x = 0.78539816339744830961566084581988 (*pi/4*) in 

taylor 0 3 x 
We expect 0.70710678118654752440084436210485. I can get 0.70710646957517809, with the above which is correct to 6 decimal places. I came across this little example in Kluge's "Abstract Computing Machines" by the way, a title I'm happy to recommend to those who share my interest in symbolic computation!

Saturday, June 1, 2013

Building an OCaml package

Building an OCaml package

The details relating to OCaml package construction are provided in the the OCaml System Release 4.00 documentation and users manual. Here's how I go about it. The key to building an OCaml package is to first compile the package modules with the -for-pack flag and then combine the compiled modules with the -pack flag. The procedure can be nicely wrapped up in an OMake macro, presented in this OMakefile schematic. 

make_package(P, libname, S) =
  for_pack (P, f) = 
    OCAMLFLAGS = -for-pack $(P)
    $(f).cmx $(f)$(EXT_OBJ):
    $(f).cmo:
    $(f).cmi:
  foreach (f, $(S))
    for_pack ($(P), $(f))
  section
    .DEFAULT: $(OCamlPackage $(P), $(S))
  .DEFAULT: $(OCamlLibraryInstall _,$(BIN_DIR), $(libname), $(P))

S=$(rootname $(file $(basename $(glob *.ml)) ))

make_package(Foo, foo, $(S))