(* port of recurse.cc *)

(*
This is the code that verifies a disjunct of nonlinear inequalities.
The are given as a list   (tf:tfunction list).  If tf = [f1;....;fk], then
the list represents the inequality (f1 < 0 \/ f2 < 0 .... fk < 0).

The end user should only need to define a cell option,
and then call recursive_verifier, which recursively bisects the domain
until a partition of the domain is found on which verifier_cell gives
a pass on each piece of the partition.

*)

flyspeck_needs "../port_interval/types.hl";;
flyspeck_needs "../port_interval/report.hl";;
flyspeck_needs "../port_interval/interval.hl";;
flyspeck_needs "../port_interval/univariate.hl";;
flyspeck_needs "../port_interval/line_interval.hl";;
flyspeck_needs "../port_interval/taylor.hl";;

module Recurse = struct

open Interval;;
open Univariate;;
open Line_interval;;
open Taylor;;

type cellOption = {
  only_check_deriv1_negative : bool;
  is_using_dihmax : bool;
  is_using_bigface126 : bool;
  width_cutoff : float;
  allow_sharp : bool;
  allow_derivatives : bool;
  mutable iteration_count : int;
  iteration_limit : int;
  recursion_depth : int;
};;

(* cell verification is complex, and we use exceptions to
    exit as soon as the status has been determined.   *)

type  cell_status = 
  | Cell_pass 
  | Cell_counterexample 
  | Cell_inconclusive of (float list * float list * float list * float list * tfunction list);;

exception Return of  cell_status;;

let return c = raise (Return c);;


(* error checking and reporting functions *)

let string_of_domain x =
  let n = mth in
  Printf.sprintf "{%f, %f, %f, %f, %f, %f}" (n x 0) (n x 1) (n x 2) (n x 3) (n x 4) (n x 5);;

let string3 (x,z,s) =  (string_of_domain x ^"\n"^ string_of_domain z ^ "\n" ^ s);;

let boolify _ = true;;

let report_current = boolify o Report.report_timed o string3;;

let report_error = boolify o Report.report_error o string3;;

let report_fatal = boolify o Report.report_fatal o string3;;

(* let t = [0.1;0.2;0.3;0.4;0.5;0.6] in report_error (t,t,"ok");; *)

let periodic_count = 
  let end_count = ref 0 in
    fun () ->
      let _ = end_count := !end_count + 1 in
	(0 = ( !end_count mod 80000));;

let check_limit opt depth = 
  let _ = opt.iteration_count <- opt.iteration_count + 1 in
   ( opt.iteration_count < opt.iteration_limit or opt.iteration_limit = 0 ) &&
     (depth < opt.recursion_depth);;

let sgn x = if (x.lo > 0.0) then 1 else if (x.hi < 0.0) then -1 else 0;;

let rec same_sgn x y = (x = []) or (sgn (hd x) = sgn (hd y) && same_sgn (tl x) (tl y));;


(* a bit tricky because unstable exceptions are common early on,
    and evaluations are very expensive. 
   We don't want a single unstable disjunct to ruin everything.

   When boxes are small, then we throw away unstable disjuncts and continue w/o them.
   When the boxes are still big, we return inconclusive to force a bisection.

   Starting with this procedure, we can throw an exception to return early,
   as soon as we are able to determine the cell_status.  All these exceptions
   get caught at the last line of verify_cell.
*)

let rec set_targets (x,z,x0,z0,tf,tis,opt,maxwidth,has_unstable) =
    try (
      if (tf = []) then
	let _ = not(has_unstable) or return (Cell_inconclusive (x,z,x0,z0,map snd tis)) in
	  List.rev tis
      else (
	let target = evalf (hd tf)  x z in

	let _ =  not( opt.only_check_deriv1_negative) or return
	  (if upper_partial target 0 < 0.0 then   Cell_pass
	   else if  lower_partial target 0 > 0.0 then Cell_counterexample
	   else Cell_inconclusive(x,z,x0,z0,tf)) in
	  
	let _ =  upper_bound target >= 0.0 or return Cell_pass in
	  
	  set_targets(x,z,x0,z0,tl tf,(target,hd tf)::tis,opt,maxwidth,has_unstable));	
    )
    with Unstable -> (
      if (2.0 *. maxwidth > opt.width_cutoff) 
      then set_targets(x,z,x0,z0,tl tf,tis,opt,maxwidth,true)  (* proclaim unstable *)
      else set_targets(x,z,x0,z0,tl tf,tis,opt,maxwidth,has_unstable) (* drop silently *);
    );
