section"Example: Lorenz attractor"
theory Lorenz_Approximation
imports
  "../../Numerics/Euler_Affine"
(*  "../../Numerics/Euler_Affine_Code" *)
  "../../Numerics/Optimize_Float"
  Initials
  "~~/src/HOL/Library/Parallel"
begin
text {* \label{sec:lorenz} *}
(*
double F0(double x, double y, double z)
{
  return lu*x - k1*(x + y)*z;
}


/* Returns the second component of the vector field F */
double F1(double x, double y, double z)
{
  return lss*y + k1*(x + y)*z;
}


/* Returns the third component of the vector field F */
double F2(double x, double y, double z)  
{
  return ls*z + (x + y)*(k2*x + k3*y);
}
*)
(* plot/tuckers_coeffs.py says:
lu = 11.8277234512 = FloatR 6658416365912628 -49
lss = -22.8277234512 = FloatR -6425432926773530 -48
ls = -2.66666666667 = FloatR -6004799503160661 -51
k1 = 0.288554928412 = FloatR 5198143472295670 -54
k2 = 2.18277234512 = FloatR 4915166360050299 -51
k3 = -1.28277234512 = FloatR -5777093055467152 -52
*)

abbreviation "lu == FloatR 6658416365912628 (- 49)"
abbreviation "lss == FloatR (- 6425432926773530) (- 48)"
abbreviation "ls == FloatR (- 6004799503160661) (- 51)"
abbreviation "k1 == FloatR 5198143472295670 (- 54)"
abbreviation "k2 == FloatR 4915166360050299 (- 51)"
abbreviation "k3 == FloatR (- 5777093055467152) (- 52)"

subsection {* Lorenz equations in Affine arithmetic *}

text {* TODO: lorenz -> lorenz-affine, lorenz-real -> lorenz*}
abbreviation "lorenz_real \<equiv> \<lambda>(x::real, y::real, z::real).
  (lu*x + -k1*(x + y)*z, lss*y + k1*(x + y)*z, ls*z + (x + y)*(k2*x + k3*y))"

approximate_affine lorenz lorenz_real

schematic_lemma lorenz_fderiv_sch:
  fixes lu lss ls k1 k2 k3 shows
  "(lorenz_real has_derivative
    ((\<lambda>(y1, y2, y3) (dy1, dy2, dy3). (?f1 y1 y2 y3 dy1 dy2 dy3, ?f2 y1 y2 y3 dy1 dy2 dy3, ?f3 y1 y2 y3 dy1 dy2 dy3)) x))
    (at x within X)"
  by (auto intro!: derivative_eq_intros simp: split_beta)

abbreviation "lorenz_d_real \<equiv> (\<lambda>(y1, y2, y3).
      \<lambda>(dy1, dy2, dy3).
         (lu * dy1 + - (k1 * (y1 + y2) * dy3 + k1 * (dy1 + dy2) * y3),
          lss * dy2 + (k1 * (y1 + y2) * dy3 + k1 * (dy1 + dy2) * y3),
          ls * dy3 + ((y1 + y2) * (k2 * dy1 + k3 * dy2) + (dy1 + dy2) * (k2 * y1 + k3 * y2))))"

lemma lorenz_fderiv: "(lorenz_real has_derivative lorenz_d_real x) (at x within X)"
  using lorenz_fderiv_sch
  by (simp add: sign_simps)

approximate_affine lorenz_d lorenz_d_real

