Set up PR

GaussTuranExampleTemp/example.jl

@@ -0,0 +1,63 @@
+using Optim
+using PreallocationTools
+using Integrals
+using DoubleFloats
+
+include("high_order_derivatives.jl")
+
+# This fix is needed to avoid crashing: https://github.com/JuliaLang/julia/pull/54201
+include("linear_algebra_fix.jl")
+
+FT = Double64
+n = 5
+s = 1
+a = FT(0.0)
+b = FT(1.0)
+
+"""
+    Example functions ϕⱼ from [1]
+"""
+function f(x::T, j::Int)::T where {T<:Number}
+    pow = if j % 2 == 0
+        j / 2 - 1
+    else
+        (j - 1) / 2 - 1 / 3
+    end
+    x^pow
+end
+
+deriv! = @generate_derivs(2)
+I = Integrals.GaussTuran(
+    f,
+    deriv!,
+    a,
+    b,
+    n,
+    s;
+    optimization_options = Optim.Options(;
+        x_abstol = FT(1e-250),
+        g_tol = FT(1e-250),
+        show_trace = true,
+        show_every = 100,
+        iterations = 10_000,
+    ),
+)
+
+# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| first N functions fⱼ
+max_int_error_upper =
+    maximum(abs(I(x -> f(x, j)) - I.cache.rhs_upper[j]) for j = 1:I.cache.N)
+@show max_int_error_upper
+# max_int_error_upper = 3.0814879110195774e-32
+
+# Integration error |I(fⱼ) - ₐ∫ᵇfⱼ(x)dx| last n functions fⱼ
+max_int_error_lower = maximum(
+    abs(I(x -> f(x, j)) - I.cache.rhs_lower[j-I.cache.N]) for
+    j = (I.cache.N+1):(I.cache.N+I.cache.n)
+)
+@show max_int_error_lower
+# max_int_error_lower = 2.8858134286698342e-30
+
+# Example with eˣ
+exp_error = abs(I(Base.exp) - (Base.exp(1) - 1))
+@show exp_error;
+# exp_error = -6.436621416953051e-18

GaussTuranExampleTemp/high_order_derivatives.jl

@@ -0,0 +1,57 @@
+using ForwardDiff
+
+"""
+    Generate a function `derivs!(out, f, x)` which computes the 0th up to the max_orderth derivative
+    of the scalar-to-scalar function f at x and stores them in ascending derivative order in `out`.
+    Hence `out` must be at least of length `max_order + 1`. 
+"""
+macro generate_derivs(max_order::Int)
+    # Create nested dual number of required depth (arg_0, …, arg_{max_order})
+    arg_assignments = [:(arg_0 = x)]
+    for i = 1:max_order
+        arg_name = Symbol("arg_$i")
+        prev_arg_name = Symbol("arg_$(i-1)")
+        push!(
+            arg_assignments,
+            :($arg_name = ForwardDiff.Dual{Val{$i}}($prev_arg_name, one($prev_arg_name))),
+        )
+    end
+
+    # Unpack the results
+    arg_max = Symbol("arg_$max_order")
+    res_unpacks = [:(res_0 = f($arg_max))]
+    for i = 1:max_order
+        res_name = Symbol("res_$i")
+        res_prev_name = Symbol("res_$(i-1)")
+        push!(res_unpacks, :($res_name = only($res_prev_name.partials)))
+    end
+
+    # Assign the results
+    out_assignments = Expr[]
+    for i = 0:max_order
+        res = Symbol("res_$i")
+        res_temp = Symbol("$(res)_temp_0")
+        push!(out_assignments, :($res_temp = $res))
+        # Create temporary variables to get to
+        # res_{i}.value.value.value…
+        for j = 1:(max_order-i)
+            res_temp = Symbol("$(res)_temp_$j")
+            res_temp_prev = Symbol("$(res)_temp_$(j-1)")
+            push!(out_assignments, :($res_temp = $res_temp_prev.value))
+        end
+        res_temp = Symbol("$(res)_temp_$(max_order - i)")
+        push!(out_assignments, :(out[$(i + 1)] = $res_temp))
+    end
+
+    # Construct the complete function definition
+    func_def = quote
+        function derivs!(out, f, x::T)::Nothing where {T<:Number}
+            $(arg_assignments...)
+            $(res_unpacks...)
+            $(out_assignments...)
+            return nothing
+        end
+    end
+
+    return func_def
+end