(* kernel.sml
 *
 * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
 * All rights reserved.
 *
 * QUESTION: should we 
 *)

structure Kernel : sig

    type coefficient = Rational.rat

  (* polynomial represented as list of coefficients, where ith element is
   * coefficient for x^i.
   *)
    type polynomial = coefficient list

    type kernel

  (* kernel name *)
    val name : kernel -> string
    val toString : kernel -> string

  (* are two kernels the same *)
    val same : kernel * kernel -> bool

  (* hash value *)
    val hash : kernel -> word

  (* kernel support *)
    val support : kernel -> int

  (* representation of i'th derivative of the kernel *)
    val curve : kernel * int -> {
	    isCont : bool,
	    segs : polynomial list	(* piece-wise polynomial that defines *)
					(* the curve over the whole support *)
	  }

    val evaluate : polynomial * int -> Rational.rat

  (* some standard kernels *)
    val tent : kernel		(* linear interpolation *)
    val ctmr : kernel		(* Catmull-Rom interpolation *)
    val bspln3 : kernel		(* cubic bspline reconstruction, doesn't interpolate *)
    val bspln5 : kernel		(* quintic bspline reconstruction, doesn't interpolate *)

  end = struct

    structure R = Rational
    structure A = Array

    val maxDiffLevels = 15		(* support upto 15 levels of differentiation *)

    type coefficient = R.rat

    val zero = R.zero
    val one = R.fromInt 1

  (* polynomial represented as list of coefficients, where ith element is
   * coefficient for x^i.
   *)
    type polynomial = coefficient list

    fun differentiate [] = raise Fail "invalid polynomial"
      | differentiate [_] = [zero]
      | differentiate (_::coeffs) = let
	  fun lp (_, []) = []
	    | lp (i, c::r) = R.*(R.fromInt i, c) :: lp(i+1, r)
	  in
	    lp (1, coeffs)
	  end

  (* evaluate a polynomial at an integer coordinate (used to test continuity) *)
    fun evaluate (poly, x) = let
	  val x = R.fromInt x
	  fun eval (sum, [], xn) = sum
	    | eval (sum, c::r, xn) = eval(R.+(sum, R.*(c, xn)), r, R.*(x, xn))
	  in
	    eval (zero, poly, one)
	  end

    type curve = {
	isCont : bool,
	segs : polynomial list	(* piece-wise polynomial that defines *)
				(* the curve over the whole support *)
      }

    datatype kernel = K of {
	name : string,
	id : Stamp.stamp,	(* unique ID *)
	support : int,		(* number of samples to left/right *)
	curves : curve option array (* cache of curves indexed by differentiation level *)
      }

  (* determine if a list of polynomials represents a continuous piece-wise polynomial *)
    fun isContinuous polys = let
	  val s = List.length polys div 2
	  fun chk (i, y_i, []) = (R.zero = y_i)
	    | chk (i, y_i, f_i::r) = let
		val y = evaluate(f_i, i)
		in
		  (y_i = y) andalso chk(i+1, evaluate(f_i, i+1), r)
		end
	  in
	    chk (~s, R.zero, polys)
	  end

  (* kernel name *)
    fun name (K{name, ...}) = name
    val toString = name

  (* kernel support *)
    fun support (K{support, ...}) = support

    fun hash (K{id, ...}) = Stamp.hash id
    fun same (K{id=a, ...}, K{id=b, ...}) = Stamp.same (a, b)

  (* representation of i'th derivative of the kernel *)
    fun curve (K{curves, ...}, k) = (case A.sub(curves, k)
	   of SOME curve => curve
	    | NONE => let
	      (* compute the (k+1)'th derivative, given the k'th *)
		fun diff (k, {isCont, segs}) = let
		      val segs' = List.map differentiate segs
		      val isCont' = isCont andalso isContinuous segs'
		      in {
			isCont = isCont',
			segs = segs'
		      } end
		fun lp (j, curve) = if (j < k)
		      then (case A.sub(curves, j+1)
			 of NONE => let
			      val curve' = diff(j+1, curve)
			      in
				A.update(curves, j+1, SOME curve');
				lp (j+1, curve')
			      end
			  | SOME curve' => lp(j+1, curve')
			(* end case *))
		      else curve
		in
		  lp (0, valOf(A.sub(curves, 0)))
		end
	  (* end case *))

  (* some standard kernels *)
    local
      val op / = R./
      fun r i = R.fromInt i
      fun mkKernel {name, support, segs} = let
	    val curves = Array.array(maxDiffLevels+1, NONE)
	    val curve0 = {
		    isCont = isContinuous segs,
		    segs = segs
		  }
	    in
	      A.update (curves, 0, SOME curve0);
	      K{name=name, id=Stamp.new(), support=support, curves=curves}
	    end
    in
    val tent : kernel = mkKernel{	(* linear interpolation *)
	    name = "tent",
	    support = 1,
	    segs = [
		[r 1, r 1],		(* -1 .. 0 *)
		[r 1, r ~1]		(*  0 .. 1 *)
	      ]
	  }
    val ctmr : kernel = mkKernel{	(* Catmull-Rom interpolation *)
	    name = "ctmr",
	    support = 2,
	    segs = [
		[r 2, r 4,  5/2,  1/2],	(* -2 .. -1 *)
		[r 1, r 0, ~5/2, ~3/2],	(* -1 .. 0 *)
		[r 1, r 0, ~5/2,  3/2],	(*  0 .. 1 *)
		[r 2, r ~4, 5/2, ~1/2]	(*  1 .. 2 *)
	      ]
	  }
    val bspln3 : kernel = mkKernel{	(* cubic bspline reconstruction, doesn't interpolate *)
	    name = "bspln3",
	    support = 2,
	    segs = [
		[ 4/3, r 2,  r 1,   1/6 ],	(* -2 .. -1 *)
		[ 2/3, r 0,  r ~1, ~1/2 ],	(* -1 .. 0 *)
		[ 2/3, r 0,  r ~1,  1/2 ],	(*  0 .. 1 *)
		[ 4/3, r ~2, r 1,  ~1/6 ]	(*  1 .. 2 *)
	      ]
	  }
    val bspln5 : kernel = mkKernel{	(* quintic bspline reconstruction, doesn't interpolate *)
	    name = "bspln5",
	    support = 3,
	    segs = [
		[ 81/40, 27/8, 9/4,  3/4,  1/8, 1/120 ],	(* -3 .. -2 *)
		[ 17/40, ~5/8, ~7/4, ~5/4, ~3/8, ~1/24 ],	(* -2 .. -1 *)
		[ 11/20, r 0,  ~1/2, r 0,  1/4, 1/12 ],		(* -1 .. 0 *)
		[ 11/20, r 0,  ~1/2, r 0,  1/4, ~1/12 ],	(*  0 .. 1 *)
		[ 17/40, 5/8,  ~7/4, 5/4, ~3/8, 1/24 ],		(*  1 .. 2 *)
		[ 81/40, ~27/8, 9/4, ~3/4, 1/8, ~1/120 ]	(*  2 .. 3 *)
	      ]
	  }
    end

  end