;;

let rec delete_false acc tis =
  if (tis=[]) then List.rev acc 
  else if (lower_bound (fst (hd tis)) > 0.0) then delete_false acc (tl tis) 
  else delete_false (hd tis::acc) (tl tis);;

(* If the function is monotone, then we can push the calculation to the boundary.
    There is a theorem that justifies the return Cell_pass steps. Here's a sketch.
     We are doing a verification on a large box x0 z0 and 
   x z represents on of many cells in the box.
   Suppose that there is a counterexample in the large box x0 z0.
   The set of c/e is a compact set.
   Write the inequalities as f1 < 0 \/ ... \/ fk < 0.
   So at a c/e all f1 >= 0 && ... && fk >=0.
   Pick out a particular counterexample by maximizing f1+...+fk over the
   set of c/e's, then among this set pick a c/e with largest y1 coordinate,
   then lexicographically continuing until the largest y6 coordinate.
   This fixes a particular point y in the domain x0 z0.

   I claim that has_monotone will detect y, so that if the inequality is false,
   we find out about it.

   In the partition of x0 z0 (of which x z is one rectangle), pick the rectangle
   containing y
   whose lower endpoint x1 is as large as possible, then continue lexicographically
   through to x6.

   I claim that has_monotone will detect y, when it comes to this box x z.
   This is clear by construction.  If we have monotonic increasing, then y is at the
   right edge of the box.  By construction the box x z is the rightmost containing y
   so the box x z meets the right edge of the large containing box x0 z0.
   Thus, it is enough to recursively search over the right edge of the box. Etc.

 *)

let rec has_monotone tis x z x0 z0 is found  = match is with
  | [] -> (x,z,x0,z0,found)
  | j::js when (mth x j >= mth z j) ->
      has_monotone tis x z x0 z0 js found
  | j::js -> 
    let allpos = List.fold_left (fun a ti -> a && lower_partial (fst ti) j >= 0.0) true tis in
    let allneg = List.fold_left (fun a ti -> a && upper_partial (fst ti) j < 0.0) true tis in
    if (allpos) then
      let _ = (mth z j >= mth z0 j) or return Cell_pass in
      let setj u = table (fun i -> (if i=j then mth z j else mth u i))  in
	has_monotone tis  (setj x) (setj z) (setj x0) (setj z0) js true
    else if (allneg) then
      let _ = (mth x j <= mth x0 j) or  return Cell_pass in
      let setj u = table (fun i -> (if i=j then mth x j else mth u i))  in
	has_monotone tis  (setj x) (setj z) (setj x0) (setj z0) js true
    else has_monotone tis x z x0 z0 js found;;

(* loop as long as monotonicity keeps making progress.  *)

let rec going_strong(x,z,x0,z0,tf,opt) =
  let (y,w) = center_form (x,z) in
  let maxwidth = maxl w in
  let tis = set_targets(x,z,x0,z0,tf,[],opt,maxwidth,false) in
  let epsilon_width = 1.0e-8 in
  let _ = (maxwidth >= epsilon_width) or return
	(if (opt.allow_sharp) then (Report.inc_corner_count(); Cell_pass) 
				   else Cell_counterexample) in
  let tis = delete_false [] tis in
  let _ = (List.length tis > 0) or return Cell_counterexample in
  let (x0,z0) = if (List.length tis < List.length tf) then (x,z)  else (x0,z0) in
  let (x,z,x0,z0,strong) = 
    if (opt.allow_derivatives) then has_monotone tis  x z x0 z0 iter6 false 
    else (x,z,x0,z0,false) in
    if (strong) then going_strong(x,z,x0,z0,map snd tis,opt) else (x,z,x0,z0,maxwidth,tis);;

let guess_optimal_corners (x,z,ti,tf) = 
  let mixedsign = List.fold_left (fun a i -> a or (sgn(mth ti.l.df i)=0)) false iter6 in
  let yyn = table (fun i-> mth (if sgn (mth  ti.l.df i) >0 then x else z) i)  in
  let yyu = table (fun i-> mth (if sgn (mth  ti.l.df i) >0 then z else x) i)  in
  let cn = line_estimate tf yyn in
  let cu = line_estimate tf yyu in
    (mixedsign,yyn,yyu,cn,cu);;

