(* 
%use "Distributions/stencil.fsh";;
*)

%use "Distributions/block.fsh" ;;

(*

Stencilling produces an array of neighbourhoods. Inner stencilling
shrinks the shape of the outer array to ensure that the neighbourhoods
do not trigger array bounds. For example,

>-|> stencil_inner_sh {1,1:int_shape}  {1,1:int_shape} {9,9:int_shape};;
J : #[[int]]
#J = { ~7,~7 : ~3,~3 : int_shape }

Inner stencilling is also used for broadcasting, as follows. If the
proclaimed "mask" is bigger than the array in a given dimension then
the actual mask size will be that of the original array, with all
entries broadcast.

>-|> stencil_inner_sh {10,10:int_shape}  {1,1:int_shape} {9,9:int_shape};;
J : #[[int]]
#J = { ~9,~9 : ~9,~9 : int_shape }

These can be mixed, as follows: 

>-|> stencil_inner_sh {1,10:int_shape}  {1,1:int_shape} {9,9:int_shape};;
J : #[[int]]
#J = { ~7,~9 : ~3,~9 : int_shape }

Outer stencilling does not change the shape of the outer array, and so
must choose values for the location sbeyond the boundaries. One
alternative would be to use null values, but these are not yet part of
the language. Another is to use a function to create new boundary
values as a function of the existing boundary. The idea is to work out
from the innermost dimensions. For a matrix the diagrammatic
description is as follows.

	-------------------------
	|			|
	-------------------------
	|			|
	-------------------------
	| | |		      | |
	|-|-|-----------------|-|
	| | |		      | |
	|-|-|-----Matrix -----|-|
	| | |		      | |
	|-|-|-----------------|-|
	| | |		      | |
	-------------------------
	|			|
	-------------------------

This clumsy diagram indicates that first the rows (as innermost
dimension) are extended, and then the columns. Note that at each stage
new entries are only neighbours to a single entry of the earlier
matrix. Thus a function f:a -> a suffices for an array of type [a].

An example is 

>-|> stencil_outer (fun x -> 0)
        {1,1:int_shape}
        {1,1:int_shape}
        (fill {4,4:int_shape} with 
        [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
;;
J : [[int]]
#J = { ~4,~4 : ~3,~3 : int_shape }

The programs are all poly-dimensional. 
*)

(* Some auxiliary functions *)


(* wides : #[a] -> #[b] -> [int]

produces an array which marks (by 0) the dimensions in which the
mask shape (the first argument) is wider than the array.
*)

let wides mask_sh sh =   
  new #ws = {numdim sh : int_shape} 
  and ndx = 0 in 
  let f0 m sh1 = skip 
  and f1 (n:size) h sh1 = (
	ws[ndx] := 
	  (if n <= (lendim sh1)
	   then 1 
	   else 0) ;
	incr ndx ;
	h (preddim sh1) )
  in 
  folddim f0 f1 mask_sh sh
  return ws 
;;

(* the inner stencil *)

let stencil_inner_sh =
folddim f0 f1 
  where f0 left right sh = zerodim sh 
  and f1 (n:size) h right sh = 
    let p = n + (lendim right) 
    and q = lendim sh 
    and base = h (preddim right) (preddim sh)
    in 
      if p < q 
      then succdim (q-p)
           (extendShape base (succdim (p+1) (zeroShape base)))
      else succdim q
           (extendShape base (succdim q (zeroShape base)))
;;


let stencil_inner_pr (target:var [a]) (source:var a) = 
new ndx0 = wides (zeroShape #target) #source 
and #ndx = #ndx0 in 
let h (ndx1:var [int]) (ndx2:var[int]) =
  ndx := zipop plus_int ndx2 (				// the inner index
	 zipop times_int ndx1 ndx0			// the outer index
	 ) in 
let pr0 ndx1 ndx2 y = (
      h ndx1 ndx2 ;
      y:= entry ndx source )		// get the entry 
and pr ndx1 = idoall (pr0 ndx1) 
in 
  idoall pr target 
end
;;


let stencil_inner left right = 
  proc2fun (stencil_inner_pr) 
	(stencil_inner_sh left right)
;;

(* add_boundary *)

let skip3 x y z = skip ;;

let update_boundary f left right (target:var a) = 
  let g n h lt rt (x:var b) = 
    let m1 = lendim #x
    and	kl = lendim lt
    and kr = lendim rt
  in 
  for kl <= i < m1-kr do
    h (preddim lt) (preddim rt) sub(x,i)
  done ;				// the inner dimensions 
  for i < kl do
    map_pr f sub(x,kl-i-1) sub(x,kl-i) 	// the left extension
  done ;
  for m1-kr <= i < m1 do 		// the right extension 
    map_pr f sub(x,i) sub(x,i-1)
  done
in 
  primrec g skip3 (numdim #target) left right target
;;

let add_boundary_pr f left right target source = 
  (map_pr_offset left identity target source ;
  update_boundary f left right target)
;;
 

let add_boundary_sh left right sh = 
  add_shapes (add_shapes left right) sh
;;

let add_boundary f left right = 
  proc2fun (add_boundary_pr f left right) (add_boundary_sh left right)
;;


(* the outer stencil *)

let stencil_outer f left right = 
  (stencil_inner left right) . (add_boundary f left right)
;;

(*

let st_test = 
  stencil_inner 
	{1,1:int_shape}
	{1,1:int_shape}
	(fill {4,4:int_shape} with 
	[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
;;


let mat = add_boundary (plus_int 1) {1,1:int_shape} {1,1:int_shape} 
		(fill {2,3:int_shape} with [0,1,2,3,4,5])
;;

let three_d = add_boundary  (plus_int 1){1,1,1:int_shape} {1,1,1:int_shape} 
	(fill {2,3,2:int_shape} with [0,1,2,3,4,5,0,1,2,3,4,5])
;;

let sto2 = 
  stencil_outer (fun x -> 0)
	{1,1:int_shape}
	{1,1:int_shape}
	(fill {4,4:int_shape} with 
	[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
;;

let 


*)