(* 
%use "Difference_equations/difference.fsh";;
*)

(* Difference equations 

diff_solver_inner : ([a] -> a) ->
                #[b] ->
                  #[c] ->
                    (var bool -> var [a] -> var [a] -> bool) -> [a] -> [a]

diff_solver_inner eqn left right converged initial acts as follows.  

left and right specify the offsets from the central point of a mask.
eqn is used to produce new entries from the associated masks. 
converged determines whether the new array is "the same as" the old array 
	(the boolean argument is usually ignored).
initial is the initial value of the array.

diff_solver_inner uses inner stencils, and so does not change the
boundary values during computation.

diff_solver_outer uses an outer stencil to supply the necessary
values.

See life.fsh for  examples. 

*)

%use "Distributions/stencil.fsh" ;; 


let diff_solver_inner_pr
	diff_eqn 	// acts on masks to produce new values 
	left 		// the mask shape (size to left of centre) 
	right 		// the mask shape (size to right of centre)
	near 		// intermediate output/criterion for termination
	final		// the result 
	(initial:var [a]) =	// initial conditions

new 
  initial_latest = true // initial is the latest value
in 
let pr y = (map_offset left diff_eqn y) . 
	(stencil_inner left right)
in 	
  final := initial ; 	// set the boundary of final
  pr final initial ;		
  while not (near !initial_latest !initial !final) do 	
   initial_latest := not initial_latest ;	 
   if initial_latest		
   then pr final initial	
   else pr initial final	
  done 			
end
;;

let diff_solver_inner diff_eqn left right near 	
  = proc2fun 
	(diff_solver_inner_pr diff_eqn left right near)
	identity 
;;


let diff_solver_outer_pr 
	boundary_fn	// produces boundary values 
	diff_eqn 	// acts on masks to produce new values 
	left 		// the mask shape (size to left of centre) 
	right 		// the mask shape (size to right of centre)
	near 		// intermediate output/criterion for termination
	final		// final value
	(initial:var [a]) =	// initial conditions
new 
  initial_latest = true // initial is the latest value
in 
let f = (map diff_eqn) . (stencil_outer boundary_fn left right)
in 
  final := f initial ;		
  while not (near !initial_latest !initial !final) do 	
   initial_latest := not initial_latest ;
   if initial_latest		
   then final := f !initial	
   else initial := f !final  
  done 			
end
;;

let diff_solver_outer boundary_fn diff_eqn left right near = 
  proc2fun (diff_solver_outer_pr boundary_fn diff_eqn left right near) 
	   identity
;;


(*

let initial = map (fun x -> int2float x) 
  (fill {5,5:int_shape} with 
  [ 	0,1,2,3,4,
	5,6,7,8,9,
   	10,11,12,13,14,
	15,16,17,18,19,
   	20,21,22,23,24
  ] )
;;


let max_float x y = if x >=. y then x else y ;;
let max_array_float = reduce max_float 0.0 ;;

let l1_norm x y = max_array_float (zipop f x y) 
  where f x y = fabs (x-.y) 
;;

let mean x = sum_float x /. (int2float (number_of_entries #x)) ;;

let dso epsilon = 
  diff_solver_outer 
	(fun x -> 0.0)
	(fun x -> sqrt (mean x))
	{1,1:int_shape}
 	{1,1:int_shape} 
	(fun b -> fun x -> fun y -> l1_norm x y <. epsilon)
	initial
;;


let dso2 = dso 0.1
;;

let dso3 = dso2;;

*)