let rec drop_numerically_false acc (x,z,x0,z0,tis) =
    match tis with 
      | [] -> (x0,z0,List.rev tis)
      | (ti,tf) :: tiss -> 
	  if (ti.l.f.lo <= 0.0) then drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss)
	  else 
	    let (mixedsign,yyn,yyu,cn,cu)= guess_optimal_corners(x,z,ti,tf) in
	      if (mixedsign) then drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss)
	      else if (min cn.f.lo cu.f.lo > 0.0 && 
			 same_sgn ti.l.df cn.df && same_sgn ti.l.df cu.df) 
	      then drop_numerically_false acc (x,z,x,z,tiss)
	      else drop_numerically_false ((ti,tf)::acc) (x,z,x0,z0,tiss);;

let rec keep_numerically_true best (x,z,x0,z0,tis,evalue) = match tis with
  | [] -> 
       if (evalue < 0.0 && List.length best >0) then (x,z,best)
       else (x0,z0,tis)
  | (ti,tf)::tiss ->
      if (ti.l.f.hi > 0.0) then keep_numerically_true best (x,z,x0,z0,tiss,evalue)
      else
	let (mixedsign,yyn,yyu,cn,cu) = guess_optimal_corners(x,z,ti,tf) in
	let evalue' = max cn.f.hi cu.f.hi  in
	  if ((evalue' < evalue) &&
		same_sgn ti.l.df cn.df && same_sgn ti.l.df cu.df && not(mixedsign)) then
	    keep_numerically_true [(ti,tf)] (x,z,x0,z0,tiss,evalue')
	  else keep_numerically_true best (x,z,x0,z0,tiss,evalue);;    

(*
This procedure is mostly guided by heuristics that don't require formal
verification. In particular, no justification is required for tossing out inequalities
(since they appear as disjuncts, we can choose which one to prove).

Formal verification is required whenever a Cell_passes is issued,
and whenever the domain (x,z) is restricted.

The record (x0,z0) of the current outer boundary must be restricted to (x,z)
whenever an inequality is tossed out.
*)

let rec verify_cell (x,z,x0,z0,tf,opt) =
  try (
  let _ = not(periodic_count ()) or report_current (x,z,"periodic report") in
  let _ = not(opt.only_check_deriv1_negative) or (List.length tf <= 1) or 
	     failwith "verify_cell: incompatible options" in
	(* XX skip the implementation of rad2, delta126, delta135, for now. *)
  let (x,z,x0,z0,maxwidth,tis) = going_strong(x,z,x0,z0,tf,opt) in
  let (x0,z0,tis) =   if (maxwidth < opt.width_cutoff && opt.allow_derivatives)  
  then  drop_numerically_false [] (x,z,x0,z0,tis)   else (x0,z0,tis) in
  let _ = (List.length tis >0) or return Cell_counterexample in
  let (x0,z0,tis) = 
    if ((maxwidth<opt.width_cutoff)&&(List.length tis >1) && opt.allow_derivatives) 
    then keep_numerically_true [] (x,z,x0,z0,tis,0.0) else (x0,z0,tis) in
    Cell_inconclusive (x,z,x0,z0,map snd tis);
  )
  with Return c -> c;;

let rec recursive_verifier (depth,x,z,x0,z0,tf,opt) = 
  let _ = check_limit opt depth or report_fatal(x,z,Printf.sprintf "depth %d" depth) in
    match verify_cell(x,z,x0,z0,tf,opt)  with
      | Cell_counterexample -> false
      | Cell_pass -> true
      | Cell_inconclusive(x,z,x0,z0,tf) ->
	  (let ( ++ ), ( / ) = up(); upadd, updiv in
	   let w2 = List.map2 upsub z x in 
	   let maxwidth2 = maxl w2 in
	   let j_wide =  try( find (fun i -> mth w2 i = maxwidth2) iter6) with
	     | _ -> failwith "recursive_verifier find" in
	   let yj = (mth x j_wide ++  mth z j_wide) / 2.0 in
	   let delta b v =table (fun i-> if (i=j_wide && b) then yj else mth v i) in
           recursive_verifier(depth+1, delta false x ,delta true z ,x0,z0,tf,opt) &&
           recursive_verifier(depth+1, delta true x ,delta false z ,x0,z0,tf,opt));;


 end;;