validated_integ3 drops the remainder each step inherits from the previous steps #208

            # Include remainder of initial condition
            x1N[i].pol.coeffs[1].rem = z#rem0[i]

using TaylorModels, TaylorModels.ValidatedInteg, IntervalArithmetic, TaylorSeries, Printf
setprecision(BigFloat, 256)

# Geodesic flow on the unit sphere: x' = v, v' = -(|v|^2/|x|^2) x.  From x0 = e_z with v0 = (v1, v2, 0), c = |v0|,
# the exact solution is x(t) = cos(c t) e_z + sin(c t) v0 / c, for every point of the initial box below.
function spray!(dx, x, p, t)
    q = (x[4]^2 + x[5]^2 + x[6]^2) / (x[1]^2 + x[2]^2 + x[3]^2)
    dx[1] = x[4]; dx[2] = x[5]; dx[3] = x[6]
    dx[4] = -q * x[1]; dx[5] = -q * x[2]; dx[6] = -q * x[3]
    nothing
end
const A, B = 0.7, 0.7 + 2pi / 256
const X0 = [interval(0.0), interval(0.0), interval(1.0),
            interval(cos(B), cos(A)), interval(sin(A), sin(B)), interval(0.0)]

# Disable the remainder absorption, so the Taylor variables keep addressing their own initial points
# (a function enclosure); the pointwise check below is only meaningful then.
@eval TaylorModels.ValidatedInteg _abs_rems!(vTMN::Vector{TaylorModelN{N,T,S}}) where {N,T,S} = nothing

# Optionally restore the dropped remainder: the package's own `_validation3!` with one line replaced.
function use_fix!(fix::Bool)
    src = split(read(joinpath(pkgdir(TaylorModels), "src", "valid_integ", "validated_integ3.jl"), String), '\n')
    i0 = findfirst(l -> startswith(l, "function _validation3!"), src)
    i1 = i0 + findfirst(==("end"), src[i0+1:end])
    fun = join(src[i0:i1], '\n')
    fix && (fun = replace(fun, "x1N[i].pol.coeffs[1].rem = z#rem0[i]" => "x1N[i].pol.coeffs[1].rem = rem0[i]"))
    Core.eval(TaylorModels.ValidatedInteg, Meta.parse(fun))
end

function check(label)
    sol = validated_integ3(spray!, X0, 0.0, 0.5, 2, 12, 1e-15, JetSpace(4, ["ξ$i" for i in 1:6]);
                           maxsteps = 1000, parse_eqs = false)
    tTM, qTM = expansion_point(sol), get_xTM(sol)
    fails, worst = 0, BigFloat(0)
    for n in 2:length(sol), a in (-0.9, 0.0, 0.9), b in (-0.9, 0.0, 0.9), fr in (0.25, 0.75)
        dom = domain(qTM[1, n]); dt = inf(dom) + fr * (sup(dom) - inf(dom))
        xi = interval.([0.0, 0.0, 0.0, a, b, 0.0])
        v1 = BigFloat(mid(X0[4])) + a * BigFloat(radius(X0[4]))
        v2 = BigFloat(mid(X0[5])) + b * BigFloat(radius(X0[5]))
        c, t = sqrt(v1^2 + v2^2), BigFloat(tTM[n]) + BigFloat(dt)
        exact = [sin(c * t) * v1 / c, sin(c * t) * v2 / c, cos(c * t)]
        for i in 1:3
            e = evaluate(evaluate(qTM[i, n], interval(dt)), xi)
            excess = max(BigFloat(inf(e)) - exact[i], exact[i] - BigFloat(sup(e)), BigFloat(0))
            fails += excess > 0; worst = max(worst, excess)
        end
    end
    @printf("%-14s steps=%d failures=%d/162 worst excess=%.3e  remainder(x1) first step %.2e, last step %.2e\n",
            label, length(sol) - 1, fails, Float64(worst),
            radius(remainder(qTM[1, 2])), radius(remainder(qTM[1, end])))
end

use_fix!(false); check("as released")
use_fix!(true);  check("with the fix")