OCaml "data language" wrapper around a few basic float data types and bits from GSL

odl.ml

(*
  Very WIP, assuming it goes anywhere.

  Requires:
  - Batteries
  - Ulib (The stdlib extension, not the Unicode library - see oasis-db)
  - gsl-ocaml

  Compile with: ocamlbuild -use-ocamlfind -pkgs gsl,ulib,batteries odl.cma
*)

open Batteries_uni

module Array = struct
  include Array

  (* Modify [b] in place *)
  let modify2 f a b =
    modifyi (
      fun i xb -> f a.(i) xb
    ) b
end

module Matrix = struct
  include Ulib.Umatrix

  let modifyij f m =
    iterij (
      fun i j x ->
        m.(i).(j) <- f i j x
    ) m
  
  let modify f m =
    modifyij (fun _ _ x -> f x) m

  let modify2 f a b =
    modifyij (
      fun i j xb ->
        f a.(i).(j) xb
    ) b
end

type 'a t =
  | S of float
  | A of float array
  | M of float array array

let i x = S (float_of_int x)
let s x = S x
let l x = A (Array.of_list x)
let a x = A x
let m x = M x
let init e =
  A (Array.of_enum e)
let init2 e =
  M (Array.of_enum (Enum.map Array.of_enum e))

let apply_op op a b =
  match a, b with
  | S s1, S s2 -> S (op s1 s2)
  | A a, S s -> A (Array.map (fun x -> op x s) a)
  | S s, A a -> A (Array.map (fun x -> op s x) a)
  | M m, S s -> M (Matrix.map (fun x -> op x s) m)
  | S s, M m -> M (Matrix.map (fun x -> op s x) m)
  | A a1, A a2 -> A (Array.map2 (fun x y -> op x y) a1 a2)
  | A a, M m -> M (Array.map (fun ma -> Array.map2 (fun x y -> op x y) a ma) m)
  | M m, A a -> M (Array.map (fun ma -> Array.map2 (fun x y -> op x y) ma a) m)
  | M m1, M m2 -> M (Matrix.map2 (fun x y -> op x y) m1 m2)

let apply_one_op op x =
  match x with
  | S s -> S (op s)
  | A a -> A (Array.map (fun x -> op x) a)
  | M m -> M (Matrix.map (fun x -> op x) m)

let modify_op op a b =
  match a, b with
  | S s1, S s2 -> invalid_arg "Scalars can not be modified"
  | A a, S s -> Array.modify (fun x -> op x s) a
  | S s, A a -> Array.modify (fun x -> op s x) a
  | M m, S s -> Matrix.modify (fun x -> op x s) m
  | S s, M m -> Matrix.modify (fun x -> op s x) m
  | A a1, A a2 -> Array.modify2 (fun x y -> op x y) a1 a2
  | A a, M m -> Array.modify (fun ma -> Array.map2 (fun x y -> op x y) a ma) m
  | M m, A a -> Array.modify (fun ma -> Array.map2 (fun x y -> op x y) ma a) m
  | M m1, M m2 -> Matrix.modify2 (fun x y -> op x y) m1 m2

let modify_one_op op x =
  match x with
  | S s -> invalid_arg "Scalars can not be modified"
  | A a -> Array.modify (fun x -> op x) a
  | M m -> Matrix.modify (fun x -> op x) m

let map f x =
  match x with
  | S s -> S (f s)
  | A a -> A (Array.map f a)
  | M m -> M (Matrix.map f m)

let amap f x =
  match x with
  | S s -> invalid_arg "Scalar can not be mapped"
  | A a -> A (f a)
  | M m -> M (Array.map f m)

let amap2 f x y =
  match x, y with
  | S s, _
  | _, S s -> invalid_arg "Scalar can not be mapped"
  | A a1, A a2 -> A (f a1 a2)
  | A a, M m -> M (Array.map (f a) m)
  | M m, A a -> M (Array.map (fun ma -> f ma a) m)
  | M m1, M m2 -> M (Array.map2 f m1 m2)

let smap f x =
  match x with
  | S s -> invalid_arg "Scalar can not be mapped"
  | A a -> S (f a)
  | M m -> A (Array.map f m)

let smap2 f x y =
  match x, y with
  | S s, _
  | _, S s -> invalid_arg "Scalar can not be mapped"
  | A a1, A a2 -> S (f a1 a2)
  | A a, M m -> A (Array.map (f a) m)
  | M m, A a -> A (Array.map (fun ma -> f ma a) m)
  | M m1, M m2 -> A (Array.map2 f m1 m2)

let modify f x =
  match x with
  | S s -> invalid_arg "Scalar can not be modified"
  | A a -> Array.modify f a
  | M m -> Matrix.modify f m

let apply_rev x =
  match x with
  | S s -> invalid_arg "Scalar can not be reversed"
  | A a -> A (Array.rev a)
  | M m -> M (Array.rev m)

let modify_rev x =
  match x with
  | S s -> invalid_arg "Scalar can not be reversed"
  | A a -> Array.rev_in_place a
  | M m -> Array.rev_in_place m

