(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(* Section: Local Fan Main Estimate                                           *)
(* Chapter: Local Fan                                                         *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2012-05-05                                                           *)
(* Date: 2013-03, revised                                                     *)
(* ========================================================================== *)

(* 
Check completeness (informally) 
of case anaysis in proof of Kepler conjecture, Main Estimate.

The purpose is to transform the list init_cs into terminal_cs
in a finite sequence of regulated moves.

This proves the arrow S_init ==> S_term described in the
appendix to the Local Fan chapter.

*)


module Check_completeness = struct

(*

We don't need the full flyspeck build, really just lib.ml.

#directory "/Users/thomashales/Desktop/googlecode/hol_light";;
#use "hol.ml";;

let flyspeck_dir = 
  (try Sys.getenv "FLYSPECK_DIR" with Not_found -> Sys.getcwd());;

loadt (flyspeck_dir^"/../glpk/sphere.ml");;
loadt (flyspeck_dir^"/general/lib.hl");;

*)


(*
****************************************
BASICS
****************************************
*)




let false_on_fail b x = try (b x) with _ -> false;;

let fail_on_false b x s = (b x or failwith "s");;

let filter_some xs = 
  mapfilter (function |Some x -> x |None -> failwith "none" ) xs;;

let length = List.length;;

let nub = Lib.uniq;;

let sortuniq = setify;;

(*
let sortuniq xs = uniq (Lib.sort (<) xs);;
*)

let psort u = if fst u <= snd u then u else (snd u,fst u);;

let rec cart a b = 
  match a with
    | [] -> []
    | a::rest -> (map (fun x -> (a,x)) b) @ cart rest b;;

let cart2 a = cart a a;;



(* from parse_ineq *)
let unsplit d f = function
  | (x::xs) ->  List.fold_left (fun s t -> s^d^(f t)) (f x) xs
  | [] -> "";;

(* let join_comma  = unsplit "," (fun x-> x);;*)
let join_semi  = unsplit ";" (fun x-> x);;
 let join_lines  = unsplit "\n" (fun x-> x);;
(* let join_space  = unsplit " " (fun x-> x);; *)

let string_of_list f xs = 
  let ss = map f xs in
   "["^(join_semi ss)^"]";;

let string_of_pair (i,j) = 
    Printf.sprintf "(%d,%d)" i j;;

let string_of_triple (i,j,s) = 
  Printf.sprintf "(%d,%d,%3.2f)" i j s;;

let string_of_f f k = 
  let ks = filter (fun (i,j) -> i<j) (cart2 (0--(k-1))) in
  let g (i,j) = (i,j,f i j) in
    string_of_list string_of_triple (map g ks);;

let arc = Sphere_math.arclength;;  (* in glpk directory *)
let sqrt = Pervasives.sqrt;;
let cos  = Pervasives.cos;;

(* a few constants.  We do tests of exact equality on floats,
   but this should work because we never do arithmetic on floats
   after initialization (except in one spot where a machine_eps
   is thrown in), so there is no source of roundoff error
   after initialization. *)

let zero = 0.0;;
let cstab = 3.01;;
let sqrt8 = Pervasives.sqrt(8.0);;

(* 0.0 insert on djz on 2013-4-20 because when long edge is cstab get 0.11 + 0.1 (cstab-cstab) = 0.11. *)
let djz = 0.11 +. 0.0 *. 0.1 *.(cstab -. sqrt8);; (* ear-value on cstab edge *) 
let target = 1.541;;
let two = 2.0;;
let twoh0 = 2.52;;
let sqrt1553 = sqrt(15.53);;
let arc1553 = arc two two sqrt1553;;
let three = 3.0;;
let four = 4.0;;
let fourh0 = 2.0 *. twoh0;;
let upperbd = 6.0;;  (* 6.0 > 4 * h0, upper bound on diags in BB *)

let d_tame i = 
  if (i <= 3) then zero
  else if (i=4) then 0.206 
  else if (i=5) then 0.4819
  else if (i=6) then 0.712
  else target;;

let tame_table_d r s = 
  if (r + 2 * s <= 3) then 0.0
  else
    let r' = float_of_int r in
    let s' = float_of_int s in
      0.103 *. (2.0 -. s') +. 0.2759 *. (r' +. 2.0 *. s' -. 4.0);;


(*
look only at global minumum points for tau*(s,v) (for s fixed)
and such that tau*(s,v)<=0.
Let index(v) = the number of edges st |vi-vj|=a(i,j). 
and such that index(v) is minimized among global minimum points in Bs.
Call this set MM(s) (index-minimizers).

The semantics that guide us will be counterexample propagation;
if a stable constraint system contains a point, does the
output of the function contain a nonempty stable constraint system?

Formally, we will prove statements of the form for various 
f:ACS -> (ACS) list
|- nonempty_M cs ==> (?cs'. mem cs'(f cs) /\ nonempty_M cs').

Termination of algorithms:
Every slice decreases k.
Every subdivision restricts the B-field through a finite set of choices
  (retaining the M-field)
Every deformation sets a new M-field constraint.

As a matter of programming style.
Procedures that refer to particular cs's should be avoided,
except the list terminal list.  Pass everything else in by argument.
Select from terminal list with filter_terminal

a_cs,b_cs,am_bs,bm_cs are symmetric functions for args i,j < k.
Their values for args k and large are not relevant.
*)

type constraint_system = 
{
  k_cs : int;
  d_cs : float;
  a_cs : int->int ->float;
  b_cs: int->int->float;
  am_cs: int->int->float;
  bm_cs: int->int->float;
  js_cs : (int*int) list;  (* psorted, both < k_cs *)
  lo_cs :int list;
  str_cs : int list;
  history_cs : string list;
};;

(*
****************************************
DEBUG
****************************************
*)


let string_of_cs cs =
  let s = string_of_list string_of_int in
  let k = string_of_int cs.k_cs in
  let d = string_of_float cs.d_cs in
  let a = string_of_f cs.a_cs (cs.k_cs) in
  let b = string_of_f cs.b_cs (cs.k_cs) in
  let am = string_of_f cs.am_cs (cs.k_cs) in
  let bm = string_of_f cs.bm_cs (cs.k_cs) in
  let js = string_of_list string_of_pair cs.js_cs in
  let lo = s cs.lo_cs in
  let str = s cs.str_cs in
    Printf.sprintf 
      "{\n k=%s;\n d=%s;\n a=%s;\n b=%s;\n am=%s;\n bm=%s;\n js=%s;\n lo=%s;\n str=%s;\n}\n"
      k d a b am bm js lo str;;

let pp_print_cs f cs = pp_print_string f (string_of_cs cs);;


(* report(string_of_cs quad_std_cs);; *)
(* #install_printer pp_print_cs;; *)

let ks_cs cs = (0-- (cs.k_cs -1));;

let ks_cart cs = cart2 (ks_cs cs);;

let extensional_equality cs cs' = 
  let bk = cs.k_cs = cs'.k_cs in
  let bd = cs.d_cs = cs'.d_cs in
  let m r r' = forall (fun (i,j) -> r i j = r' i j) (ks_cart cs) in
  let ba = m cs.a_cs cs'.a_cs in
  let bb = m cs.b_cs cs'.b_cs in
  let bam = m cs.am_cs cs'.am_cs in
  let bbm = m cs.bm_cs cs'.bm_cs in
  let ps js = map psort js in
  let bj = set_eq (ps cs.js_cs) (ps cs'.js_cs) in
  let bstr = set_eq (cs.str_cs) cs'.str_cs in
  let blo = set_eq cs.lo_cs cs'.lo_cs in
    bk && bd && ba && bb && bam && bbm && bstr && blo && bj;;



let debug_cs = ref [];;


let catch_failure cs_term h = 
  fun cs ->
    try
      h cs
    with Failure _ -> 
      if exists (extensional_equality cs_term) !debug_cs 
      then [cs_term] else [];;

let rec rev_assoc_e e a l =
  match l with
    (x,y)::t -> if e a y then x else rev_assoc_e e a t
  | [] -> failwith "find";;

let rec build_path e dict cs cs' buf = 
  if (e cs cs') then cs::buf else
    let v = rev_assoc_e e cs' dict in
      build_path e dict cs v (cs'::buf);;

build_path (=) [(0,4);(4,9);(9,3);(9,7)] 0 7 [];;

let rec mk_assoc h cs_term dict cs_init = 
  let e = extensional_equality in
  if can (Lib.find  (e cs_term)) cs_init then dict
  else 
    let out = map h cs_init in
    let zs = zip cs_init out in
    let f(z,css) = map (fun cs-> (z,cs)) css in
    let zss = List.flatten (map f zs) in
      mk_assoc h cs_term (zss @ dict) (List.flatten out);;

let debug_trace h cs_init cs_term =
  let dict = mk_assoc h cs_term [] [cs_init] in
    build_path extensional_equality dict cs_init cs_term [];;


let report_cs cs = report(string_of_cs cs);;


let failcs (cs,s) = 
  report_cs cs; debug_cs:= [cs]; failwith s;;


(*
****************************************
MORE BASIC OPERATIONS
****************************************
*)



(*
There are B-fields (hard constraints) and M-fields (soft constraints).
The B-fields define the domain of the constraint system s.  
The M-fields partition the minimizers in B_s.

k,a,b, are the primary B-fields
d,j are secondary B-fields that define the function tau.

lo,str,am,bm are M-fields.

js represented as pairs (i,j) with i<j and adjacent.
a,b are bounds defining BB_s.

am,bm do not affect the stable constraint system s.
The deformation lemmas are all partitioning the M-fields.

Subdivision partitions the B-fields.
Transfer resets the M-fields and changes the B-field.
*)


let modify_cs cs =
  fun ?k:(k=cs.k_cs) ?a:(a=cs.a_cs) ?b:(b=cs.b_cs) 
    ?am:(am=cs.am_cs) ?bm:(bm=cs.bm_cs)
    ?d:(d=cs.d_cs)
      ?js:(js=cs.js_cs) ?lo:(lo=cs.lo_cs) 
      ?str:(str=cs.str_cs) ?h:(history=cs.history_cs) () ->
 {
  k_cs = k;
  a_cs = a;
  b_cs = b;
  am_cs = am;
  bm_cs = bm;
  d_cs = d;
  js_cs = js;
  lo_cs = lo;
  str_cs = str;
  history_cs = history;
};;

let hist cs s = modify_cs cs ~h:(s::cs.history_cs) ();;

(* usage : modify_cs hex_std_cs ~k:4 ();; *)

let reset_to_unadorned cs = 
  modify_cs cs ~am:cs.a_cs ~bm:cs.b_cs ~lo:[] ~str:[] ();;


let ink k i = (i+1) mod k;;

let dek k i = (i+k-1) mod k;;

let inc cs i = ink cs.k_cs i;;

let dec cs i = dek cs.k_cs i;;

let inka k a = 
    map (ink k) a;;

let mk_j cs i = psort (i,inc cs i);;

let memj cs (i,j) = mem (psort (i,j)) cs.js_cs;;

let compact u = length (sortuniq u) = length u;;

let ks_cs cs = (0-- (cs.k_cs -1));;

let ks_cart cs = cart2 (ks_cs cs);;

let alldiag cs =   
  filter (fun (i,j) -> (i<j) && not (inc cs i = j) && not (inc cs j = i))
    (ks_cart cs);;

let allpair cs = 
  filter (fun (i,j) -> (i<j)) (ks_cart cs);;

let alledge cs = 
  let ks = ks_cs cs in
  let ed = map (fun i -> psort(i,inc cs i)) ks in
    ed;;

let all = (fun t -> true);;

let htmax cs p = if (mem p cs.lo_cs) then two else twoh0;;

let htmin = two;;

let arcmin cs (i,j) = arc (htmax cs i) (htmax cs j) (cs.am_cs i j);;

let arcmax cs (i,j) = arc htmin htmin (cs.bm_cs i j);;

let m_count cs = 
  let m i = 
    let j = inc cs i in 
      (cs.b_cs i j > twoh0 or cs.a_cs i j > two) in
    List.length (filter m (ks_cs cs));;


(*
****************************************
BASIC OPS ON CONSTRAINT SYSTEMS
****************************************
*)

(*
let is_stable cs =
  let f (i,j) = 
    (i = j) or (2.0 <= cs.a_cs i j) in
  let ks =  ks_cs cs in 
  let k2s = cart2 ks in
  let bm = forall f k2s or failwith "is_st:bm" in 
  let bs = zero_diag cs.a_cs ks or failwith "is_st:bs" in 
  let g i = (cs.b_cs i (inc cs i) <= cstab) in
  let bs' = forall g ks or failwith "is_st:bs'" in 
  let fj (i,j) = not(memj cs (i,j)) or 
    (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
  let bj = forall fj k2s or failwith "is_st:bj" in 
    bm && bs && bs' && bj;;

let is_tri_stable cs = 
  let bk = (3 = cs.k_cs) or failwith "ts:3" in
  let  f (i,j) = 
    (i = j) or (2.0 <= cs.a_cs i j) in
  let ks =  ks_cs cs in 
  let k2s = cart2 ks in
  let bm = forall f k2s or failwith "ts:bm" in 
  let bs = zero_diag cs.a_cs ks  or failwith "ts:bs" in
  let g i = (cs.b_cs i (inc cs i) < four) in
  let bs' = forall g ks or failwith "ts:bs'" in 
  let fj (i,j) = not(memj cs (i,j)) or 
    (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
  let bj = forall fj k2s or failwith "ts:bj" in 
    bk && bm && bs && bs' && bj;;
*)

let is_unadorned cs = 
  let f (i,j) = (cs.a_cs i j = cs.am_cs i j) && (cs.b_cs i j = cs.bm_cs i j) in
  let k2s = ks_cart cs in
  let bm = forall f k2s (* or failwith "unadorned" *) in
    bm && (cs.lo_cs=[]) && (cs.str_cs=[]);;

let zero_diag a xs = forall (fun i -> ((=) zero) (a i i)) xs;;

let stable_constraint_system_test1 cs = 
  let ks = ks_cs cs in
  let k2s = ks_cart cs in
  let f (i,j) = 
    (cs.a_cs i j = cs.a_cs j i) &&
      (cs.b_cs i j = cs.b_cs j i) &&
      (cs.am_cs i j = cs.am_cs j i) &&
      (cs.bm_cs i j = cs.bm_cs j i) &&
      (cs.a_cs i j <= cs.am_cs i j) &&
      (cs.am_cs i j <= cs.bm_cs i j) &&
      (cs.bm_cs i j <= cs.b_cs i j) in
  let fj (i,j) = (i < j) && ((inc cs i = j) or (inc cs j = i)) in
  let ls =  cs.lo_cs @  cs.str_cs in
  let _ = ((3 <= cs.k_cs) && (cs.k_cs <= 6) ) or failwith ( "test1:k_bound") in
  let _ =  (subset ls ks && subset cs.js_cs k2s) or failwith "test1:I_lo, I_str, J" in
  let _ = (cs.d_cs < 0.9) or failwith "test1:d_bound" in
  let _ = forall f k2s or failwith ("test1:ab_relations") in
  let _ = (compact cs.js_cs) or failwith "test1:J set" in 
  let _ = (forall fj cs.js_cs) or failwith "test1:J edge" in 
  let _ = (length cs.js_cs + cs.k_cs <= 6) or failwith "test1:card_J" in
    ();;  

let stable_constraint_system_test2 cs =
  let ks =  ks_cs cs in 
  let k2s = ks_cart cs in
  let f2 (i,j) = 
    (i = j) or (2.0 <= cs.a_cs i j) in
  let _ = forall f2 k2s or failwith "test2:a packing" in 
  let _ = zero_diag cs.a_cs ks or failwith "test2:a diagonal" in 
  let b_edge i = cs.b_cs i (inc cs i) in
  let b_bound i = if (cs.k_cs > 3) then (b_edge i <= cstab) else (b_edge i < four) in
  let _ = forall b_bound ks or failwith "test2:b edge bounds" in 
  let j_bound (i,j) = not(memj cs (i,j)) or 
    (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
  let _ = forall j_bound k2s or failwith "test2:J bounds" in 
  let edge (i,j) = (j = inc cs i) or (i = inc cs j) in  (* bug corrected 2/23/2013, was (i = inc c i) *)
  let fm (i,j) = edge(i,j) && (i < j) && (two < cs.a_cs i j or twoh0< cs.b_cs i j) in
  let m = length (filter fm k2s) in
  let _ = (m + cs.k_cs <= 6) or failwith "test2:m" in
    ();;

let assert_stable_cs cs = 
  try 
    stable_constraint_system_test1 cs;
    stable_constraint_system_test2 cs
  with Failure s -> failcs(cs,s);;


(* f is an edge preserving map *)

let map_cs f g cs = 
  let ks = ks_cs cs in
  let _ = forall (fun i-> f(g i) = i) ks or failwith "map_cs:not inverse" in
  let _ = forall (fun i-> g(f i) = i) ks or failwith "map_cs:not inverse" in
  let a' a i j = a (f i) (f j) in
  let g2 (i,j) = psort (g i , g j ) in
  {
    k_cs = cs.k_cs;
    d_cs = cs.d_cs;
    a_cs = a' (cs.a_cs);
    b_cs = a' (cs.b_cs);
    am_cs = a' (cs.am_cs);
    bm_cs = a' (cs.bm_cs);
    lo_cs = map g cs.lo_cs;
    str_cs = map g cs.str_cs;
    js_cs = map g2 cs.js_cs;
    history_cs = cs.history_cs;
  };;

let opposite_cs cs = 
  let f i = (cs.k_cs - (i+1)) in
    map_cs f f cs;;

let rotate_cs cs = map_cs (inc cs) (dec cs) cs;;

let rec rotatek_cs k cs = 
  funpow k rotate_cs cs;;

(* anadorned isomorphisms  *)

let unadorned_iso_strict_cs cs cs' = 
  let cs = reset_to_unadorned cs in
  let cs' = reset_to_unadorned cs' in
  let bk = cs.k_cs = cs'.k_cs in
  let bd = cs.d_cs = cs'.d_cs in
  let m r r' = forall (fun (i,j) -> r i j = r' i j) (ks_cart cs) in
  let ba = m cs.a_cs cs'.a_cs in
  let bb = m cs.b_cs cs'.b_cs in
  let ps js = map psort js in
  let bj = set_eq (ps cs.js_cs) (ps cs'.js_cs) in
    bk && bd && ba && bb && bj;;

let unadorned_iso_proper_cs cs cs' = 
  Lib.exists (fun i -> unadorned_iso_strict_cs (rotatek_cs i cs) cs') (ks_cs cs);;

let unadorned_iso_cs cs cs' = 
  (cs.k_cs = cs'.k_cs) && 
 (  unadorned_iso_proper_cs cs cs' or unadorned_iso_proper_cs (opposite_cs cs) cs');;



(*
****************************************
FUNCTION BUILDERS
****************************************
*)

let cs_adj adj diag k i j = 
  let s t = string_of_int t in
  if (i<0) or (j<0) or (i>=k) or (j>=k) 
  then failwith ("adj out of range"^(s k)^(s i)^(s j)) 
  else if (i=j) then zero
  else if (j = ink k i) or (i = ink k j) then adj
  else diag;;

let a_pro pro adj diag k i j = 
  if (i<0) or (j<0) or (i>=k) or (j>=k) then failwith "pro out of range"
  else if (i=j) then zero
  else if (i=0 && j=1) or (j=0 && i=1) then pro
  else if (j = ink k i) or (i = ink k j) then adj
  else diag ;;

(*  OLD BUG IF data not sorted:
let funlist data d i j = 
  if (i=j) then zero
  else assocd (psort (i,j)) data d;;
*)

let funlist data = 
  let data' = map (fun ((i,j),d) -> (psort(i,j),d)) data in
    fun d i j ->
      if (i=j) then zero
      else assocd (psort (i,j)) data' d;;

let override a (p,q,d) i j = 
  if psort (p,q) = psort (i,j) then d else a i j;;

let overrides a data =
  let fd = funlist data in
      fun i j ->
	fd (a i j) i j;;



(*
****************************************
TABLES OF AUGMENTED CONSTRAINT SYSTEMS
****************************************
*)

let mk_unadorned (k,d,a,b,h) = {
  k_cs = k;
  d_cs = d;
  a_cs = a;
  b_cs = b;
  am_cs = a;
  bm_cs = b;
  js_cs = [];
  lo_cs = [];
  str_cs = [];
  history_cs = [h];
};;

let hex_std_cs = mk_unadorned (6, d_tame 6,
  cs_adj two twoh0 6,
  cs_adj twoh0 upperbd 6,"hex std init");;

let pent_std_cs = mk_unadorned (5,d_tame 5,
  cs_adj two twoh0 5,
  cs_adj twoh0 upperbd 5,"pent std init");;

let quad_std_cs = mk_unadorned (4,d_tame 4,
  cs_adj two twoh0 4,
  cs_adj twoh0 upperbd 4,"quad std init");;

let tri_std_cs = mk_unadorned (3,d_tame 3,
  cs_adj two twoh0 3,
  cs_adj twoh0 upperbd 3,"tri std init");;

let pent_diag_cs = mk_unadorned (
   5,
   0.616,
   cs_adj two sqrt8 5,
   cs_adj twoh0 upperbd 5,"pent diag init");;

let quad_diag_cs = mk_unadorned (
   4,
   0.467,
   cs_adj two three 4,
   cs_adj twoh0 upperbd 4,"quad diag init");;

let pent_pro_cs = mk_unadorned (
   5,
   0.616,
   a_pro twoh0 two twoh0 5,
   a_pro sqrt8 twoh0 upperbd 5,"pent pro init");;

let quad_pro_cs = mk_unadorned (
   4,
   0.477,
   a_pro twoh0 two sqrt8 4,
   a_pro sqrt8 twoh0 upperbd 4,"quad pro init");;

let init_cs = [
  hex_std_cs ;
  pent_std_cs ;
  quad_std_cs ;
  tri_std_cs ;
  pent_diag_cs ;
  quad_diag_cs ;
  pent_pro_cs ;
  quad_pro_cs ;
];;

map assert_stable_cs init_cs;;
(forall is_unadorned init_cs) or failwith "init_cs:unadorned";;

(* now for the terminal cases done by interval computer calculation *)

let terminal_hex =
mk_unadorned (
   6,
   d_tame 6,
   cs_adj two cstab 6,
   cs_adj two upperbd 6,"terminal hex");;

let terminal_pent = 
mk_unadorned (
   5,
   0.616,
   cs_adj two cstab 5,
   cs_adj two upperbd 5,"terminal pent");;

(* two cases for the 0.467 bound: all top edges 2 or both diags 3 *)

let terminal_adhoc_quad_5691615370 = (* doesn't exist: rhombus side 2, diags >= 3 *)
mk_unadorned (
   4,
   0.467,
   cs_adj two three 4,
   cs_adj two upperbd 4,"terminal");;

let terminal_adhoc_quad_9563139965B = 
mk_unadorned (
   4,
   0.467,
   cs_adj two three 4,
   cs_adj twoh0 three 4,"terminal");;

let terminal_adhoc_quad_4680581274 = mk_unadorned(
 4,
 0.616 -. 0.11,
 funlist [(0,1),cstab ; (0,2),cstab ; (1,3),cstab ] two,
 funlist [(0,1),cstab ; (0,2),upperbd ; (1,3),upperbd ] two,"terminal 4680581274");;

let terminal_adhoc_quad_7697147739 = mk_unadorned(
 4,
 0.616 -. 0.11,
 funlist [(0,1),sqrt8 ; (0,2),cstab ; (1,3),cstab ] two,
 funlist [(0,1),sqrt8 ; (0,2),upperbd ; (1,3),upperbd ] two,"terminal 7697147739");;

(* special case of tri_492...
let terminal_tri_3456082115 = mk_unadorned(
 3,
 0.5518 /. 2.0,
 funlist [(0,1), cstab; (0,2),twoh0; (1,2),two] two,
 funlist [(0,1), 3.22; (0,2),twoh0; (1,2),two] two,"terminal 3456082115");;
*)

let terminal_tri_7720405539 = mk_unadorned(
 3,
 0.5518 /. 2.0 -. 0.2,
 funlist [(0,1),cstab; (0,2),twoh0; (1,2),two] two,
 funlist [(0,1),3.41; (0,2),twoh0; (1,2),two] two,"terminal 7720405539");;
 
let terminal_tri_2739661360 = mk_unadorned(
 3,
 0.5518 /. 2.0 +. 0.2,
 funlist [(0,1),cstab; (0,2),cstab; (1,2),two] two,
 funlist [(0,1),3.41; (0,2),cstab; (1,2),two] two,"terminal 2739661360");;
 
(* range increased by combining with previous case *)

let terminal_tri_9269152105 = mk_unadorned(
 3,
 0.5518 /. 2.0 ,
 funlist [(0,1),cstab; (0,2),cstab; (1,2),two] two,
 funlist [(0,1),3.62; (0,2),cstab; (1,2),two] two,"terminal 9269152105");;

let terminal_tri_4922521904 = mk_unadorned(
 3,
 0.5518 /. 2.0 ,
 funlist [(0,1),cstab; (0,2),twoh0; (1,2),two] two,
 funlist [(0,1),3.339; (0,2),twoh0; (1,2),two] two,"terminal 4922521904");;

(*
Note: May 25, 2012.
The interval ineq requires bounds on (0,2) of 3.41 3.634,
but we have Delta[3.634,2,2,3.01,2.52,3.01]<0, so the upper bound 
always holds.
In the range 3.01--3.41, we can slice and use
upper echelon inequalities _772 and _273.
I haven't implemented this, but it would be an easy extension of
the upper_echelon procedure.

OK. Implemented as upper_echelonC.
We no longer need quad_163.
We use 5512912661 deformation: 3.15/h0 instead
*)

(*
let terminal_quad_1637868761 = mk_unadorned(
 4,
 0.5518,
 funlist [(0,1),two; (1,2), cstab; 
                  (2,3),twoh0; (0,3),two; (0,2),3.41; (1,3),cstab] two,
 funlist [(0,1),two; (1,2), cstab;
                  (2,3),twoh0; (0,3),two; (0,2),3.634; (1,3),upperbd] two,"terminal 1637868761");;
*)

(*
let terminal_quad_1637868761 = mk_unadorned(
 4,
 0.5518,
 funlist [(0,1),two; (1,2), cstab; 
                  (2,3),twoh0; (0,3),two; (0,2),cstab; (1,3),cstab] two,
 funlist [(0,1),two; (1,2), cstab;
                  (2,3),twoh0; (0,3),two; (0,2),upperbd; (1,3),upperbd] two,"terminal 1637868761");;
*)


let ear_cs = 
{
  k_cs = 3;
  d_cs = 0.11;
  a_cs = funlist [(0,1),sqrt8] two;
  b_cs = funlist [(0,1),cstab] twoh0;
  am_cs = funlist [(0,1),sqrt8] two;
  bm_cs = funlist [(0,1),cstab] twoh0;
  js_cs = [0,1];
  lo_cs = [];
  str_cs = [];
  history_cs= ["ear 3603097872"];
};;

(* two consequences of usual ear *)
let ear_jnull = modify_cs ear_cs ~js:[] ();;

let ear_stab = mk_unadorned(
  3,
  djz,
  funlist [(0,2),cstab] two,
  funlist [(0,2),cstab] twoh0,
  "ear_stab"
);;

let terminal_ear_3603097872 = ear_cs;;

let terminal_tri_5405130650 = 
{
  k_cs = 3;
  d_cs =  0.477 -.  0.11;
  a_cs = funlist [(0,1),sqrt8;(0,2),twoh0;(1,2),two] two;
  b_cs = funlist [(0,1),cstab;(0,2),sqrt8;(1,2),twoh0] twoh0;
  am_cs = funlist [(0,1),sqrt8;(0,2),twoh0;(1,2),two] two;
  bm_cs = funlist [(0,1),cstab;(0,2),sqrt8;(1,2),twoh0] twoh0;
  js_cs = [0,1];
  lo_cs = [];
  str_cs = [];
  history_cs=["terminal 5405130650"];
};;

let terminal_tri_5766053833 = 
{
  k_cs = 3;
  d_cs =  0.696 -. 2.0 *. 0.11; 
  a_cs = funlist [(0,1),sqrt8;(1,2),sqrt8] two;
  b_cs = funlist [(0,1),cstab;(1,2),cstab] two;
  am_cs = funlist [(0,1),sqrt8;(1,2),sqrt8] two;
  bm_cs = funlist [(0,1),cstab;(1,2),cstab] two;
  js_cs = [(0,1);(1,2)];
  lo_cs = [];
  str_cs = [];
  history_cs=["terminal 5766053833"];
};;

let terminal_tri_5026777310a = mk_unadorned(
 3,
  0.6548 -. 2.0 *. 0.11,
 funlist [(0,1),sqrt8;(1,2),sqrt8] two,
 funlist [(0,1),cstab;(1,2),cstab] twoh0,"terminal 5026777310a");;

let terminal_tri_7881254908 = mk_unadorned(
 3,
  0.696 -. 2.0 *. 0.11,
 funlist [(0,1),sqrt8;(1,2),sqrt8] twoh0,
 funlist [(0,1),cstab;(1,2),cstab] twoh0,"terminal 7881254908");;

let terminal_std_tri_OMKYNLT_3336871894 = mk_unadorned(
 3,
 zero,
 funlist [] two,
 funlist [] two,"terminal 3336871894");;


(*
(* 1107929058 *)
let terminal_std_tri_OMKYNLT_2_1  = mk_unadorned(
 3,
  tame_table_d 2 1,
 funlist [(0,1),twoh0] two,
 funlist [(0,1),twoh0] two,"terminal 1107929058");;

let terminal_std_tri_7645170609 = mk_unadorned(
 3,
  tame_table_d 2 1,
 funlist [(0,1),sqrt8] two,
 funlist [(0,1),sqrt8] two,"terminal 7645170609");;

(* 1532755966 *)
let terminal_std_tri_OMKYNLT_1_2  = mk_unadorned(
 3,
  tame_table_d 1 2,
 funlist [(0,1),two] twoh0,
 funlist [(0,1),two] twoh0,"terminal 1532755966");;

let terminal_std_tri_7097350062a = mk_unadorned(
 3,
  tame_table_d 1 2,
 funlist [(0,1),twoh0;(0,2),sqrt8] two,
 funlist [(0,1),twoh0;(0,2),sqrt8] two,"terminal 7097350062a");;

let terminal_std_tri_2900061606 = mk_unadorned(
 3,
  tame_table_d 1 2 +. (tame_table_d 2 1 -. 0.11),
 funlist [(0,1),twoh0;(0,2),cstab] two,
 funlist [(0,1),twoh0;(0,2),cstab] two,"terminal 2900061606");;

let terminal_std_tri_2200527225 = mk_unadorned(
 3,
  tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
 funlist [(0,1),two;] sqrt8,
 funlist [(0,1),two;] sqrt8,"terminal 2200527225");;

let terminal_std_tri_3106201101 = mk_unadorned(
 3,
  tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
 funlist [(0,1),two;(0,2),cstab] sqrt8,
 funlist [(0,1),two;(0,2),cstab] sqrt8,"terminal 3106201101");;

let terminal_std_tri_9816718044 = mk_unadorned(
 3,
  tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
 funlist [(0,1),two] cstab,
 funlist [(0,1),two] cstab,"terminal 9816718044");;

let terminal_std_tri_1080462150 = mk_unadorned(
 3,
  tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
 funlist [] twoh0,
 funlist [] twoh0,"terminal 1080462150");;

let terminal_std_tri_4143829594 = mk_unadorned(
 3,
  tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
 funlist [(0,1),cstab] twoh0,
 funlist [(0,1),cstab] twoh0,"terminal 4143829594");;

let terminal_std_tri_7459553847 = mk_unadorned(
 3,
  tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
 funlist [(0,1),twoh0] cstab,
 funlist [(0,1),twoh0] cstab,"terminal 7459553847");;

let terminal_std_tri_4528012043 = mk_unadorned(
 3,
  tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
 funlist [] cstab,
 funlist [] cstab,"terminal 4528012043");;
*)

(* added 2013-06-04 *)

let terminal_tri_4010906068 = mk_unadorned(
 3,
  0.476,
 funlist [] twoh0,
 funlist [] cstab,"terminal 4010906068");;

let terminal_tri_6833979866 = mk_unadorned(
 3,
  0.2759,
 funlist [(0,1),two] twoh0,
 funlist [(0,1),twoh0] cstab,"terminal 6833979866");;

let terminal_tri_5541487347 = mk_unadorned(
 3,
  0.103,
 funlist [(0,1),twoh0] two,
 funlist [(0,1),sqrt8] twoh0,"terminal 5541487347");;




(* use unit_cs as a default terminal object

let unit_cs = terminal_std_tri_OMKYNLT_3336871894;;
 *)


let terminal_cs = [
 terminal_hex;
 terminal_pent; 
 terminal_adhoc_quad_5691615370; 
 terminal_adhoc_quad_9563139965B; 
 terminal_adhoc_quad_4680581274; 
 terminal_adhoc_quad_7697147739; 

(* upper echelon related *)
 terminal_tri_7720405539; 
 terminal_tri_2739661360; 
 terminal_tri_9269152105; 
 terminal_tri_4922521904; 

(* ear related *)
 terminal_ear_3603097872; 
 ear_jnull; 
 terminal_tri_5405130650; 


 terminal_tri_5026777310a; 
 terminal_tri_7881254908; 
 terminal_std_tri_OMKYNLT_3336871894; 

(* added 2013-06-04 *)
terminal_tri_4010906068;
terminal_tri_6833979866;
terminal_tri_5541487347;

(* terminal_tri_3456082115; special case of tri_492... *)
(* terminal_quad_1637868761;  *)
(* ear_stab; test 2013-06-03 *)
(* terminal_tri_5766053833; *) (* removed 2013-06-03 *)


(* removed 2013-06-04
 terminal_std_tri_OMKYNLT_2_1 ; (* 1107929058 *)
 terminal_std_tri_7645170609; 
 terminal_std_tri_OMKYNLT_1_2 ; (* 1532755966 *)
 terminal_std_tri_7097350062a; 
 terminal_std_tri_2900061606; 
 terminal_std_tri_2200527225; 
 terminal_std_tri_3106201101; 
 terminal_std_tri_9816718044; 
 terminal_std_tri_1080462150; 
 terminal_std_tri_4143829594;
 terminal_std_tri_7459553847;
 terminal_std_tri_4528012043;
*)

];;

let filter_terminal f = filter f terminal_cs;;

map ( assert_stable_cs) (filter_terminal all);;
(forall ( is_unadorned) (filter_terminal all)) or failwith "terminal_cs:unadorned";;

let is_ear cs = 
  (cs.k_cs = 3) && (length cs.js_cs = 1) && (unadorned_iso_cs cs  ear_cs);;


(*
****************************************
OPERATIONS AND TRANSFORMATIONS
****************************************
*)

(* transfer move an M-field to a larger 
   B-fields resetting M-fields. b5 not used *)

let transfer_to = 
  let machine_eps = 1e-08 in 
  let b1 cs cs' = (cs.k_cs =  cs'.k_cs) in
  let b2 cs cs' = (cs.d_cs <= cs'.d_cs +. machine_eps) in
  let c3 cs cs' (i,j) = 
    cs.am_cs i j >= cs'.a_cs i j  && cs.bm_cs i j <= cs'.b_cs i j  in
  let b3 cs cs' = forall (c3 cs cs') (allpair cs) in 
  let b4 cs cs' = 
    if is_ear cs then is_ear cs'
    else subset (map psort cs'.js_cs) (map psort cs.js_cs) in
  let b5 cs cs' = 
    subset cs'.lo_cs cs.lo_cs && 
      subset cs'.str_cs cs.str_cs in
    fun cs cs' -> 
      is_unadorned cs' && 
      b1 cs cs' && b2 cs cs' && b3 cs cs' && b4 cs cs' ;;

let proper_transfer_cs cs cs' = 
  Lib.exists (fun i -> transfer_to (rotatek_cs i cs) cs') (ks_cs cs);;

let equi_transfer_cs cs cs' = 
  (cs.k_cs = cs'.k_cs) && 
 (  proper_transfer_cs cs cs' or proper_transfer_cs (opposite_cs cs) cs');;

let equi_transfer_to_list terminals = 
  let f1 = map (C equi_transfer_cs) terminals in
    fun cs -> exists (fun f -> f cs) f1;;

let rec transfer_union a b = 
  match a with
      [] -> b
    | a::aas -> if exists (equi_transfer_cs a) b 
      then transfer_union aas b
      else 
	let b' = filter (not o (C equi_transfer_cs a)) b in
	  transfer_union aas (a::b');;

(* optimized version of equi_transfer_to_list *)

let x_equi_transfer_to_list terms = 
  let machine_eps = 1e-08 in
  let _ = not (can (Lib.find (not o is_unadorned)) terms) or failwith "eq:unadorned" in
  let has_ear = exists is_ear terms in
  let rotates cs = map (fun i -> rotatek_cs i cs) (ks_cs cs) in
  let orotates cs = map (fun i -> rotatek_cs i (opposite_cs cs)) (ks_cs cs) in
  let props = List.flatten (map rotates terms) in
  let oprops = List.flatten (map orotates terms) in
  let terms =  (props @ oprops)  in
  let list_le = forall2 (<=) in 
  let eval_ls f cs = map (fun(i,j)-> f i j) (allpair cs) in
  let dump_e cs = (cs.k_cs,cs.d_cs+.machine_eps,		    
		   eval_ls cs.a_cs cs,eval_ls cs.b_cs cs,
		   cs.js_cs) in
  let dump_m cs = (cs.k_cs,cs.d_cs,eval_ls cs.am_cs cs,eval_ls cs.bm_cs cs,
		   sortuniq (map psort cs.js_cs)) in
  let dump' = setify (map (dump_e) terms) in
    fun cs ->
      if is_ear cs then has_ear
      else
	let (k,d,am,bm,js) = dump_m cs in
	  exists (fun (k',d',a',b',js') ->
		    (k=k') && (d<=d') && (list_le a' am) &&
		      (list_le bm b') && (subset (js') (map psort js))) dump';;

(* changes B-field, M remains a minimizer on smaller domain. Type 1 restriction. 
   The only places b gets modified is with restrict_type1_cs, restrict_type2_cs, subdivision, and 
   initialization of pents.  
*)

let restrict_type1_cs cs = 
  (* preconditions added 2013-05-15 *)
  let k = cs.k_cs in
  let jcount = List.length (filter (fun i -> memj cs (i,inc cs i)) (ks_cs cs)) in
  let _ = jcount = 0 or 3 < k or failwith "type1 cs 1" in
    (modify_cs cs ~b:cs.bm_cs ());;

(* restrict shrinks to B-field down to the M-field,
   leaving M-field intact. *)

let restrict_type2_cs cs =
  (* 2013-05-14, added preconditions *)
  let k = cs.k_cs in
  let jcount = List.length (filter (fun i -> memj cs (i,inc cs i)) (ks_cs cs)) in
  let _ = jcount = 0 or failwith "type2 cs" in
  let _ = forall (fun i -> cs.am_cs i i = zero) (ks_cs cs) or failwith "type2 zero" in
  let cs' = modify_cs cs ~a:cs.am_cs ~b:cs.bm_cs () in
  let _ = m_count cs' + k <= 6 or failwith "type2 cs 2" in
    cs';;

(* half_slice works on B-fields, resetting M-fields. 
   This version does not create an ear. 
   It does just one side. *)

let half_slice cs p q dv = 
  (* 2013-05-14, long diag preconditions added *)
  let k = cs.k_cs in
  let _ = (k = 4) or (cs.bm_cs p q <= cstab) or failwith "half_slice long diag" in
  let _ = (cs.bm_cs p q < 4.0) or failwith "half_slice long diag" in
  let p = p mod k in
  let q = q mod k in
  let q' = if (q<p) then q+k else q in
  let k' = 1+q'-p in
  let av = cs.am_cs p q in
  let bv = cs.bm_cs p q in
  let a = override cs.a_cs (p,q,av) in
  let b = override cs.b_cs (p,q,bv) in
  let _ = (k' >2) or failwith "half_slice underflow" in
  let _ = (k'  <  k) or failwith "half_slice overflow" in
  let r i = (i + k - p) mod k in
  let s i' = (i' + p) mod k in
  let shift f i' j' = f (s i') (s j') in
  let cd1 = mk_unadorned (k',dv,shift a, shift b,"slice:"^(join_semi cs.history_cs)) in
  let js = map (fun (i,j) -> psort (r i,r j)) (intersect (cart2 (p--q)) cs.js_cs) in
  let cd2 = modify_cs cd1 ~js:js () in
    cd2;;

(* calculate the value of the new d from the old one,
   using standard assumptions about the desired terminal inequalities.
   This does cases where no "ear" is involved.
   This inequality is to be used when every edge has bounds in
   one of the three intervals [2,2h0], [2h0,sqrt8], [sqrt8,3.01] *)

let calc_d_std_earless cs p q = 
  let k = cs.k_cs in
  let p = p mod k in
  let q = q mod k in
  let d = cs.d_cs in
  let q' = if (q<p) then q+k else q in
  let k' = 1+q'-p in
  let _ = (k'=3) or raise Not_found in
  let av = cs.am_cs p q in
  let bv = cs.bm_cs p q in
  let a = override cs.a_cs (p,q,av) in
  let b = override cs.b_cs (p,q,bv) in
  let edge3(i,j) = (sqrt8 <= a i j && b i j <= cstab)  in
  let edge2(i,j) = (twoh0 <= a i j && b i j <= sqrt8) && (not(edge3(i,j))) in
  let edge1(i,j) = (two <= a i j  && b i j <= twoh0) && (not (edge2 (i,j))) in
  let edge3s(i,j) = (a i j = cstab && b i j = cstab) in
  let edgeL(i,j) = (twoh0 <= a i j && b i j <= cstab) in
  let edge (i, j) = edge1 (i ,j) or edge2 (i, j) or   edge3 (i, j) or edgeL(i,j) in
  let (p1,p2,p3) = (p,inc cs p,funpow 2 (inc cs) p) in
  let triedge = [(p1,p2);(p2,p3);(p3,p1)] in
  let _ = forall edge triedge or failcs (cs, "slice_dstd") in 
  let c_edge1 = length (filter (edge1) triedge) in
  let c_edge2 = length (filter (edge2) triedge) in
  let c_edge3= length (filter (edge3) triedge) in
  let c_edge3s = length (filter edge3s triedge) in
  let c_edgeL = length (filter edgeL triedge) in
  let eps = tame_table_d 2 1 -. 0.11 in 
  let d12 = tame_table_d 1 2 in
  let m3 = eps *. float_of_int (c_edge3) in
  let flat1 = (c_edge1=2) && (c_edge2=1), tame_table_d 2 1 in
  let flat2s = (c_edge1=2) && (c_edge3s=1), djz in
  let flat2 = (c_edge1=2) && (c_edge3=1), 0.11 in
  let atype = (c_edge1=1) && (c_edge2+c_edge3 =2), d12 +.  m3 in 
  let btype = (c_edge1=0) && (c_edge2+c_edge3=3), tame_table_d 0 3 +. 3.0 *. eps in
  let typeL = (c_edge1=1) && (c_edgeL=2),0.2619 in
  let (_,dpq) = List.find (fst) [flat1;flat2s;flat2;atype;btype;typeL] in
    (dpq, d-. dpq);;

let sort_slice_order cs p q = 
  let inc2 = funpow 2 (inc cs) in
  let edge1(i,j) = (two <= cs.a_cs i j  && cs.b_cs i j <= twoh0) in
  if (cs.k_cs > 4) 
  then 
    (if (inc2 p = q) then (p,q,false) else (q,p,true))
  else 
    let p1 = inc cs p in
    let _ = (inc cs p1 = q) or  failwith "sso:not a diag" in
      (if (edge1(p,p1) && edge1(p1,q)) then (p,q,false) else (q,p,true));;
  
let calc_d_std cs p q = 
  let k = cs.k_cs in
  let swap (a,b) = (b,a) in 
  let inc3 = funpow 3 (inc cs) in
    if (k=6) && (inc3 p = q) 
    then 
      let d = cs.d_cs /. two in (d,d)
    else
      let (p,q,swapped) = sort_slice_order cs p q in
      let d = 
	try 
	  calc_d_std_earless cs p q 
	with Not_found  -> failwith "calc_d_std" in
	(if swapped then swap d else d);;

let slice_cs cs p q dvpq dvqp mk_ear = 
  let machine_eps = 1e-08 in  (* will go away in formalization *)
  let _ = dvpq +. dvqp +. machine_eps >= cs.d_cs or failwith "slice_cs:bad d" in
  let cpq = half_slice cs p q dvpq in
  let cqp = half_slice cs q p dvqp in
  let addj cs = modify_cs cs ~js:((0,(cs.k_cs - 1))::cs.js_cs) () in
  let cpq = if (mk_ear) then addj cpq else cpq in
  let cqp = if (mk_ear) then addj cqp else cqp in
  let _ = not(mk_ear) or is_ear cpq or is_ear cqp or failwith "slice_cs:ear" in
    [cpq;cqp];;

let slice_std cs p q =
  let (dvpq,dvqp) = calc_d_std cs p q in
  let mk_ear = false in
  let rv = slice_cs cs p q dvpq dvqp mk_ear in
  let _ = forall (fun cs' -> cs'.k_cs < cs.k_cs) rv in
    rv;;


(*
****************************************
DEFORMATIONS
****************************************
*)


(* deformation acts on M-fields, keeping the domain fixed. *)

let deform_ODXLSTC_cs p cs = 
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  let ks = ks_cs cs in
  let ksp = subtract ks [p] in
  let diag = subtract ks [p0;p1;p2] in
  let _ = mem p ks or failwith "odx:out of range" in
  let _ = not(memj cs (p0,p1)) or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p cs.lo_cs) or raise Unchanged in
  let m q =  (cs.a_cs p1 q = cs.bm_cs p1 q) in
  let _ = forall (not o m) ksp or raise Unchanged in
  let n q = (fourh0 < cs.b_cs p1 q) in
  let _ = forall n diag or raise Unchanged in
  let cs1 = modify_cs cs ~lo:(sortuniq (p::cs.lo_cs)) () in
  let csp q = 
    if (cs.am_cs p q > cs.a_cs p q) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in 
    cs1::(filter_some(map csp ksp));;

let deform_IMJXPHR_cs p cs = 
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  let ks = ks_cs cs in
  let ksp = subtract ks [p] in
  let diag = subtract ks [p0;p1;p2] in
  let _ = mem p ks or failwith ("imj:out of range"^(string_of_int p)) in
  let _ = mem p cs.str_cs or raise Unchanged in 
  let _ = not(memj cs (p0,p1)) or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p cs.lo_cs) or raise Unchanged in
  let m q =  (cs.a_cs p1 q = cs.bm_cs p1 q) in
  let _ = forall (not o m) diag or raise Unchanged in
  let _ = (m p0 <> m p2) or raise Unchanged in  (* boolean xor *)
  let p' = if (m p0) then p0 else p2 in
  let ksp' = subtract ksp [p'] in
  let n q = (fourh0 < cs.b_cs p1 q) in
  let _ = forall n diag or raise Unchanged in
  let cs1 = modify_cs cs ~lo:(sortuniq (p::cs.lo_cs)) () in
  let csp q =
    if (cs.am_cs p q > cs.a_cs p q) then None
    else Some(modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in 
    cs1::(filter_some(map csp ksp'));;

let deform_NUXCOEA_cs p cs =     
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  let ks = ks_cs cs in
  let ksp = subtract ks [p] in
  let diag = subtract ks [p0;p1;p2] in
  let _ = mem p ks or failwith ("nux:out of range"^(string_of_int p)) in
  let _ = mem p cs.str_cs or raise Unchanged in 
  let _ = not(memj cs (p0,p1)) or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let m q =  (cs.a_cs p1 q = cs.bm_cs p1 q) in
  let _ = forall (not o m) diag or raise Unchanged in
  let _ = (m p0 <> m p2) or raise Unchanged in  (* boolean xor *)
  let p' = if (m p0) then p0 else p2 in
  let ksp' = subtract ksp [p'] in
  let n q = (fourh0 < cs.b_cs p1 q) in
  let _ = forall n diag or raise Unchanged in
  let csp q =
    if (cs.am_cs p q > cs.a_cs p q) then None
    else Some(modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in 
    (filter_some(map csp ksp'));;

(* 
apply M-field deformation at (p,p+1)=(p1,p2) 
This assumes not lunar.
*)

let deform_1834976373_A1_single p cs =
  let p1 = p in
  let p2 = inc cs p in
  let ks = ks_cs cs in
  let alldiag' = alldiag cs in 
  let _ = mem p ks or failwith "1834:out of range" in
  let _ = arcmax cs (p1, p2) < arc1553 or
    raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m (i,j) =  (cs.a_cs i j = cs.bm_cs i j) in
  let _ = forall (not o m) alldiag' or raise Unchanged in
  let _ = not (m (p1,p2)) or raise Unchanged in
  let m' (i,j) = (cs.b_cs i j = cs.am_cs i j) in
  let _ = not (m'(p1,p2)) or raise Unchanged in
  let n (i,j) = (fourh0 < cs.b_cs i j) in
  let _ = forall n alldiag' or raise Unchanged in
  let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
  let cspq (p,q) = 
    if (cs.am_cs p q > cs.a_cs p q) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in 
  let cspq' (p,q) = 
    if (cs.bm_cs p q < cs.b_cs p q) then None
    else Some (modify_cs cs ~am:(override cs.am_cs (p,q,cs.b_cs p q)) ()) in 
    cs1::cs2::
      (filter_some(cspq' (p1,p2) :: (cspq (p1,p2)) ::(map cspq alldiag')));;

(* 
apply M-field deformation at straight p=p1, double edge (p0,p1) (p1,p2) 
This assumes not lunar.
*)

let deform_1834976373_A1_double p cs =
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  let ks = ks_cs cs in
  let alldiag = alldiag cs in
  let _ = mem p ks or failwith "1834-double:out of range" in
  let _ = (arcmax cs (p1,p2) +. arcmax cs (p0,p1) < arc1553) or
    raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(memj cs (p0,p1)) or raise Unchanged in
  let _ = not(mem p0 cs.str_cs) or raise Unchanged in
  let _ = mem p1 cs.str_cs or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m (i,j) =  (cs.a_cs i j = cs.bm_cs i j) in
  let _ = forall (not o m) alldiag or raise Unchanged in
  let _ = not (m (p1,p2) && m(p0,p1)) or raise Unchanged in
  let m' (i,j) = (cs.b_cs i j = cs.am_cs i j) in
  let _ = not (m'(p1,p2) && m'(p0,p1)) or raise Unchanged in
  let n (i,j) = (fourh0 < cs.b_cs i j) in
  let _ = forall n alldiag or raise Unchanged in
  let cs1 = modify_cs cs ~str:(sortuniq (p0::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
  let cspq (p,q) = 
    if (cs.am_cs p q > cs.a_cs p q) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in 
  let csmin =
    if (cs.am_cs p1 p2 > cs.a_cs p1 p2 or cs.am_cs p0 p1 > cs.a_cs p0 p1)
    then None
    else Some (modify_cs cs ~bm:(overrides cs.bm_cs 
      [(p0,p1),cs.a_cs p0 p1; (p1,p2),cs.a_cs p1 p2]) ()) in
  let csmax = 
    if (cs.bm_cs p1 p2 < cs.b_cs p1 p2 or cs.bm_cs p0 p1 < cs.b_cs p0 p1) 
    then None
    else Some (modify_cs cs ~am:(overrides cs.am_cs 
      [(p0,p1),cs.b_cs p0 p1; (p1,p2),cs.b_cs p1 p2]) ()) in
    cs1::cs2::(filter_some(csmin :: csmax:: (map cspq alldiag)));;


(* 
apply M-field deformation at (p,p+1)=(p1,p2) 
*)

(*
correction -- originally I had constraints on the straightness
of p0. By the beta lemma for nonreflexive local fans (or rather by the
monotonocity of dih in the opposite edge), the
azimuth angle is decreasing at p0 as p1 pivots in the direction of p2.
So cs0 is not used.
*)

let deform_4828966562 p0 p1 p2 cs =
  let ks = ks_cs cs in
  let diag = subtract ks [p0;p1;p2] in
  let _ = mem p1 ks or failwith "482:out of range" in
  let _ = (three <= cs.a_cs p0 p2) or     raise Unchanged in
  let _ = (cs.b_cs p1 p2 <= twoh0) or raise Unchanged in
  let _ = (cs.b_cs p0 p1 <= cstab) or raise Unchanged in
  let _ = (cs.a_cs p1 p2 < cs.bm_cs p1 p2) or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
(*  let _ = not(mem p0 cs.str_cs) or raise Unchanged in *)
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m q =  (cs.a_cs p1 q = cs.bm_cs p1 q) in
  let _ = forall (not o m) diag or raise Unchanged in
  let n q = (fourh0 < cs.b_cs p1 q) in
  let _ = forall n diag or raise Unchanged in
(*  let cs0 = modify_cs cs ~str:(sortuniq (p0::cs.str_cs)) () in  *)
  let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
  let cspq q = 
    if (cs.am_cs p1 q > cs.a_cs p1 q) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (p1,q,cs.a_cs p1 q)) ()) in 
    (* cs0:: *)  cs1::cs2::
      (filter_some( (cspq p2) ::(map cspq diag)));;

let deform_4828966562A p cs =
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  deform_4828966562 p0 p1 p2 cs;;

let deform_4828966562B p cs =
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  deform_4828966562 p2 p1 p0 cs;;

(*
Obsolete comment.
Use obtuse angle at p1 to avoid node straigtenings.
Uses "1117202051" and "4559601669" for obtuseness criteria.
//This one could cause infinite looping because we are deleting
//attributes str_cs at p0 p1.
Only use as a last resort.

Update: looping prevented.
It is a bad choice of words to call these deformations.
We are not actually
deforming, we are listing the constraints that prevent deformation.
The deformation might be a diagonal at its lower bound,
a straight node at p1, or an edge (p1,p2) at its lower bound.

If there are other conditions in the stable constraint system, that is ok.
If mem p2 cs.str_cs, then we don't need to eliminate it, because
we never deform it.  We are just saying that if an element of M_s
has a straight node at p2, then it must also have some other constraint
that prevents deformation.

So I have modified the code so not to remove p0 p2 from str_cs.

*)

(* check for obtuse angle at p1, based on delta_x4 calcs in main_estimate_ineq.hl,
    117202051, 2449601669, 45596016696. *)

let obtuse_crit cs p0 p1 p2 = 
    ((cs.bm_cs p0 p1= two) && (mem p2 cs.lo_cs)) or
      (mem p0 cs.lo_cs && mem p2 cs.lo_cs) or 
      ((cs.bm_cs p0 p1=two) && (mem p0 cs.lo_cs) );;

(*
Explanation of obtuse_sloc in the next deformation lemma.
Preconditions: hexagon (k=6) and three straights: p0,p4 and either p3 or p5,
  making an effective triangle.
Precondition: bm p1 p2 <= twoh0
The spherical law of cosines arg has been replaced with ineq 8430954724,
which does it as an ineq delta4_y < 0.

Here is the old argument:
new condition for obtuseness by spherical law of cosines.
if p3 p4 are both straight, then the three segments
give cos c <= cos(3.0 *. arc 2.52 2.52 2.) < -0.76 < -0.6.
If there are a total of 3 straights, then we have an effective triangle arclengths a,b,c=opposite.
assume p1 p2 are not straight and the side p1 p2 [2,2.52].
Then (0.206,0.685)=(cos(arc 2. 2. 2.52),cos(arc 2.52 2.52 2.)).
By the sloc, cos c - cos a cos b <= cos c + |cos a cos b|
  < cos c + 0.6 < 0. and we are obtuse.

If we try to make p5 straight, we have three consec straights p4,p5,p0.
making a chain of 4 edges. arc[2,2,2.52]*4 > pi, making a circular fan.
Hence the preconditions imply that p3 (not p5) is straight.
It is stated as it is, because we don't know if p0,p1,p2 are a inc sequence
or dec sequence, so only p4 is well-defined.
*)


let deform_4828966562_obtuse p0 p1 p2 cs =
  let ks = ks_cs cs in
  let diag = subtract ks [p0;p1;p2] in
  let _ = mem p1 ks or failwith "482:out of range" in
  let _ = (cstab <= cs.a_cs p0 p2) or     raise Unchanged in
  let _ = (cs.b_cs p1 p2 <= twoh0) or raise Unchanged in
  let _ = (cs.b_cs p0 p1 <= twoh0) or raise Unchanged in
  let obtuse = obtuse_crit cs p0 p1 p2 in
  let obtuse_sloc = 
    let p4 = funpow 3 (inc cs) p1 in
      (cs.k_cs = 6) &&
	(length (sortuniq cs.str_cs) = 3) && (subset [p0;p4] cs.str_cs) &&
	(cs.bm_cs p1 p2 <= twoh0) && not (mem p1 cs.str_cs) &&
	not (mem p2 cs.str_cs) in
  let _ = obtuse or obtuse_sloc or raise Unchanged in    
  let _ = (cs.a_cs p1 p2 < cs.bm_cs p1 p2) or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let m q =  (cs.a_cs p1 q = cs.bm_cs p1 q) in
  let _ = forall (not o m) diag or raise Unchanged in
  let n q = (fourh0 < cs.b_cs p1 q) in
  let _ = forall n diag or raise Unchanged in
  let cs_str1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cspq q = 
    if (cs.am_cs p1 q > cs.a_cs p1 q) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (p1,q,cs.a_cs p1 q)) ()) in 
    cs_str1:: filter_some (map cspq (p2::diag));;

let deform_4828966562A_obtuse p cs =
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  deform_4828966562_obtuse p0 p1 p2 cs;;

let deform_4828966562B_obtuse p cs =
  let p0 = dec cs p in
  let p1 = p in
  let p2 = inc cs p in
  deform_4828966562_obtuse p2 p1 p0 cs;;

(*

range calculations for 6843920790.
2.0 *. arc 2.52 2.52 2.0 > arc 2.0 2.0 2.38;;
2.0 *. arc 2.0 2.0 2.52 < arc1553;;
As the documentation for this inequality indicates, it is
applied to the triangle p4 p1 p2.

This also applies to some quads and triangles, 
implemented as separate deformations.

Edited: May 23, 2012.
need to consider more diagonals.  For instance (p1,p3) may bottom out.
*)


let deform_6843920790 p1 cs =
  let _ = (cs.k_cs = 5) or raise Unchanged in
  let ks = ks_cs cs in
  let p0 = dec cs p1  in
  let p2 = inc cs p1  in
  let p3 = inc cs p2  in
  let p4 = inc cs p3  in
  let _ = mem p1 ks or failwith "684:out of range" in
  let _ = (cs.am_cs p1 p2 = cstab && cstab = cs.bm_cs p1 p2)
    or raise Unchanged in
  let _ = (cstab <= cs.am_cs p1 p4) or raise Unchanged in
  let _ = (cstab <= cs.am_cs p4 p2) or raise Unchanged in
  let f (i,j) = (cs.bm_cs i j <= twoh0) in
  let _ = forall f [(p0,p1);(p2,p3);(p3,p4);(p4,p0)] or raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m (i,j) =  (cs.a_cs i j < cs.bm_cs i j) in
  let _ = forall m (alldiag cs) or raise Unchanged in
  let _ = m (p1,p2) or raise Unchanged in
  let n (i,j) = (fourh0 < cs.b_cs i j) in
  let _ = forall n (alldiag cs) or raise Unchanged in
  let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
    (* diags at p4 can be excluded, their lengths are fixed under def. *)
  let non_p4diag = [(p3,p1);(p3,p0);(p0,p2)] in 
  let cspq (i,j) = 
    if (cs.am_cs i j > cs.a_cs i j) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in 
     cs1::cs2::
      (filter_some(map cspq ((p1,p2)::non_p4diag)));;

(*
deformation 5512912661 functions in an almost identical manner
to 6843920790 on quads.  We take the union of the two deformations
here.
*)

let deform_6843920790_quad p1 p2 p3 p4 cs =
  let _ = (cs.k_cs = 4) or raise Unchanged in
  let ks = ks_cs cs in
  let adj (i,j) = (inc cs i = j) or (inc cs j = i) in
  let _ = forall adj [(p1,p2);(p2,p3);(p3,p4);(p4,p1)] or failwith "684q" in
  let _ = subset [p1;p2;p3;p4] ks or failwith "684:range" in
  let arc238 = arc two two 2.38 in
    (* ineq 684 conditions. *)
  let a2 = (two <= cs.am_cs p1 p2 && cs.bm_cs p1 p2 <= cstab) in
  let b2 = (arc238 <= arcmin cs(p1,p4) && cs.bm_cs p1 p4 <= cstab) in
  let c2 = (cstab <= cs.am_cs p2 p4 && arcmax cs(p2,p3)+.arcmax cs(p3,p4) <= arc1553) in
    (* ineq 5512912661 conditions. *)
  let a2' = (arc238 <= arcmin cs (p1,p2) && cs.bm_cs p1 p2 <= cstab) in
  let b2' = (two <= cs.bm_cs p1 p4 && cs.bm_cs p1 p4 <= twoh0) in
  let arc315 = arc two two (3.15/.1.26) in
  let c2' = (arc315 <= arcmin cs ( p2, p4) && 
	      ((arcmax cs(p2,p3)+.arcmax cs(p3,p4) <= arc1553) or
	       (arcmax cs (p2,p1) +.arcmax cs (p1,p4) <= arc1553))) in
  let _ = (a2 && b2 && c2) or (a2' && b2' && c2') or  raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m (i,j) =  (cs.a_cs i j < cs.bm_cs i j) in
  let _ = forall m (alldiag cs) or raise Unchanged in
  let _ = m (p1,p2) or raise Unchanged in
  let n (i,j) = (fourh0 < cs.b_cs i j) in
  let _ = forall n (alldiag cs) or raise Unchanged in
  let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
  let cspq (i,j) = 
    if (cs.am_cs i j > cs.a_cs i j) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in 
     cs1::cs2::
      (filter_some( [cspq (p1,p2);cspq (p1,p3)] ));;

let deform_6843920790_tri p1 cs =
  let p2 = inc cs p1 in
  let p3 = inc cs p2 in
  let _ = (cs.k_cs = 3) or raise Unchanged in
  let ks = ks_cs cs in
  let _ = mem p1 ks or failwith "684:tri" in
  let _ = subset [p1;p2;p3] ks or failwith "684:range" in
  let arc238 = arc two two 2.38 in
  let a2 = (two <= cs.am_cs p1 p2 && cs.bm_cs p1 p2 <= cstab) in
  let range(p,q) = (arc238 <= arcmin cs(p,q) && arcmax cs (p,q) <= arc1553) in
  let _ = (a2 && range(p3,p1) && range(p3,p2)) or  raise Unchanged in
  let _ = not(memj cs (p1,p2)) or raise Unchanged in
  let _ = not(mem p1 cs.str_cs) or raise Unchanged in
  let _ = not(mem p2 cs.str_cs) or raise Unchanged in
  let m (i,j) =  (cs.a_cs i j < cs.bm_cs i j) in
  let _ = m (p1,p2) or raise Unchanged in
  let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
  let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
  let cspq (i,j) = 
    if (cs.am_cs i j > cs.a_cs i j) then None
    else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in 
     cs1::cs2::
      (filter_some( [cspq (p1,p2)]));;

let deform_684_quadA p1 cs =
  let f j = funpow j (inc cs) in
  let p2 = f 1 p1 in
  let p3 = f 2 p1 in
  let p4 = f 3 p1 in
    deform_6843920790_quad p1 p2 p3 p4 cs;;

let deform_684_quadB p1 cs =
  let f j = funpow j (dec cs) in
  let p2 = f 1 p1 in
  let p3 = f 2 p1 in
  let p4 = f 3 p1 in
    deform_6843920790_quad p1 p2 p3 p4 cs;;



let deformations postfilter k = 
  let r f p cs = filter postfilter (f p cs) in
  let m d = map (r d) (0--(k-1)) in
  let rtl d cs = map restrict_type1_cs (d cs) in 
  let u =   [deform_ODXLSTC_cs;
	     deform_IMJXPHR_cs;
	     deform_NUXCOEA_cs;
	     deform_4828966562A_obtuse;
	     deform_4828966562B_obtuse;
	     deform_1834976373_A1_single;
	     deform_1834976373_A1_double;
	     deform_4828966562A;
	     deform_4828966562B;
	    deform_6843920790;
	    deform_684_quadA;
	    deform_684_quadB;
	    deform_6843920790_tri] in
    map rtl (List.flatten (map m u));;

let name_of_deformation k i =
  let names_hex = 
    ["odx";"imj";"nux";"482ao";"482bo";"1834s";
     "1834d";"482a";"482b";"684";"684qA";"684qB";"684t"] in
  let offset = i mod k in
  let s = i/k in
    (List.nth names_hex s) ^ "-" ^ (string_of_int offset);;


(*
****************************************
 SUBDIVISON
****************************************
*)

(* division acts on B-fields, shrinks domain.
   The global minima Ms remain global minima on smaller domain.
   We can keep adornments.
   Avoid Unchanged.  Pass through if c is out of range. 
   Restricts the domain if possible. *)

(*

Five domain possibilities, starting with [a,am,bm,b]

c <= a      --> unchanged
a < c <= am -->  cs2:[c,am,bm,b]
am < c < bm -->  cs1:[a,am,c,c]; cs2:[c,c,bm,b]
bm <= c < b -->  cs1:[a,am,bm,c]
b <= c      --> unchanged.

*)

let subdivide_cs p q c cs =
  let _ = (cs.k_cs > 3) or failwith "subdiv: triangle" in
  let _ = (0 <= p && p< cs.k_cs) or failwith "p out of range divide_cs" in
  let _ = (0 <= q && q< cs.k_cs) or failwith "q out of range divide_cs" in
  let _ = not(p = q) or failwith "subdiv: p q must be unequal" in
  let _ = not( memj cs (p,q)) or failwith "subdiv: (p,q) in J" in
  let e = alledge cs in
  let (p,q) = psort(p,q) in
  let _ = not(mem (p,q) e) or (cs.a_cs p q > two) or (cs.b_cs p q > twoh0) or 
    failwith "subdiv: cannot subdivide standard edge" in
  let _ = mem (p,q) (allpair cs) or failwith "subdivide range" in
  let _ = not (c = cs.am_cs p q) or not(mem(p,q) e) or failwith "subdiv: c is am on edge" in
  let a = cs.a_cs p q in
  let b = cs.b_cs p q in
  let am = cs.am_cs p q in
  let bm = cs.bm_cs p q in
  let cs1 = modify_cs cs ~b:(override cs.b_cs (p,q,c)) () in
  let cs2 = modify_cs cs ~a:(override cs.a_cs (p,q,c)) () in
  let _ =  (a < c  && c < b) or 
    (report ("warning:unchanged subdivide "^(string_of_float c));
		debug_cs:= cs::!debug_cs; true) in 
    if not (a < c && c < b) then [cs]
    else if bm <= c then [cs1]
    else if c <= am then [cs2]
    else (* am < c < bm *)
      let cs1' = modify_cs cs1 ~bm:(override cs.bm_cs (p,q,c)) () in
      let cs2' = modify_cs cs2 ~am:(override cs.am_cs (p,q,c)) () in
	[cs1';cs2'];;

let between_cs c cs (i,j) = 
  (cs.a_cs i j < c && c < cs.b_cs i j );;

let can_subdivide f_diag c cs =
  exists (between_cs c cs) (f_diag cs);;

let find_subdivide_edge f_diag c cs =
  let diag = f_diag cs in
  let ind = index true (map (between_cs c cs) (diag)) in
    List.nth diag ind;;

let subdivide_all_c_diag f_diag c = 
  let p = partition (can_subdivide f_diag c) in
  let rec sub init term =
    match init with
      | [] -> term 
      | cs::css -> 
	    let (i,j) = find_subdivide_edge f_diag c cs in
	    let kss = subdivide_cs i j c cs in
	    let (u,v) = p kss in
	      sub (u @ css) (v @ term) in
    fun init ->
      let (u,v) = p init in
	sub u v;;


(* claims serve as documentation for
   how the calculations are progressing from init_cs to terminal_cs.

   To do: make this into a formal deductive system.
   With rules of inference given by deformation rules, etc.

   a = what we assert to prove at that point, 
   b=what we leave for later 

   reduce remaining if r transfers to a.
     remove from b those that transfer to terminal.
     put residual b back in remaining.
  *)


let remaining = ref [];;
remaining := init_cs;; (* arrow *)

let claim_arrow(a,b) = 
  let r = !remaining in
  let transfers_to_a cs = exists (equi_transfer_cs cs) a in
  let r' = filter (not o transfers_to_a) r in
  let transfers_from_b cs = exists (C equi_transfer_cs cs) b in
  let r'' = filter (not o transfers_from_b) r' in
  let _ = remaining := (b @ r'') in
    !remaining;;

(*
****************************************
HEXAGONS
****************************************
*)

(* flow on hexagons.
  hex-std-
   subdivide all diags at stab.
   apply deformations,
     in cyclic repetition.
      setting aside all cs st ?bm diag <= stab.
   checking they are assert_stable_cs.
      
*)

(* claim  in this section goes as follows *)

(* todo: construct from generic subdivide *)

let hex_std_preslice_02,hex_std_preslice_03 = 
  let c1 = subdivide_cs 0 2 cstab hex_std_cs in
  let c2 = subdivide_cs 0 3 cstab hex_std_cs in
    hist (hd c1) "hex_std_preslice 02" , hist (hd c2) "hex_std_preslice 03";;

let hex_std_postslice = (* claim A  *)
  let cs = hex_std_preslice_02 in
  let d'= 0.616 in 
  let d = cs.d_cs -. d' in 
  let vv = slice_cs cs 2 0 d' d  false in
  let vv02 = map (C hist "slice_hex_02") vv in
  let cs = hex_std_preslice_03 in
  let d' = cs.d_cs /. 2.0 in
  let vv = slice_cs cs 0 3 d' d' false in
  let vv03 = transfer_union (map (C hist "slice_hex_03") vv) [] in
    vv02 @ vv03;;

let hex_assumptions = [hex_std_preslice_02;hex_std_preslice_03];;

claim_arrow([hex_std_cs],hex_assumptions);; (* claim B *)
claim_arrow(hex_assumptions,hex_std_postslice);; (* claim A *)



(* end claim *)

let ok_for_more_hex = 
  let is_hex cs = (cs.k_cs =6) in
  let terminals = filter_terminal is_hex in
    fun cs ->
      let _ = (cs.k_cs = 6) or failwith "ok:6" in
      let _ = assert_stable_cs cs  in
      let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
      let generic_at i = 
	let p0 = i in
	let p1 = inc cs i in
	let p2 = inc cs p1 in
	let p3 = inc cs p2 in
	  not (subset[p0;p1;p2;p3] cs.lo_cs && subset[p1;p2] cs.str_cs) in
      let bg = forall generic_at (ks_cs cs) in
      let bunfinished = not (exists (transfer_to cs) terminals) in
	bstr && bg && bunfinished;;

let hex_deformations = deformations ok_for_more_hex 6;;

let name_of_deformation_hex = name_of_deformation 6;;

let filtered_subdivide postfilter init  = 
  let sub = subdivide_all_c_diag (alldiag) cstab init in
    filter postfilter sub;;

let rec apply_first dl cs = 
  match dl with
      [] -> failwith ( "all deformations fail")
    | d::dls ->
	try
	  d cs
	with Unchanged -> apply_first dls cs;;

let preslice_ready cs = 
  if (cs.k_cs=4) && (cs.d_cs=0.467) 
  then
    (cs.bm_cs 0 2 = three) && (cs.bm_cs 1 3 = three)
  else  
    exists (fun (i,j) -> cs.bm_cs i j <= cstab ) (alldiag cs);;

let rec general_loop2 df postfilter c active stab_diags =
    if (c <= 0) or length stab_diags > 0 then (active,stab_diags) 
    else match active with
	[] -> ([],stab_diags) 
      | cs::css -> 
	    try 
	      let kss = apply_first df cs in
	      let (u,v) = partition preslice_ready kss in
	      let u' = filter postfilter u in
		general_loop2 df postfilter (c-1) (v @ css) (u' @  stab_diags)
	    with 
               Failure s -> (active,stab_diags);;

let execute_hexagons() =   (* claim B *)
  let postfilter = not o (equi_transfer_to_list hex_assumptions) in 
  let hex_ultra = filtered_subdivide postfilter [hex_std_cs] in
  let hex_loop = general_loop2 hex_deformations postfilter in
    (([],[]) =hex_loop 200000 hex_ultra []);;

(* if true, it means that all hexagons have been reduced
   to the one terminal hexagon, together with two cases of hex_std_preslice
*)

(* time execute_hexagons ();; 117secs. 
   working in svn:2819,2824m,2829m,2830m,2832m *)


(*
****************************************
PENTAGONS
****************************************
*)

(* flow on pentagons is the same as for hexagons *)


let ok_for_more_pent =
  let terminals = filter_terminal (fun cs -> cs.k_cs = 5) in
    fun cs ->
      let _ = (cs.k_cs = 5) or failwith "ok:5" in
      let _ = assert_stable_cs cs in
      let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
      let sph_tri_ineq i = 
	if (length cs.str_cs < 2) then true
	else 
	  let p0 = i in
	  let p1 = inc cs p0 in
	  let p2 = inc cs p1 in
	  let p3 = inc cs p2 in
	  let p4 = inc cs p3 in
	    if not(subset [p1;p2] cs.str_cs) then true
	    else
	      let e03min = arcmin cs (p0,p1) +. 
		arcmin cs (p1,p2) +. arcmin cs (p2,p3) in
	      let e03max = arcmax cs (p3,p4) +. arcmax cs (p4, p0) in
		e03min <= e03max in
      let bunfinished = not (exists (transfer_to cs) terminals) in
	bstr && bunfinished && forall sph_tri_ineq (ks_cs cs);;


let pent_deformations = deformations ok_for_more_pent 5;;

let name_of_deformation_pent = name_of_deformation 5;;

(* these are the cases handled in the pent section *)

let pent_init = 
  filter (fun cs ->(5=cs.k_cs)) (init_cs @ hex_std_postslice);;

(* it doesn't matter how our assumptions are generated,
   as long as they are eventually discharged. *)

let pent_composite_cs = mk_unadorned (
  5,
  0.616,
  a_pro two two cstab 5,
  a_pro cstab twoh0 upperbd 5,"pent composite") ;;

let pent_assumptions = 
  let pent_comp_rediag_cs (p,q) = 
    let cs = pent_composite_cs in
      modify_cs cs
	~b:(override cs.b_cs (p,q,cstab))
	~bm:(override cs.bm_cs (p,q,cstab)) () in
  let alld = alldiag pent_std_cs in
  let ffh ((p,q), cs) = subdivide_cs p q cstab cs in
  let preslices = List.flatten (map ffh (cart alld pent_init)) in
  let pent_comb = map pent_comp_rediag_cs alld in
  let cstab_preslices = filter preslice_ready (pent_comb @ preslices) in
  let union_cstab_preslices = transfer_union cstab_preslices [] in 
  let pa =  map (C hist "preslice") union_cstab_preslices in
    pa;;

(* we introduce pent_composite_ultra_cs 
   which combines all the pent_init cases,
   then we run execute_pentagons to get rid of it too.
*)

let pent_composite_ultra_cs = (* claim C *)
  let pl = filtered_subdivide 
    (not o (equi_transfer_to_list pent_assumptions)) 
    (pent_composite_cs::pent_init) in
    transfer_union pl [];;

(* claim C *)
claim_arrow(pent_init,pent_composite_ultra_cs @pent_assumptions);;

(* pent section main claim *)

let execute_pentagons() =  (* claim D *)
  let pent_filter = not o (equi_transfer_to_list pent_assumptions) in 
  let pent_loop = general_loop2 pent_deformations pent_filter in
    (([],[])=pent_loop 200000 (pent_composite_ultra_cs) []);;


(* 
   time execute_pentagons();;

   if true, it means that all pentagons have been reduced
   to the one terminal pentagon, together with cases of pent_assumptions
   worked 36sec. in svn:2821, May 20, 2012.
   svn 2834m.
*)

claim_arrow(pent_composite_ultra_cs,pent_assumptions);; (* claim D *)


(*
****************************************
ECHELON QUADS
****************************************
*)

(* These are the quadrilaterals that need subdivision in
the 'ultra' range of diagonals >= cstab.
There are two ways of doing this, either by picking the smaller
diagonal (echelon B) or picking the "better" diagonal, even if
not the shortest.

This should be used as the last resort on a quadrilateral.
It fails unless the slice reduces to terminal cases.
*)

let delta_am_diag2_neg cs d = 
  let y01 = cs.am_cs 0 1 in
  let y03 = cs.am_cs 0 3 in
  let y21 = cs.am_cs 2 1 in
  let y23 = cs.am_cs 2 3 in
    Sphere_math.delta_y d y01 y03 d y23 y21 < 0.0;;

let delta_am_diag_neg cs p1 d = (* d from p1 to p3, am from p2 to p4 *)
  let p2 = inc cs p1 in
  let p3 = inc cs p2 in
  let p4 = inc cs p3 in
  let y12 = cs.am_cs p1 p2 in
  let y14 = cs.am_cs p1 p4 in
  let y32 = cs.am_cs p3 p2 in
  let y34 = cs.am_cs p3 p4 in
  let am = cs.am_cs p2 p4 in
    Sphere_math.delta_y d y12 y14 am y34 y32 < 0.0;;

let check_echelon_precondition cs = 
  try
    let _ = (cs.k_cs = 4) or failwith "echelon1" in
    let fixed (i,j) = (cs.am_cs i j = cs.bm_cs i j) in
    let _ = forall fixed (alledge cs) or failwith "echelon:edge" in
    let _ = forall (fun(i,j) ->cs.am_cs i j >= cstab) (alldiag cs) or 
      failwith "ech:diag" in
    let _ = (cs.js_cs = []) or failwith "ech:j" in 
      true
  with Failure _ -> false;;

let upper_echelonA =
  let assumptions = filter_terminal all in
  let transfer = x_equi_transfer_to_list assumptions in
    fun (p1,(dval,diag)) cs ->
      let _ = check_echelon_precondition cs or raise Unchanged in
      let _ = delta_am_diag_neg cs p1 diag or raise Unchanged in
      let p2 = inc cs p1 in
      let p3 = inc cs p2 in
      let dval' = cs.d_cs -. dval in
      let css = filter 
	(fun cs -> cs.bm_cs p1 p3 <= diag) (subdivide_cs p1 p3 diag cs) in
      let css' = List.flatten 
	(map (fun cs -> slice_cs cs p1 p3 dval dval' false) css) in
      let css'' = filter (not o transfer) css' in
      let _ = (css''=[]) or raise Unchanged in
	[];;

let upper_echelonC = 
  let assumptions = filter_terminal all in
  let transfer = x_equi_transfer_to_list assumptions in
    fun p1 cs ->      
      let diag341 = 3.41 in
      let _ = check_echelon_precondition cs or raise Unchanged in
      let p2 = inc cs p1 in
      let p3 = inc cs p2 in
      let p4 = inc cs p3 in
      let edge_is c (i,j) = (cs.am_cs i j = c) in
      let rhs = [(p1,p2);(p2,p3)] in
      let lhs = [(p3,p4);(p4,p1)] in
      let has_one c e = (length(filter(edge_is c) e) = 1) in
      let _ = (has_one two rhs) or raise Unchanged in
      let _ = (has_one two lhs) or raise Unchanged in
      let _ = (has_one cstab rhs) or raise Unchanged in
      let _ = (has_one twoh0 lhs) or (has_one cstab lhs) or raise Unchanged in
      let css = subdivide_cs p1 p3 diag341 cs in
      let (css1,css2) = partition (fun cs -> cs.bm_cs p1 p3 <= diag341) css in
      let _ = (length css1 = 1) or failwith "ech:C" in
      let cs1 = restrict_type2_cs (hd css1) in
      let dval = 0.4759 in
      let dval' = cs.d_cs -. dval in
      let css' = slice_cs cs1 p1 p3 dval dval' false in
      let css'' = filter (not o transfer) css' in
      let _ = (css'' =[]) or raise Unchanged in
	css2;;

let upper_echelonB = 
  let assumptions = filter_terminal all in
  let transfer = x_equi_transfer_to_list assumptions in
    fun diag cs ->
      let _ = check_echelon_precondition cs or raise Unchanged in
      let _ = delta_am_diag2_neg cs diag or raise Unchanged in
      let css = filter 
	(fun cs -> cs.bm_cs 0 2 <= diag or cs.bm_cs 1 3 <= diag)
	(subdivide_all_c_diag alldiag diag [cs]) in
      let css02,css13 = partition 
	(fun cs -> cs.bm_cs 0 2 <= diag) css in
      let dcases = if (diag=3.41) then [0.0759;0.4759] else [0.2759] in
      let op(a,b) = (b,a) in
      let dv = map (fun d -> (d,cs.d_cs -. d)) dcases in
      let dv' = (dv @ map op dv) in
      let f (p,q) (d,d') cs = 
	let css' = slice_cs cs p q d d' false in
	let css'' = filter (not o transfer) css' in
	let _ = (css''=[]) or raise Unchanged in
	  () in
	try
	  ignore(map (apply_first (map (f(0,2)) dv')) css02);
	  ignore(map (apply_first (map (f(1,3)) dv')) css13);
	  []
	with Failure _ -> raise Unchanged;;

let upper_echelon  = 
  (* we get the case data from echelon data in terminal_cs. *)
  let dataA = cart (0--3)
    [(0.0759,3.41);(0.4759,3.41);(0.2759,3.339);(0.2759,3.62)] in
  let dataB = [3.41;3.339;3.62] in
  let dataC = (0--3) in
  let cases = map upper_echelonA dataA @ map upper_echelonB dataB @
    map upper_echelonC dataC in
    fun cs -> 
      apply_first cases cs;;



(*
****************************************
QUADRILATERALS
****************************************
*)

(* this handles quad_pro_cs modulo the return cs's *)

let (quad_477_preslice_short,quad_477_preslice_long) = (* claim E *)
  let preslices = subdivide_all_c_diag alldiag cstab [quad_pro_cs] in
  let vv =  (map (C hist "preslice pro") preslices) in
  let p = filter (fun cs -> cs.b_cs 0 2 = cstab) 
    (subdivide_cs 0 2 cstab quad_pro_cs) in
  let ww = transfer_union (p @ vv) [] in
    partition preslice_ready ww;;

(* This is the case where ears are required -- finally! *)

let execute_quad_to_ear() = (* claim E' *)
  let transfer = x_equi_transfer_to_list (filter_terminal all) in
  let _ = (length quad_477_preslice_short = 1) or failwith "handle 477" in
  let cs = hd quad_477_preslice_short in
  let mk_ear = true in
  let inc2 = funpow 2 (inc cs) in
  let f i =  slice_cs cs i (inc2 i) (0.11) (0.477 -. 0.11) mk_ear in
  let ind = filter (can f) (0--3) in
  let g i =  forall (transfer) (f i) in
    exists g ind;;


(* let quad_remaining = !remaining;;
remaining :=quad_remaining;;
 *)

claim_arrow([quad_pro_cs],quad_477_preslice_long);;  (* claim E *)

let triquad_assumption = 
[
  mk_unadorned (
    4,
    0.513,
    overrides (cs_adj two cstab 4) [((1,2),twoh0);((2,3),twoh0)],
    overrides (cs_adj twoh0 upperbd 4) [((1,2),cstab);((2,3),cstab)],
    "triquad assumption"
  );
  mk_unadorned (
    4,
    0.513,
    overrides (cs_adj two cstab 4) [((1,2),twoh0);((0,3),twoh0)],
    overrides (cs_adj twoh0 upperbd 4) [((1,2),cstab);((0,3),cstab)],
    "triquad assumption"    
  );
  mk_unadorned (
    3,
    0.4278,
    cs_adj twoh0 cstab 3,
    cs_adj cstab upperbd 3,
    "tri assumption"
  );
(*
  mk_unadorned (
    3,
    0.103,
    funlist [(0,1),twoh0] two,
    funlist [(0,1),sqrt8] twoh0,
    "td 2 1");
  mk_unadorned (
    3,
    0.0,
    cs_adj two cstab 3,
    cs_adj twoh0 cstab 3,
    "td 3 0");
*)
];;

claim_arrow([],triquad_assumption);;  (* claim - no justification required *)

(* 477 has been handled already. 
  We make the case involving the ear part of the terminal
   set so that we don't need to deal with ears any further.
  can make it an assumption *)

let terminalj_cs = quad_477_preslice_short @ (filter_terminal all);;

let terminal_quad =  
  triquad_assumption @ terminalj_cs;;

(* see calc in main_estimate.hl. If opposite edges are too short
   the angle cannot be straight.  *)
let can_be_straight_2485876245b cs p1 = 
  if (not (cs.k_cs = 4)) then true
  else 
    let p0 = dec cs p1 in
    let p2 = inc cs p1 in
    let p3 = inc cs p2 in
      not (cstab <= cs.am_cs p1 p3 &&
	     cs.bm_cs p1 p2 <= cstab &&
	     cs.bm_cs p1 p0 <= cstab &&
	     cs.bm_cs p2 p3 <= twoh0 &&
	     cs.bm_cs p3 p0 <= twoh0);;

let ok_for_more_tri_quad assumptions = 
  let terminal_transfer = x_equi_transfer_to_list assumptions in
  fun cs -> 
  let _ = (mem cs.k_cs [3;4]) or failwith "ok:3,4" in
  let b467_2485876245a = 
    if (cs.d_cs=0.467) && (cs.k_cs = 4) &&
      forall (fun (i,j)-> cs.bm_cs i j <= twoh0) [(0,1);(1,2);(2,3);(3,0)] &&
      forall (fun (i,j)-> cs.am_cs i j >= three) [(0,2);(1,3)]
    then ( cs.str_cs = []) else true in
  let deltay cs = Sphere_math.delta_y (* if neg then geometric impossibility *)
	(cs.bm_cs 0 1) (cs.bm_cs 0 3) (cs.am_cs 0 2)
	(cs.bm_cs 2 3) (cs.bm_cs 1 2) (cs.am_cs 1 3) >= 0.0 in
  let b4 = (cs.k_cs = 4) in
  let b4a = if b4 then deltay cs else true in
  let is_477 = b4 && (cs.d_cs=0.477) &&
      forall (fun (i,j)-> cs.bm_cs i j <= twoh0) [(0,1);(1,2);(2,3);(3,0)] &&
      forall (fun (i,j)-> cs.am_cs i j >= cstab) [(0,2);(1,3)] in
  let b477 =
    if is_477
    then
      let b477a =  cs.str_cs = [] in
       let b477c = 
	not(exists (equi_transfer_cs cs) quad_477_preslice_short) in
	(b477a && b477c) 
    else true in      
  let sph_tri_ineq p1 = 
    let p0 = dec cs p1 in
    let p2 = inc cs p1 in
    let p3 = inc cs p2 in
      if (not (mem p1 cs.str_cs)) then true (* in tri str=[] *)
      else
	arcmin cs (p0,p1) +. arcmin cs (p1,p2) <
	  arcmax cs (p0,p3) +. arcmax cs (p3,p2) in
  let sph_tri_ineq2 p1 = (* rather ad hoc to kill a case *)
    let p0 = dec cs p1 in
    let p2 = inc cs p1 in
    let p3 = inc cs p2 in
      if (not (mem p1 (intersect cs.str_cs cs.lo_cs))) then true
      else
	cs.am_cs p0 p1 < cs.bm_cs p0 p3 or
	  cs.am_cs p1 p2 < cs.bm_cs p2 p3 in 
  let _ = assert_stable_cs cs in
  let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
  let bunfinished = not (terminal_transfer cs) in
    bstr && b4a && bunfinished && b467_2485876245a && 
      b477 && forall sph_tri_ineq (cs.str_cs) &&
      forall sph_tri_ineq2 (intersect cs.str_cs cs.lo_cs) &&
     forall (can_be_straight_2485876245b cs) (cs.str_cs);;

let special_quad_init = 
  quad_diag_cs :: quad_477_preslice_long;;

let execute_special_quads = (* claim F *)
  let quad_deformations = 
    deformations (ok_for_more_tri_quad terminal_quad) 4 in
  let terminal_transfer = x_equi_transfer_to_list terminal_quad in
    fun () ->
  let quad_loop = general_loop2 quad_deformations
    (not o terminal_transfer) in
    (([],[])=quad_loop 200000  special_quad_init []);;

(*
execute_special_quads();;
*)

claim_arrow(special_quad_init,[]);;  (* claim F *)

(* execute_special_quads is true, it means that 
   all special quads have been reduced
   to terminal quads, 
   worked in svn:2821, May 20, 2012.
   svn:2824.
*)


(* now the general case remains
   policies and procedure:
   - no more use of js_cs and ears.
   - determine d by slice_dpq procedure.

   The order is extremely important!
   1- if the cs transfers to a terminal or is not "ready" for more, 
   then were are done with it.
   2- if stab is between lower and upper bounds of a diag, then subdivide it
   3- if sqrt8 is between lower and upper bounds of a diag, then sub
   4- if there is a diag [2h0,sqrt8], then slice it.
   5- if there is a diag [sqrt8,cstab], then slice it.
   6- if a deformation applies, then apply it.
   7- do "upper echelon" treatment of quads (subdivisions on edges > 3.01)

*)

let ok_for_more assumptions = 
  let ok3 = ok_for_more_tri_quad assumptions in
    fun cs ->
  ((cs.k_cs = 6) && (ok_for_more_hex cs)) or 
  ((cs.k_cs=5) && (ok_for_more_pent cs)) or
  ((mem cs.k_cs [3;4]) && (ok3 cs));;


let handle_general_case skip8 assumptions = 
  let ok = ok_for_more assumptions in
  let dff k = deformations (fun cs -> true) k in
  let defs = [[];[];[];dff 3;dff 4;dff 5;dff 6] in
  let inrange cs a b (p,q) = (a <= cs.am_cs p q && cs.bm_cs p q <= b) in
  let allunderstable cs = 
    if(cs.k_cs < 4) then []
    else filter (fun (i,j)-> (cs.b_cs i j <= cstab)) (allpair cs) in
  fun cs -> 
  let alld = alldiag cs in
  (* 1- *)
  if not(ok cs) then []
 (*  else if exists (equi_transfer_cs cs) terminal_quad then [] *)
    (* 2,3- *)
  else if not(skip8) && (can_subdivide allunderstable sqrt8 cs) 
  then (subdivide_all_c_diag allunderstable sqrt8 [cs])
  else if (can_subdivide allunderstable twoh0 cs) 
  then (subdivide_all_c_diag allunderstable twoh0 [cs])
  else 
    try (* 4- *)
      let (p,q) = Lib.find  (inrange cs twoh0 sqrt8) alld in
	slice_std cs p q 
    with Failure _ ->
      try (* 5 *)
	let (p,q) = Lib.find  (inrange cs sqrt8 cstab) alld in
	  slice_std cs p q 
      with Failure _ ->
	try (* 6 *)
	  if (can_subdivide alldiag cstab cs) 
	  then (subdivide_all_c_diag alldiag cstab [cs])
	  else
	    let k = cs.k_cs in
	    let dl = List.nth defs k in
	      apply_first dl cs 
	with Failure _ -> 
	  try (* 7 *)
	    upper_echelon cs 
	  with Failure _ -> failcs(cs, "no handler found");;


let handle_loop skip8 assumptions = 
  let handle_one = handle_general_case skip8 assumptions in
  let rec handle_loop_rec c ls =
    if (c<=0) then ls 
    else match ls with
      | [] -> []
      | cs::css -> 
	  let v = handle_one cs in
	  let (a,b) = partition (fun cs -> cs.k_cs > 4) (v @ css) in
	  let (b',b'') = partition (fun cs -> cs.k_cs = 3) b in
	    handle_loop_rec (c-1) ( a @ b' @ b'') in
    handle_loop_rec;;


let r_init = !remaining;;  

let execute_triquad () =  (* claim G*)
  let b1 = ([] = (handle_loop false terminal_quad 50000)  r_init) in
  let b2 = ([] = (handle_loop true terminalj_cs 50000)  triquad_assumption) in
    b1 && b2;;

claim_arrow(r_init,terminal_quad);; (* claim G *)
claim_arrow(triquad_assumption,terminalj_cs);; (* claim G *)
claim_arrow(terminalj_cs,[]);; (* claim E *)

let all_cases_done = (!remaining = []);;

let execute() = 
  (time execute_hexagons() &&
    time execute_pentagons() &&
    time execute_quad_to_ear() &&
    time execute_special_quads() &&
    time execute_triquad() &&
    all_cases_done);;


end;;


(* scratch area *)

(*
let handle = handle_general_case true terminalj_cs;;

debug_cs := [];;
execute_triquad();;
let cs1 = (hd !debug_cs);;
report_cs cs1;;
let cs2 = upper_echelon cs1;;
report_cs (hd cs2);;
let cs3 = hd(handle (hd cs2));;
report_cs cs3;;
1;;

let hl = time (handle_loop true terminalj_cs 50000)  tri_cases_left;;
let cs_term = hd (!debug_cs);;
let cs_init = List.nth quad_cases_left 1;;
report_cs cs_term;;
let handle1 = handle_general_case true terminalj_cs;;
let handle1' =  catch_failure cs_term handle1;;
let kl = debug_trace handle1' cs_init cs_term;;
List.length kl;;
let cs' = List.nth kl 4;;
report_cs cs';;
handle1 cs';;
slice_std cs' 1 3 ;;
calc_d_std cs' 1 3;;

let kl = handle1 cs1;;
map report_cs kl;;
let cases_left = 
  filter (not o (x_equi_transfer_to_list (filter_terminal all)) terminal_quad;;
*)