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

subsection {* Oil reservoir in Affine arithmetic *}
text {*\label{sec:exampleoil}*}
abbreviation "oil_real \<equiv> \<lambda>(y::real, z::real). (z, z*z + -3 * inverse (inverse 1000 + y*y))"
approximate_affine oil oil_real

lemma oil_deriv_ok: fixes y::real
shows "1 / 1000 + y*y = 0 \<longleftrightarrow> False"
proof -
  have "1 / 1000 + y*y > 0"
    by (auto intro!: add_pos_nonneg)
  thus ?thesis by auto
qed

abbreviation "oil_d_real \<equiv> \<lambda>(y::real, z) (dy, dz). (dz, 2 * dz * z  + 6 * (inverse (inverse 1000 + y*y) * (dy * (y * inverse (inverse 1000 + y*y)))))"

lemma oil_fderiv: "((\<lambda>(y::real, z::real). (z, z * z + -3 * inverse (inverse 1000 + y * y))) has_derivative oil_d_real x)
   (at x within X)"
  by (auto intro!: derivative_eq_intros simp: oil_deriv_ok split_beta inverse_eq_divide)

approximate_affine oil_d oil_d_real

abbreviation "oil_dd_real \<equiv>
  \<lambda>(x::real, y::real). \<lambda>(dx::real, dy::real). \<lambda>(dx'::real, dy'::real).
    (0*x, 2 * dy * dy' + (6 * (inverse (inverse 1000 + x * x) *
    (dx * (dx' * inverse (inverse 1000 + x * x) +
    (- x) * (inverse (inverse 1000 + x * x) * (2 * (dx' * x)) * inverse (inverse 1000 + x * x))))) +
    (- 12) * (inverse (inverse 1000 + x * x) * (dx' * (x * (inverse (inverse 1000 + x * x) *
    (dx * (x * inverse (inverse 1000 + x * x)))))))))"
approximate_affine oil_dd oil_dd_real

lemma oil_fderiv2: "((\<lambda>x. oil_d_real x dx) has_derivative oil_dd_real x dx) (at x)"
  by (auto intro!: derivative_eq_intros simp: split_beta' oil_deriv_ok)

abbreviation "oil_ivp \<equiv> \<lambda>optns args. uncurry_options oil optns (hd args) (tl args)"
abbreviation "oil_d_ivp \<equiv> \<lambda>optns args. uncurry_options oil_d optns (hd args) (hd (tl args)) (tl (tl args))"
abbreviation "oil_dd_ivp \<equiv> \<lambda>optns args. uncurry_options oil_dd optns (hd args) (hd (tl args)) (hd (tl (tl args))) (tl (tl (tl args)))"

interpretation oil!: aform_approximate_ivp
  oil_ivp
  oil_d_ivp
  oil_dd_ivp
  oil_real
  oil_d_real
  oil_dd_real
  apply unfold_locales
  apply (rule Joints2_JointsI)
  unfolding list.sel
  apply (rule oil, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (rule oil_fderiv)
  apply (rule oil_fderiv2)
  apply (rule Joints2_JointsI)
  apply (rule oil_d, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (rule Joints2_JointsI)
  apply (rule oil_dd, assumption, assumption)
  apply (drule length_set_of_apprs, simp)--"TODO: prove in affine-approximation"
  apply (auto intro!: continuous_intros simp: split_beta' oil_deriv_ok)
  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 "oil_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 (-1, 0) 2,
    checkpoint_sections = [],
    stop_sections = [Sctn (-1, 0) 10],
    safe_appr = \<lambda>_. True,
    snap_sections = 1,
    min_section_distance = 1/2/2/2,
    next_section_distance_factor = FloatR 1 30,      (* dont do intersections *)
    next_section_turn_distance_factor = FloatR 1 30, (* dont do intersections *)
    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 oil_optns'::"nat \<Rightarrow> nat \<Rightarrow> _"
  where "oil_optns' p w = oil_optns \<lparr>
      irect_mod := p,
      collect_granularity := FloatR 2 (- int w),
      max_tdev_thres := FloatR 1 10, (* do not split *)
      adaptive_atol := FloatR 1 (- int w)\<rparr>"

definition "oil_limit w = poincare_distance_d_dres (oil_ivp) (oil_d_ivp) (oil_dd_ivp) (oil_optns' 0 w) [aform_of_point (10, 0)]"

export_code oil_limit nat_of_integer int_of_integer in SML module_name Lorenz file "oil_core.ML"

end