schematic_lemma lorenz_fderiv2_sch:
  fixes lu lss ls k1 k2 k3 shows
  "((\<lambda>x. lorenz_d_real x dx) has_derivative
    ((\<lambda>(y1, y2, y3) (dy1, dy2, dy3) (dy1', dy2', dy3').
      (?f1 y1 y2 y3 dy1 dy2 dy3 dy1' dy2' dy3',
      ?f2 y1 y2 y3 dy1 dy2 dy3 dy1' dy2' dy3',
      ?f3 y1 y2 y3 dy1 dy2 dy3 dy1' dy2' dy3')) x dx))
    (at x within X)"
  by (auto intro!: derivative_eq_intros simp: split_beta)

abbreviation "lorenz_dd_real \<equiv> \<lambda>(y1::real, y2::real, y3::real) (dy1, dy2, dy3) (dy1', dy2', dy3').
  (- k1 * ((dy1' + dy2') * dy3 + (dy1 + dy2) * dy3') + 0 * y1,
   k1 * ((dy1' + dy2') * dy3 + (dy1 + dy2) * dy3'),
   (dy1' + dy2') * (k2 * dy1 + k3 * dy2) + (dy1 + dy2) * (k2 * dy1' + k3 * dy2'))"

abbreviation "lorenz_dd_real' \<equiv> \<lambda>(y1::real, y2::real, y3::real) (dy1, dy2, dy3) (dy1', dy2', dy3').
  (- k1 * (dy1' + dy2') * dy3 - k1 * (dy1 + dy2) * dy3',
         k1 * ((dy1' + dy2') * dy3 + (dy1 + dy2) * dy3'),
         (dy1' + dy2') * (k2 * dy1 + k3 * dy2) +
          (dy1 + dy2) * (k2 * dy1' + k3 * dy2'))"

lemma lorenz_dd_eq: "lorenz_dd_real' = lorenz_dd_real"
  by (force simp: algebra_simps)

lemma lorenz_fderiv2: "((\<lambda>x. lorenz_d_real x dx) has_derivative lorenz_dd_real x dx) (at x within X)"
  using lorenz_fderiv2_sch[of dx x X]
  unfolding lorenz_dd_eq[symmetric]
  by (simp add: algebra_simps split_beta')

approximate_affine lorenz_dd lorenz_dd_real

abbreviation "lorenz_ivp \<equiv> \<lambda>optns args. uncurry_options lorenz optns (hd args) (tl args)"
abbreviation "lorenz_d_ivp \<equiv> \<lambda>optns args. uncurry_options lorenz_d optns (hd args) (hd (tl args)) (tl (tl args))"
abbreviation "lorenz_dd_ivp \<equiv> \<lambda>optns args. uncurry_options lorenz_dd optns (hd args) (hd (tl args)) (hd (tl (tl args))) (tl (tl (tl args)))"

interpretation aform_approximate_ivp
  lorenz_ivp
  lorenz_d_ivp
  lorenz_dd_ivp
  lorenz_real
  lorenz_d_real
  lorenz_dd_real
  apply unfold_locales
  apply (rule Joints2_JointsI)
  unfolding list.sel
  apply (rule lorenz, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (rule lorenz_fderiv)
  apply (rule lorenz_fderiv2)
  apply (rule Joints2_JointsI)
  apply (rule lorenz_d, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (rule Joints2_JointsI)
  apply (rule lorenz_dd, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (auto intro!: continuous_intros)
  done

definition "constrain x y = (if x = y then () else ())"
lemma [code]: "constrain x y = ()" by (simp add: constrain_def)

definition "lorenz_optns =
    \<lparr>
    precision = 53,
    tolerance = FloatR 1 (- 4),
    stepsize  = FloatR 1 (- 8),
    iterations = 40,
    halve_stepsizes = 10,
    widening_mod = 40,
    rk2_param = FloatR 1 (- 1),
    max_tdev_thres = FloatR 1 (- 3),
    presplit_summary_tolerance = FloatR 1 (- 2),
    collect_granularity = FloatR 1 (- 3),
    start_section = Sctn (0, 0, -1) (- (215/2/2/2)),
    checkpoint_sections = [Sctn (0, 0, 1) 27],
    stop_sections = [Sctn (0, 0, -1) (- 27), Sctn (0, 0, -1) (FloatR (- 25) (- 8))],
    safe_appr = \<lambda>_. True,
    snap_sections = 1/2/2/2/2,
    min_section_distance = 0,
    next_section_distance_factor = 25,
    next_section_turn_distance_factor = 0,
    printing_fun = (let re = print_rectangle 0 0 0; af = print_aform; _ = constrain re af;
      both = (\<lambda>c. let _ = re c in af c) in both),
    irect_mod = 8,
    result_fun = ivls_result 20 40,
    max_intersection_step = FloatR 1 (- 10),
    adaptive_atol = FloatR 1 (-30),
    method_id = 2
  \<rparr>"

definition lorenz_optns'::"nat \<Rightarrow> nat \<Rightarrow> _"
  where "lorenz_optns' p w = lorenz_optns \<lparr>
      irect_mod := p,
      collect_granularity := FloatR 2 (- int w),
      max_tdev_thres := FloatR 2 (- int w),
      adaptive_atol := FloatR 1 (- int (w * 2))\<rparr>"

definition mirror_irects
  where "mirror_irects =
    map (\<lambda>irect. case irect of [i, j, k] \<Rightarrow> [if j < 0 then - i else i , abs j, k] | irect \<Rightarrow> irect)"

definition "lorenz_irects p w i j =
  do {
    let optns = lorenz_optns' p w;
    stops \<leftarrow> poincare_irects_dres (lorenz_ivp) (lorenz_d_ivp) (lorenz_dd_ivp) optns [i, j, 27 * 2 ^ p];
    let (returns, aborted) = List.partition (\<lambda>x. snd (snd (Inf_aform x)) = 27) stops;
    let irects = irects_of_apprs optns returns;
    let irects = mirror_irects irects;
    let _ = (if length aborted > 0 then
      (let
        _ = println (STR ''= Aborted close to the origin: '');
        ivl_abort = ivl_of_apprs aborted;
        _ = print_eucl (Inf_aform ivl_abort);
        _ = print (STR '' -- '');
        _ = print_eucl (Sup_aform ivl_abort);
        _ = println (STR [])
      in ())
      else ());
    let _ = println (STR ''= Initial Rectangles of Return:'');
    let _ = print_irects irects;
    dRETURN irects
  }"

(* value [code] "lorenz_irects 8 5 198 541" *)
(* value [code] "lorenz_irects 8 8 373 588" *)


definition lorenzc_optns::"int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> _"
  where "lorenzc_optns w m a r = lorenz_optns \<lparr>
      collect_granularity := FloatR 1 w,
      max_tdev_thres := FloatR 1 w,
      adaptive_atol := FloatR 1 a,
      start_section := if r\<noteq>0 then Sctn (0, 0, -1) (- 24) else
        Sctn (0, -1) (- FloatR 1 100),
      next_section_distance_factor := if r ~= 0 then 100 else FloatR 1 100,
      next_section_turn_distance_factor := if r ~= 0 then 0 else FloatR 1 100,
      max_intersection_step := FloatR 1 (- 2),
      method_id := nat m\<rparr>"

definition "decimal x ex = (if ex \<ge> 0 then x * 10 ^ nat ex else real_divl 53 x (10^nat (-ex)))"

definition lorenzc::"[int, int, int, int,int, int, int, int,int, int] \<Rightarrow> _" where
  "lorenzc w m a r x ex y ey u eu =
  poincare_distance_d_dres (lorenz_ivp) (lorenz_d_ivp) (lorenz_dd_ivp)
  (lorenzc_optns w m a r)
  [msum_aform 53 (aform_of_point (decimal x ex, decimal y ey, 27)) (aform_of_ivl 0 (decimal u eu, decimal u eu, 0))]"

subsection {* Extraction of Code *}
text {* Code is extracted and can be used when imported using lorenz.ML *}
export_code lorenzc lorenz_irects nat_of_integer int_of_integer in SML module_name Lorenz file "lorenz_core.sml"

end
