Error propagating calculator

Error propagation is an important tool for measuring the uncertainty of a final result based on the uncertainties of the underlying measurements. It is also the bane of many first year science students' lives. A calculator that worked with uncertain quantities and did the error propagation automatically would be a useful tool, eg 22.0,5.0% + 10.0,1.0 to add 22 with an uncertainty of 5% to 10 with an uncertainty of 1.

Here I show how to accomplish this using Ocaml. Ocaml is a good choice for this task because the program can be broken down easily into its independent sub-tasks using the language facilities:

(Source code for this tutorial is available in [uncert_calc.tgz].)

The second task has already been done in SooHyoung Oh's Ocamlyacc tutorial, at least the hard bits about defining the parser. Fairly simple modifications are required to deal with the new number type.

The lexer, calcLex.mll:

{
  open CalcParse
}

let digit=['0'-'9']
let exponent = ['e' 'E'] ['+' '-'] digit+
let floating = (digit+ '.' digit* | digit* '.' digit+) exponent?

rule token = parse
  | '\n'{  EOL } 
  | floating 
    { (NUMBER (float_of_string (Lexing.lexeme lexbuf)))}
  | ',' {   COMMA }
  | '%' {   PERCENT }
  | '+' {  PLUS }
  | '-' {  MINUS }
  | '*' {  MUL }
  | '/' {  DIV }
  | '(' {  LPAREN }
  | ')' {  RPAREN }
  | _   { token lexbuf }
  | eof { raise End_of_file}

The paser, calcParse.mly:

%{
  let parse_error s = print_endline s; flush stdout
%}

%token EOL
%token LPAREN RPAREN
%token <float> NUMBER
%token COMMA PERCENT
%token PLUS MINUS MUL DIV

%left PLUS MINUS
%left MUL DIV

%start input
%type <unit> input
%%

input: /* empty */ {}
  | input line {} 
;

line:  EOL { }
  | exp EOL 
    { Printf.printf "\t%s\n" (AlgNum.to_string $1);
      flush stdout }
  | error EOL { Printf.printf "Syntax error!\n";
                flush stdout} 
;

exp: NUMBER COMMA NUMBER
      { AlgNum.of_num (Uncert.make $1 $3) }
  | NUMBER COMMA NUMBER PERCENT 
      { AlgNum.of_num (Uncert.make_tol $1 ($3 /. 100.0))}
  | exp PLUS exp  { AlgNum.add $1 $3 }
  | exp MINUS exp { AlgNum.sub $1 $3 }
  | exp MUL exp   { AlgNum.mul $1 $3 }
  | exp DIV exp   { AlgNum.div $1 $3 }
  | LPAREN exp RPAREN { $2 }
;

%%

The driving program, calcmain.ml:

let main () = 
  try 
    let lexbuf = Lexing.from_channel stdin in
    while true do 
      CalcParse.input CalcLex.token lexbuf
    done 
  with End_of_file -> exit 0

let _ = Printexc.print main ()

The first task involves dealing with the intricacies of the Ocaml module, as I have been describing in my previous few blog entries. The basic method is to create two abstract datatypes, one to represent the UncertainNumber type and the other to deal with the AlgebraicNumber datatype for which the addition etc operators are defined. This allows the extension of the calculator to other datatypes such as ComplexNumbers or Vectors by redefining the construction and accessor functions in the underlying representation and the algebraic functions in the AlgebraicNumber datatype. From the point of view of the parser the only functions that it is concerned with is AlgebraicNumber.add etc.

In nums.ml:

(* Type for an Uncertain Number.  This is a number represented by its
   mean and standard deviation. *)
module type UncertSig = sig
  type t
  (* Make an uncertain number from mean and stdev *)
  val make:  float -> float -> t
  (* Make an uncertain number from mean and relative uncertainty in
     mean *)
  val make_tol: float -> float -> t
  (* Extract the mean from an uncertain number *)
  val mean:  t -> float
  (* Extract the stdev (ie uncertainty) from an uncertain number *)
  val stdev: t -> float
  (* Convert the number to a string for printing *)
  val to_string: t -> string
