#=
testODE:
- Julia version: 
- Author: ymocquar
- Date: 2019-11-13
=#
using HOODESolver
using DifferentialEquations
using LinearAlgebra
using Plots
using Random
function twoscales_solve( par_u0::PrepareU0, order, t, nb)
end
function fctMain(n_tau)
    seed=123321
    Random.seed!(seed)
    u0 = rand(BigFloat, 4)
    B = 2rand(BigFloat, 4, 4) - ones(BigFloat, 4, 4)
    println("seed = $seed B=$B")
    nbeps=48
    nbmaxtest=3
    paseps = big"10.0"^0.125
    ordmin=3
    ordmax=14
    ordpas = 1
    t_max = big"1.0"
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    for order = ordmin:ordpas:ordmax
        y = ones(Float64, nbeps, nbmaxtest)
        x=zeros(Float64,nbeps)
        for i_eps=1:nbeps
            epsilon = big"1.0"/paseps^i_eps
            x[i_eps] = epsilon
            fct = u -> B*u
            parphi = PreparePhi(epsilon, n_tau, A, fct, B)
            println("prepareU0 eps=$epsilon n_tau=$n_tau")
            @time par_u0 = PrepareU0(parphi, order+2, u0)
            solref = getexactsol(parphi, u0, t_max)
            eps_v = convert(Float32,epsilon)
            println("epsilon = $eps_v solref=$solref")
            for indc=1:nbmaxtest
                nb = 10*10^indc
                @time parGen = PrepareTwoScalePureAB(nb, t_max, order, par_u0)
                @time sol = twoscales_pure_ab(parGen, only_end=true)
                println("solref=$solref")
                println("nb=$nb sol=$sol")
                diff=solref-sol
                println("nb=$nb dt=$(1.0/nb) normInf=$(norm(diff,Inf)) norm2=$(norm(diff))")
                y[i_eps, indc] = norm(diff,Inf)
                println("epsilon=$epsilon result=$y")
            end
        end
        labels=Array{String,2}(undef, 1, nbmaxtest)  
        for i=1:nbmaxtest
            nb = 10*10^i
            labels[1,i] = " delta t,order=$(1/nb), $order "
            mindiff = minimum(y[:,i])
            bornediff = min(mindiff*1e20,1)
            for j=1:nbeps
                if y[j,i] >= bornediff
                    y[j,i] = NaN
                end
            end
        end
        gr()
        p=Plots.plot(
    x,
    y,
    xlabel="epsilon",
    xaxis=:log,
    ylabel="error",
    yaxis=:log,
    legend=:topleft,
    label=labels,
    marker=2
)
        prec_v = precision(BigFloat)
        Plots.savefig(p,"out/r_$(prec_v)_$(order)_$(n_tau)_epsilon_lgn_2.pdf")
    end
end
fctMain(32)