module Infix = struct
  module C = struct
    (* Constructive operators *)
    let ( ~- ) x = apply_one_op ( ~-. ) x
    let ( + ) x y = apply_op ( +. ) x y
    let ( - ) x y = apply_op ( -. ) x y
    let ( * ) x y = apply_op ( *. ) x y
    let ( / ) x y = apply_op ( /. ) x y
    let ( ** ) x y = apply_op ( ** ) x y
  end

  module D = struct
    (* Destructive operators *)
    let ( ~-: ) x = modify_one_op ( ~-. ) x
    let ( +: ) x y = modify_op ( +. ) x y
    let ( -: ) x y = modify_op ( -. ) x y
    let ( *: ) x y = modify_op ( *. ) x y
    let ( /: ) x y = modify_op ( /. ) x y
    let ( **: ) x y = modify_op ( ** ) x y
  end
end

module Stats = struct
  let absdev a = smap Gsl.Stats.absdev a
  let correlation a b = smap2 Gsl.Stats.correlation a b
  let covariance a b = smap2 Gsl.Stats.covariance a b
  let kurtosis a = smap Gsl.Stats.kurtosis a
  let mean a = smap Gsl.Stats.mean a
  let min a = smap Gsl.Stats.min a
  let max a = smap Gsl.Stats.max a
  let sd a = smap Gsl.Stats.sd a
  let skew a = smap Gsl.Stats.skew a
  let variance a = smap Gsl.Stats.variance a
end

module Interp = struct
  type interp = Gsl.Interp.interp

  let make kind xs ys =
    match xs, ys with
    | A a, A b -> Gsl.Interp.make_interp kind a b
    | _ -> invalid_arg "Scalars and matrices can't be interpolated"

  let eval_array interp xs =
    let ys = Array.make (Array.length xs) 0.0 in
    Gsl.Interp.eval_array interp xs ys;
    ys

  let eval interp x =
    match x with
    | S s -> S (Gsl.Interp.eval interp s)
    | A a -> A (eval_array interp a)
    | M m -> M (Array.map (eval_array interp) m)
end

include Infix.C
include Infix.D

let ( @@ ) f x = f x

odl.mli

type 'a t = private
  | S of float
  | A of float array
  | M of float array array

val i : int -> [`scalar] t
val s : float -> [`scalar] t
val l : float list -> [`array] t
val a : float array -> [`array] t
val m : float array array -> [`matrix] t
val init : float BatEnum.t -> [`array] t
val init2 : float BatEnum.t BatEnum.t -> [`matrix] t

val apply_op : (float -> float -> float) -> 'a t -> 'b t -> 'c t
val apply_one_op : (float -> float) -> 'a t -> 'b t
val modify_op : (float -> float -> float) -> 'a t -> 'b t -> unit
val modify_one_op : (float -> float) -> 'a t -> unit

val map : (float -> float) -> 'a t -> 'a t

val amap : (float array -> float array) -> 'a t -> 'b t
val amap2 :
  (float array -> float array -> float array) -> 'a t -> 'b t -> 'c t

val smap : (float array -> float) -> 'a t -> 'b t
val smap2 : (float array -> float array -> float) -> 'a t -> 'b t -> 'c t

val modify : (float -> float) -> 'a t -> unit
val apply_rev : 'a t -> 'b t
val modify_rev : 'a t -> unit

module Infix :
  sig
    module C :
      sig
        val ( ~- ) : 'a t -> 'b t
        val ( + ) : 'a t -> 'b t -> 'c t
        val ( - ) : 'a t -> 'b t -> 'c t
        val ( * ) : 'a t -> 'b t -> 'c t
        val ( / ) : 'a t -> 'b t -> 'c t
        val ( ** ) : 'a t -> 'b t -> 'c t
      end
    module D :
      sig
        val ( ~-: ) : 'a t -> unit
        val ( +: ) : 'a t -> 'b t -> unit
        val ( -: ) : 'a t -> 'b t -> unit
        val ( *: ) : 'a t -> 'b t -> unit
        val ( /: ) : 'a t -> 'b t -> unit
        val ( **: ) : 'a t -> 'b t -> unit
      end
  end
module Stats :
  sig
    val absdev : 'a t -> 'b t
    val correlation : 'a t -> 'b t -> 'c t
    val covariance : 'a t -> 'b t -> 'c t
    val kurtosis : 'a t -> 'b t
    val mean : 'a t -> 'b t
    val min : 'a t -> 'b t
    val max : 'a t -> 'b t
    val sd : 'a t -> 'b t
    val skew : 'a t -> 'b t
    val variance : 'a t -> 'b t
  end
module Interp :
  sig
    type interp = Gsl.Interp.interp
    val make : Gsl.Interp.interp_type -> 'a t -> 'b t -> interp
    val eval_array : interp -> float array -> float array
    val eval : interp -> 'a t -> 'b t
  end

include module type of Infix.C
include module type of Infix.D

val ( @@ ) : ('a -> 'b) -> 'a -> 'b