end

module Uncert : UncertSig = struct
  type t = float * float
  let make mean stdev = (mean, stdev)
  let make_tol mean tol = (mean, mean *. tol)
  let mean (m, _) = m
  let stdev(_, s) = s
  let to_string n =
    Printf.sprintf "%f +/- %f" 
      (mean n) (stdev n)
end

(* Algebraic number type.  An algebraic number can be scaled by a
   float, or added, subtracted, multiplied or divided by another
   algebraic number *)
module type AlgSig = sig
  type elt
  type t
  (* The of_num and to_num functions are the interface between this
  abstraction layer and the underlying layer *)
  val of_num: elt -> t
  val to_num: t -> elt
  (* scale number by a float *)
  val scale:  t -> float -> t
  (* algebraic operations on the type *)
  val add: t -> t -> t
  val sub: t -> t -> t
  val mul: t -> t -> t
  val div: t -> t -> t
  (* print number *)
  val to_string: t -> string
end

module AlgNum : AlgSig 
  with type elt = Uncert.t = struct
  type elt = Uncert.t
  type t = float * float
  let of_num uncert = 
    (Uncert.mean uncert, Uncert.stdev uncert)
  let to_num (m, s) =
    Uncert.make m s
  let rt_sum_sqr x y = sqrt((x *. x) +. (y *. y))
  let scale (m, st) sc = (m *. sc, st *. sc)
  let add (m1, s1) (m2, s2) =
    (m1 +. m2, (rt_sum_sqr s1 s2))
  let sub (m1, s1) (m2, s2) = 
    (m1 -. m2, (rt_sum_sqr s1 s2))
  let mul (m1, s1) (m2, s2) =
    scale 
      (1.0, rt_sum_sqr (s1 /. m1) (s2 /. m2)) 
      (m1 *. m2)
  let div (m1, s1) (m2, s2) =
    scale 
      (1.0, rt_sum_sqr (s1 /. m1) (s2 /. m2)) 
      (m1 /. m2)
  let to_string n = 
    Uncert.to_string (to_num n)
end

The third task is to tie the newly defined abstract datatype to the parser. This bit is a tad inelegant as I cannot just use the nums.ml file directly. Instead I have to extract the signature definitions into separate Uncert.mli and AlgNum.mli files (without the module and end lines) and the implementations (ie the struct definitions) into Uncert.ml and AlgNum.ml (changing the elt type of AlgNum to Uncert.t).

Finally, write the Makefile:

all: Uncert.cmo AlgNum.cmo calcParse.cmo calcLex.cmo ucalc
    ocamlmktop Uncert.cmo AlgNum.cmo calcParse.cmo calcLex.cmo -o calc.top

clean: 
    -rm *.cm? calcParse.mli calcParse.ml calcLex.ml

ucalc: Uncert.cmo AlgNum.cmo calcmain.ml calcParse.cmo calcLex.cmo
    ocamlc Uncert.cmo AlgNum.cmo calcParse.cmo calcLex.cmo calcmain.ml -o ucalc

calcLex.cmo: calcLex.ml
    ocamlc -c calcLex.ml

calcParse.cmo: calcParse.ml Uncert.cmo AlgNum.cmo calcParse.cmi
    ocamlc -c calcParse.ml

calcParse.cmi: calcParse.mli
    ocamlc -c calcParse.mli

calcParse.mli calcParse.ml: calcParse.mly
    ocamlyacc calcParse.mly

calcLex.ml: calcLex.mll calcParse.cmi
    ocamllex calcLex.mll

AlgNum.cmo: AlgNum.ml AlgNum.cmi Uncert.cmo
    ocamlc -c AlgNum.ml

Uncert.cmo: Uncert.ml Uncert.cmi
    ocamlc -c Uncert.ml

AlgNum.cmi: AlgNum.mli
    ocamlc -c AlgNum.mli

Uncert.cmi: Uncert.mli
    ocamlc -c Uncert.mli