theory Example_Lorenz_Classic
imports
  "../Numerics/Euler_Affine"
  "../Numerics/Optimize_Float"
begin


abbreviation "lorenz_real \<equiv> \<lambda>(x::real, y::real, z::real). (10 * (y + - z), x * (28 + - z) + - y, x*y + - 8 * z * inverse 3)"
approximate_affine lorenz lorenz_real

abbreviation "lorenz_d_real \<equiv> \<lambda>(x1, x2, x3). \<lambda>(a, aa, b).
  (10 * aa + - 10 * b, a * (28 + - x3) + - x1 * b + - aa, x1 * aa + a * x2 + - 8 * b * inverse 3)::real*real*real"

lemma lorenz_fderiv:
  "(lorenz_real has_derivative lorenz_d_real x) (at x)"
  by (auto intro!: derivative_eq_intros ext simp: split_beta inverse_eq_divide)

approximate_affine lorenz_d lorenz_d_real

abbreviation "lorenz_dd_real \<equiv>
  \<lambda>(x1::real, x2::real, x3::real).
    \<lambda>(dx, dy, dz).
    \<lambda>(dx', dy', dz'). (0*x1, - (dx * dz') + - dx' * dz, dx' * dy + dx * dy')
    ::real*real*real"
approximate_affine lorenz_dd lorenz_dd_real

lemma lorenz_fderiv2: "((\<lambda>x. lorenz_d_real x dx) has_derivative
    lorenz_dd_real x dx) (at x)"
  by (auto intro!: derivative_eq_intros ext simp: split_beta')

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 lorenz!: 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 simp: split_beta')
  done

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

(* aform_of_ivl (FloatR 5 (- 2), FloatR 146 (- 6)) (FloatR 49 (- 5), FloatR 149 (- 6)) *)

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,
    min_section_distance = 1/2/2/2,
    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>"

(* width(splitting) Euler, atol, reduction *)
definition lorenz_optns'::"int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> int \<Rightarrow> _"
  where "lorenz_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, -1) (- FloatR 100 (- 6)) 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 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)))"

(* uses approximate inital sets, but this is irrelevant for the timings *)
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, 36)) (aform_of_ivl 0 (decimal u eu, decimal u eu, 0))]"

export_code lorenzc nat_of_integer int_of_integer in SML module_name Lorenz file "lorenz_core.ML"

end
