(* 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

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

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

    val evaluate : polynomial * int -> Rational.rat

  (* some standard kernels *)
    val tent : kernel		(* linear interpolation *)
    val ctmr : kernel		(* Catmull-Rom interpolation *)
    val bspl3 : kernel		(* cubic bspline reconstruction, doesn't interpolate *)
    val bspl5 : 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 = {
	isOdd : bool,
	isCont : bool,
	segs : polynomial list	(* piece-wise polynomial that defines *)
				(* the curve over the positive support *)
      }

    datatype kernel = K of {
	name : string,
	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 [_] = true
      | isContinuous (f0::r) = let
	  fun chk (i, f_i, []) = (R.zero = evaluate(f_i, i))
	    | chk (i, f_i, f_i1::r) = let
		val y_i = evaluate(f_i, i)
		val y_i1 = evaluate(f_i1, i)
		in
		  (y_i = y_i1) andalso chk(i+1, f_i1, r)
		end
	  in
	    chk (1, f0, r)
	  end

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

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

  (* 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, {isOdd, isCont, segs}) = let
		      val segs' = List.map differentiate segs
		      val isOdd' = not isOdd
		      val isCont' = if isCont andalso isContinuous segs'
			    then (not isOdd') orelse (evaluate(List.hd segs', 0) = R.zero)
			    else false
		      in {
			isOdd = isOdd',
			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 = {
		    isOdd = false,
		    isCont = isContinuous segs,
		    segs = segs
		  }
	    in
	      A.update (curves, 0, SOME curve0);
	      K{name=name, support=support, curves=curves}
	    end
    in
    val tent : kernel = mkKernel{	(* linear interpolation *)
	    name = "tent",
	    support = 1,
	    segs = [[r 1, r ~1]]
	  }
    val ctmr : kernel = mkKernel{	(* Catmull-Rom interpolation *)
	    name = "ctmr",
	    support = 2,
	    segs = [
		[r 1, r 0, ~5/2, 3/2],
		[r 2, r ~4, 5/2, ~1/2]
	      ]
	  }
    val bspl3 : kernel = mkKernel{	(* cubic bspline reconstruction, doesn't interpolate *)
	    name = "bspl3",
	    support = 2,
	    segs = [
		[ 2/3, r 0, r ~1, 1/2 ],
		[ 4/3, r ~2, r 1, ~1/6 ]
	      ]
	  }
    val bspl5 : kernel = mkKernel{	(* quintic bspline reconstruction, doesn't interpolate *)
	    name = "bspl5",
	    support = 3,
	    segs = [
		[ 11/20, r 0, ~1/2, r 0, 1/4, ~1/12 ],
		[ 17/40, 5/8, ~7/4, 5/4, ~3/8, 1/24 ],
		[ 81/40, ~27/8, 9/4, ~3/4, 1/8, ~1/120 ]
	      ]
	  }
    end

  end