text:
#= R. Stoner

 tests zonation scripts; one unzoned case without rad. damage
 run this file from COAST directory
 when debugging make sure you start julia in COAST directory as well

=#

using Revise
import Pkg
Pkg.activate(pwd())
using COAST
using COAST.Zoned
using CSV
using DataFrames
using Test
using JuMP
using Ipopt
using Interpolations
using Plots
# check if COAST was loaded
@testset "COAST.jl" begin
    @test loaded_COAST()==true
end

@testset "zonation.jl" begin
    #=

     Unzoned U-Pb rutile from Cherniak, 00
     no radiation damage

    =#
    # TODO: remove superfluous rad. damage params
    # A-Th/He rad. damage params; these aren't used here
    alpha = 0.04672
    c0 = 0.39528
    c1 = 0.01073
    c2 = -65.12969
    c3 =  -7.91715
    rmr0 = 0.79
    eta_q = 0.91
    L_dist = 8.1*1e-4 # cm (!)

    # general params
    U238 = 1.0*1e-6
    R = 1.9872*1e-3
    L = 1.0*1e-4 # m
    n_T_pts = 1000
    n_iter = 300.0
    Ea = 59.752
    logD0_a2 = log10(3.9e-10/(L^2))

    # preprocess concentrations
    (U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = ppm_to_atoms_per_volume(U238,0.0,density=3.20)

    # set up time-temperature (t-T) path
    T1 = collect(LinRange(900.0,600.0,ceil(Int,n_T_pts/2)).+273.15)
    T2 = collect(LinRange(600.0,400.0,floor(Int,n_T_pts/2)).+273.15)
    T = vcat(T1,T2)
    times1 = collect(LinRange(70.0,50.0,ceil(Int,n_T_pts/2)+1).*3.1558e7*1e6)
    times2 = collect(LinRange(49.9,0.01,floor(Int,n_T_pts/2)).*3.1558e7*1e6)
    times = vcat(times1,times2)

    mass_pb_t1 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,0.83,eta_q,L_dist,
                                         0.0,0.0,0.0,R,Ea,logD0_a2,n_iter,
                                         U238_mol,U238_V,U235_mol,U235_V,
                                         Th232_mol,Th232_V,L,times...,T...)
    pre_pb_t1 = (8*(U238_mol*exp(50.7*1e6*sec_in_yrs/τ38)-U238_mol)+ # HeFTy output
                 7*(U235_mol*exp(50.7*1e6*sec_in_yrs/τ35)-U235_mol))
    # 50.8 - analytical output
    print("Unzoned Cherniak rutile numerical test numerical is $(mass_pb_t1)! \n")
    print("Unzoned Cherniak rutile test predicted HeFTy is $(pre_pb_t1)! \n")
    @test isapprox(mass_pb_t1/pre_pb_t1,1.0; atol = 1e-2) # <1 Ma error
    # end unzoned test

    #=

     can test zonation function can handle array of U values?
     U-Pb rutile (Cherniak,00)

    =#
    # time stepping setup
    Nt = 39
    t_end = 1.0
    dt = t_end/(Nt-1)

    # set up time-temperature (t-T) path
    t1 = collect(LinRange(70.0,50.0,ceil(Int,Nt/2)).*3.1558e7*1e6)
    t2 = collect(LinRange(49.9,0.01,floor(Int,Nt/2)).*3.1558e7*1e6)
    t = vcat(t1,t2)
    T1 = collect(LinRange(900.0,600.0,ceil(Int,Nt/2)).+273.15)
    T2 = collect(LinRange(600.0,400.0,floor(Int,Nt/2)).+273.15)
    T = vcat(T1,T2)
    Ea = 250.0*1e3
    D0 = 3.9*1e-10
    L = 1e-4
    Nx = 30
    r = LinRange(0.0,L,Nx)
    U238 = 1.0
    U238_new = U238*ones(Nx)
    Nt = length(t)

    # test function uses splatting; JuMP doesn't like arrays
    Pb06forw = zonation_forward(Ea,COAST.Rjoules,D0,0.0,0.0,L,Nt,Nx,t...,T...,r...,U238_new...)
    Pbtot = [0.0]
    r = LinRange(0.0,L,Nx)

    # integrate/convert profile to total Pb concentration
    # check date w. HeFTy date
    for i in 2:Nx
        Pbtot[1] += pi*(r[i]^3-r[i-1]^3)*(Pb06forw[i]+Pb06forw[i-1])/2
    end
    Pbtot = Pbtot[1]/(pi*L^3)
    print("FD: total lead predicted is $(Pbtot)\n")
    Pb_heft = (1.0*(U238*exp(50.7*1e6*COAST.sec_in_yrs/COAST.τ38)-U238))
    print("HeFTy: total lead predicted is $(Pb_heft)\n")
    @test isapprox(Pbtot/Pb_heft,1.0; atol = 5e-2) # higher Nx and Nt lead to better match

    #=

     synthetic test of constrained optimizer (IPOPT) used by COAST
     U-Pb rutile (Cherniak,00)
     array of U values but unzoned; cooling path
     same as previous test

    =#

    # set up t-T path
    Nt = 39
    t_end = 1.0
    dt = t_end/(Nt-1)
    T1 = collect(LinRange(900.0,600.0,ceil(Int,Nt/2)).+273.15)
    T2 = collect(LinRange(600.0,400.0,floor(Int,Nt/2)).+273.15)
    T = vcat(T1,T2)
    t1 = collect(LinRange(70.0,50.0,ceil(Int,Nt/2)).*3.1558e7*1e6)
    t2 = collect(LinRange(49.9,0.01,floor(Int,Nt/2)).*3.1558e7*1e6)
    t = vcat(t1,t2)

    # Cherniak, 00 diffusion params
    # setup U and grain size
    Ea = 250.0*1e3
    D0 = 3.9*1e-10
    L = 1e-4
    Nx = 30
    r = LinRange(0.0,L,Nx)
    U238 = 1.0
    U238 = U238*ones(Nx)

    # create functions to be passed to IPOPT
    # IPOPT-friendly funcs
    jac_constraint = create_jac_constraint(Ea,COAST.Rjoules,D0,0.0,0.0,Nt,Nx,t,r,U238)
    constraint_zon = create_constraint_zon(Ea,COAST.Rjoules,D0,U238,0.0,0.0,Nt,Nx,t,r)

    # bounds temperature
    T_L = zeros(Float64,Nt)
    T_U = 2900*ones(Float64,Nt) # T_U has to be defined, set bound arbitrarily high

    # synthetic Pb values
    initPb = [0.008857378583018752, 0.008856812929792447, 0.00885511177504213, 0.008852259612534526, 0.008848230516985477, 0.008842987011603362, 0.008836478776030107, 0.008828640808589248, 0.008819390873808581, 0.008808625978040217, 0.00879621748087797, 0.008782004239101911, 0.008765782846776095, 0.008747293506359331, 0.008726199226256422, 0.008702054713998092, 0.008674259257549808, 0.008641984671522813, 0.008604064447778814, 0.008558822633270739, 0.008503808762407928, 0.008435383783207719, 0.008348058711971884, 0.008233386603361441, 0.008077941780303029, 0.007859134069521883, 0.007534875513179412, 0.007010548962254766, 0.005967321901702219, 0.0]

    # bounds Pb (i.e. constraint bounds)
    Pb_L = initPb - 0.02*initPb
    Pb_U = initPb + 0.02*initPb

    prob = createProblem(
        Nt, # should match length T_L, T_U otherwise julia may crash
        T_L, # lower bound on temperature
        T_U, # upper bound
        Nx,
        Pb_L, # lower bound of Pb concentration; i.e. constraints
        Pb_U,
        Nt*Nx, # number elements in constraint jacobian
        Nt*Nx,
        smooth_zon_objective, # objective function; sum smoothness (multiobjective not yet possible)
        constraint_zon,
        grad_smooth_zon_objective,
        jac_constraint,
    )
    addOption(prob, "hessian_approximation", "limited-memory") # hessian rarely (if ever) called; approximation is ok
    addOption(prob, "tol", 300.0) # smoothness not as much of an issue here
    addOption(prob,"print_level",3)
    prob.x = 1.1*T
    status=solveProblem(prob) # this can crash julia if prob set up incorrectly (see above)
    solve_successful = :Solve_Succeeded
    @test Ipopt.ApplicationReturnStatus[status]==solve_successful

    # test if previous test results match these
    Pbs = zeros(Nx) # preallocate
    constraint_zon(T,Pbs)

    @test minimum(isapprox.(Pbs[1:end-1]./Pb06forw[1:end-1],1.0;atol=1e-10))
    print("Test diffusion func matches main diffusion func! \n")

    #=

    Test with Smye, 2018 rutile dataset

    =#
    # Cherniak, 00 diffusion params
    # setup U and grain size

    # define constant params
    Ea = 250.0*1e3 # J/mol => 59.75 kcal/mol
    D0 = 3.9*1e-10 # m^2/s => 3.9e-6 cm^2/s
    U238 = 1.0
    Nt = 50

    # load data and preprocess
    ngrains = 1
    subdir = "/test"
    filename = "/smyedata.csv"
    raw_data = CSV.read(pwd()*subdir*filename,DataFrame,header=false)
    data = Matrix(raw_data)
    ngrains = size(data,2)÷3
    # separate data by grain
    data = [data[:,3*igrain-2:(3*igrain)] for igrain in 1:ngrains]
    data = Dict((1:ngrains).=>data)
    #data[1] = data[1][1:11,:] # smye cut out data deeper in grain (plateau?)
    # convert from um->m
    for val in values(data)
        val[:,1]/=1e6
    end

    # fill in data to center of grain
    rad = 73.0*1e-6 # maximum grain size
    Nx = zeros(Int64,ngrains) # preallocate for number of nodes in grain
    for (key, value) in data
        dr = diff(value[1:2,1])[1]
        r_added_min = maximum(value[:,1])
        r_added = collect(r_added_min:dr:(rad+dr))
        n_added = length(r_added)
        # fill in all other values as well
        Pb06U38_grain_center=fill(value[end,2],n_added)
        sigPb06U38_grain_center=fill(value[end,3],n_added)
        data_grain_center = hcat(r_added,Pb06U38_grain_center,sigPb06U38_grain_center)
        data[key] = vcat(value,data_grain_center)
        data[key] = vcat([0.0 0.0 0.0],data[key])
        Nx[key] = size(data[key],1)
    end

    # correct model to model end date
    end_model_Ma = 800.0
    Pb_corr_end_model = 1.0*(exp(COAST.lambda_38*COAST.sec_in_yrs*end_model_Ma*1e6))-1.0

    # pull out individual components from data for clarity
    U238_new = Dict{Int64,Vector{Float64}}()
    r = Dict{Int64,Vector{Float64}}()
    Pbmeas = Dict{Int64,Vector{Float64}}()
    sigmeas = Dict{Int64,Vector{Float64}}()
    for (key,value) in data
        value[2:end,2] = value[2:end,2].-Pb_corr_end_model
        push!(r,key=>value[:,1])
        push!(U238_new,key=>U238*ones(length(value[:,1])))
        push!(Pbmeas,key=>reverse(value[:,2]))
        push!(sigmeas,key=>reverse(value[:,3]))
    end

    t1 = collect(LinRange(1100.0,1000.0,ceil(Int,Nt/2)).*3.1558e7*1e6)
    dt = -maximum(diff(t1))/(3.1558e7*1e6)
    t2 = collect(LinRange(1000.0-dt,end_model_Ma,floor(Int,Nt/2)).*3.1558e7*1e6)
    t = vcat(t1,t2)
    T1 = collect(LinRange(830.0,870.0,ceil(Int,Nt/2)).+273.15)
    T2 = collect(LinRange(870.0,400.0,floor(Int,Nt/2)).+273.15)
    T = vcat(T1,T2)

    upper_bound_arr = 1100.0*ones(Nt)
    lower_bound_arr = 673.0*ones(Nt)
    lower_bound_arr[1] = 1100.0
#     lower_bound_arr[1]=1100.0
#     upper_bound_constr = reverse(data[:,2]+3*data[:,3])
#     lower_bound_constr = reverse(data[:,2]-3*data[:,3])
#     reassign_ind = 10
#     upper_bound_constr[end-reassign_ind+1:end] = reverse(data[1:reassign_ind,2]+2.0*data[1:reassign_ind,3])
#     lower_bound_constr[end-reassign_ind+1:end] = reverse(data[1:reassign_ind,2]-2.0*data[1:reassign_ind,3])
#     lower_bound_constr[lower_bound_constr.<0.0] .= 0.0

#     datasig=vcat(2*data[1:reassign_ind,3],3*data[reassign_ind+1:end,3])

#     # create functions to be passed to IPOPT
#     # IPOPT-friendly funcs
#     jac_constraint = Zoned.create_jac_constraint(Ea,COAST.Rjoules,D0,0.0,0.0,Nt,Nx,t,r,U238_new)
#     constraint_zon = Zoned.create_constraint_zon(Ea,COAST.Rjoules,D0,U238_new,0.0,0.0,Nt,Nx,t,r)

#     # better for underconstrained setups!
#     prob_smye = createProblem(
#         Nt, # should match length T_L, T_U otherwise julia may crash
#         lower_bound_arr, # lower bound on temperature
#         upper_bound_arr, # upper bound
#         Nx,
#         lower_bound_constr, # lower bound of Pb concentration; i.e. constraints
#         upper_bound_constr,
#         Nt*Nx, # number elements in constraint jacobian
#         Nt*Nx,
#         Zoned.smooth_zon_objective, # objective function; sum smoothness (multiobjective not yet possible)
#         constraint_zon,
#         Zoned.grad_smooth_zon_objective,
#         jac_constraint,
#     )

#     addOption(prob_smye, "hessian_approximation", "limited-memory") # hessian rarely (if ever) called; approximation is ok
#     addOption(prob_smye, "tol", 1.0) # smoothness not as much of an issue here
#     addOption(prob_smye,"print_level",5)
#     prob_smye.x = copy(T)

#     status=solveProblem(prob_smye)

#     # Value 2
#     # create functions to be passed to IPOPT
#     # IPOPT-friendly funcs
#     U238_new = reshape(U238_new,:,1)
#     Pbmeas = reshape(reverse(data[:,2]),:,1)
#     sigmeas = reshape(reverse(data[:,3]),:,1)
#     r = reshape(r,:,1)
    empty_dict = Dict(keys(U238_new).=>0.0.*values(U238_new))
    zon_loglike_objective = Zoned.create_zon_loglike_objective(Ea,COAST.Rjoules,D0,U238_new,empty_dict,empty_dict,Nt,Nx,t,r,Pbmeas,sigmeas,ngrains)
    grad_loglike_objective = Zoned.create_grad_loglike_objective(Ea,COAST.Rjoules,D0,U238_new,empty_dict,empty_dict,Nt,Nx,t,r,Pbmeas,sigmeas,ngrains)

    n_constr = Nt-1
    upper_bound_constr = 100*ones(n_constr)
    lower_bound_constr = 0.0*ones(n_constr)

    prob2 = createProblem(
        Nt, # should match length T_L, T_U otherwise julia may crash
        lower_bound_arr, # lower bound on temperature
        upper_bound_arr, # upper bound
        n_constr,
        lower_bound_constr, # lower bound of Pb concentration; i.e. constraints
        upper_bound_constr,
        Nt*n_constr, # number elements in constraint jacobian
        Nt*n_constr, # number elements in constraint hessian; not needed in approximation
        zon_loglike_objective, # objective function; sum smoothness (multiobjective not yet possible)
        Zoned.decrease_constraint,
        grad_loglike_objective,
        Zoned.jac_decrease_constraint,
    )
    addOption(prob2, "hessian_approximation", "limited-memory") # hessian rarely (if ever) called; approximation is ok
    addOption(prob2, "tol", 1e-4) # smoothness not as much of an issue here
    addOption(prob2,"print_level",5)
    prob2.x = copy(T)
    status=solveProblem(prob2)

    constraint_zon = create_constraint_zon(Ea,COAST.Rjoules,D0,U238_new[1],0.0,0.0,Nt,Nx[1],t,r[1])
    Pbs = zeros(Nx[1]) # preallocate
    constraint_zon(prob2.x,Pbs)

    constraint_zon = create_constraint_zon(Ea,COAST.Rjoules,D0,U238_new[2],0.0,0.0,Nt,Nx[2],t,r[2])
    Pbs2 = zeros(Nx[2]) # preallocate
    constraint_zon(prob2.x,Pbs2)

    # plot results
    # plot(t/3.156e13,prob2.x.-273.15,label="fitted t-T",legend=:topleft,lw=2)
    # plot!(t/3.156e13,T.-273.15,label="starting t-T",lw=2)
    # xlabel!("time (Ma)")
    # ylabel!("temperature (°C)")
    # title!("decreasing: GB119C−42")

    # scatter(r[1][1:end]*1e6,data[1][1:end,2],ribbon=reverse(sigmeas[1][1:end]),xflip=true)
    # plot!(r[1][1:end]*1e6,reverse(Pbs[1:end]),lw=2,label="fitted values")

end


text:
#####
##### Products between fields and grid lengths
#####

@inline Δx_uᶠᶜᵃ(i, j, k, grid, u) = @inbounds Δxᶠᶜᵃ(i, j, k, grid) * u[i, j, k]
@inline Δx_vᶜᶠᵃ(i, j, k, grid, v) = @inbounds Δxᶜᶠᵃ(i, j, k, grid) * v[i, j, k]
@inline Δx_cᶜᶜᵃ(i, j, k, grid, c) = @inbounds Δxᶜᶜᵃ(i, j, k, grid) * c[i, j, k]

@inline Δy_uᶠᶜᵃ(i, j, k, grid, u) = @inbounds Δyᶠᶜᵃ(i, j, k, grid) * u[i, j, k]
@inline Δy_vᶜᶠᵃ(i, j, k, grid, v) = @inbounds Δyᶜᶠᵃ(i, j, k, grid) * v[i, j, k]
@inline Δy_cᶜᶜᵃ(i, j, k, grid, c) = @inbounds Δyᶜᶜᵃ(i, j, k, grid) * c[i, j, k]

@inline Δy_uᶠᶜᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Δyᶠᶜᵃ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Δy_vᶜᶠᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Δyᶜᶠᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

@inline Δx_uᶠᶜᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Δxᶠᶜᵃ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Δx_vᶜᶠᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Δxᶜᶠᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

#####
##### Products between fields and grid areas
#####

@inline Ax_ψᵃᵃᶠ(i, j, k, grid, u) = @inbounds Axᵃᵃᶠ(i, j, k, grid) * u[i, j, k]
@inline Ax_ψᵃᵃᶜ(i, j, k, grid, c) = @inbounds Axᵃᵃᶜ(i, j, k, grid) * c[i, j, k]
@inline Ay_ψᵃᵃᶠ(i, j, k, grid, v) = @inbounds Ayᵃᵃᶠ(i, j, k, grid) * v[i, j, k]
@inline Ay_ψᵃᵃᶜ(i, j, k, grid, c) = @inbounds Ayᵃᵃᶜ(i, j, k, grid) * c[i, j, k]
@inline Az_ψᵃᵃᵃ(i, j, k, grid, c) = @inbounds Azᵃᵃᵃ(i, j, k, grid) * c[i, j, k]

@inline Ax_uᶠᶜᶜ(i, j, k, grid, u) = @inbounds Axᶠᶜᶜ(i, j, k, grid) * u[i, j, k]
@inline Ay_vᶜᶠᶜ(i, j, k, grid, v) = @inbounds Ayᶜᶠᶜ(i, j, k, grid) * v[i, j, k]
@inline Az_wᶜᶜᵃ(i, j, k, grid, w) = @inbounds Azᶜᶜᵃ(i, j, k, grid) * w[i, j, k]

@inline Ax_uᶠᶜᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Axᶠᶜᶜ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Ay_vᶜᶠᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Ayᶜᶠᶜ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Az_wᶜᶜᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Azᶜᶜᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

@inline Ax_cᶜᶜᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Axᶜᶜᶜ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Ay_cᶜᶜᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Ayᶜᶜᶜ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Az_cᶜᶜᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Azᶜᶜᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

@inline Ax_ζᶠᶠᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Axᶠᶠᶜ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Ay_ζᶠᶠᶜ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Ayᶠᶠᶜ(i, j, k, grid) * f(i, j, k, grid, args...)

@inline Ax_ηᶠᶜᶠ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Axᶠᶜᶠ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Az_ηᶠᶜᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Azᶠᶜᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

@inline Ay_ξᶜᶠᶠ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Ayᶜᶠᶠ(i, j, k, grid) * f(i, j, k, grid, args...)
@inline Az_ξᶜᶠᵃ(i, j, k, grid, f::F, args...) where F <: Function = @inbounds Azᶜᶠᵃ(i, j, k, grid) * f(i, j, k, grid, args...)

ext: jl


text:
loggers = Dict{String, Logger}()

function getlogger(name::String)
    names = rsplit(name, ".", 2)
    parent = (length(names) > 1) ? getlogger(names[1]) : Logging._root

    haskey(loggers, name) && return loggers[name]
    loggers[name] = Logger(name, Logging.DEBUG, STDOUT, parent)
end

function configure_loggers(args...)
    configure(Logging._root, args...)
    for name in ["circuitscape", "cfg"]
        configure_logger(name, args...)
    end
end

configure_logger(name::String, args...) = configure(getlogger(name, args...))


text:
function getInteractionHamIPicture(
    Ham_0::Array,
    Ham_I::Array,
    t::AbstractFloat,
)::Array
    U_array = evolutionOperator(Ham_0, t)
    return adjoint(U_array) * Ham_I * U_array
end

function getInteractionHamIPicture(
    Ham_0::Operator,
    Ham_I::Operator,
    t::AbstractFloat,
)::Operator
    U_op = evolutionOperator(Ham_0, t)
    return U_op' * Ham_I * U_op
end

function getInteractionHamIPictureA(
    Ham_0::Operator,
    Ham_I::Operator,
    t::AbstractFloat,
)::Array
    U_op = evolutionOperator(Ham_0, t)
    return (U_op' * Ham_I * U_op).data
end

function getInteractionHamIPictureA(
    Ham_int::Array,
    H_lambda::Array,
    H_S::Array,
    H_Sinv::Array,
    t::AbstractFloat,
)
    U_op = evolutionOperatorA(H_lambda, H_S, H_Sinv, t)
    U_opd = adjoint(U_op)
    return U_opd * Ham_int * U_op
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
using Test
using LinearAlgebra
using QuantumPropagators

@testset "random Hermitian state" begin

    # Input
    N = 1000

    X = rand(ComplexF64, (N, N))
    H = Hermitian(X)

    dt = 0.5

    Ψ₀ = rand(ComplexF64, N)
    Ψ₀ ./= norm(Ψ₀)

    @test norm(Ψ₀) ≈ 1

    # Expected
    U = exp(-im * H * dt)
    Ψ_out_expected = U * Ψ₀
    @test norm(Ψ_out_expected) ≈ 1

    Ψ = copy(Ψ₀)
    wrk = ExpPropWrk(Ψ₀)
    expprop!(Ψ, H, dt, wrk)
    Ψ_out = copy(Ψ)
    @test norm(Ψ_out) ≈ 1

    # Comparison (should be exact to machine precision)
    @test norm(Ψ_out - Ψ_out_expected) ≈ 0
end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets
using ECOS

# -----------------------------------------------------------------------------
# Trivial tests
# -----------------------------------------------------------------------------

@testset "Guessing game --------- " begin
	@polyvar x y
	p  = (x - y)^2
	Sx = @set 1 - x^2 >= 0
	Sy = @set 1 - y^2 >= 0
	opt = lazy_relax(ECOS.Optimizer, SDDPBridge)

	expected_measure_x = AtomicMeasure(variables(Sx), [
		WeightedDiracMeasure([1], 0.5),
		WeightedDiracMeasure([-1], 0.5)])
	expected_measure_y = AtomicMeasure(variables(Sy), [
		WeightedDiracMeasure([0], 1)])
	expected_value = 1

	pvalue, pstrat_x, pstrat_y = solve_game(+p, Sx, Sy, opt, iteration=0x1)
	dvalue, dstrat_y, dstrat_x = solve_game(-p, Sy, Sx, opt, iteration=0x1)

	eps = 1e-3
	@test isapprox(expected_value, pvalue, atol=eps)
	@test isapprox(pvalue, -dvalue, atol=eps)
	@test isapprox(expected_measure_x, dstrat_x, atol=eps)
	@test isapprox(expected_measure_y, pstrat_y, atol=eps)
end


text:
lock(reentrantLock) do
    # do whatever
end


text:
# This list was obtained through AWK with Knitro 11.1 (and then pared down) with the following two commands:
# grep "#define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$2,$3) }'
# grep "#  define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$3,$4) }'
# For version up until 10.3, use following commands :
# grep "#define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$2,$3) }'
# grep "#  define" knitro.h | awk '{ printf("const %s = Int32(%s)\n",$3,$4) }'

const KNTRUE = Int32(1)
const KNFALSE = Int32(0)
const KN_INFINITY = Float64(1.0e20)
const KN_PARAMTYPE_INTEGER = Int32(0)
const KN_PARAMTYPE_FLOAT = Int32(1)
const KN_PARAMTYPE_STRING = Int32(2)
const KN_OBJGOAL_MINIMIZE = Int32(0)
const KN_OBJGOAL_MAXIMIZE = Int32(1)
const KN_OBJTYPE_CONSTANT = Int32(-1)
const KN_OBJTYPE_GENERAL = Int32(0)
const KN_OBJTYPE_LINEAR = Int32(1)
const KN_OBJTYPE_QUADRATIC = Int32(2)
const KN_CONTYPE_CONSTANT = Int32(-1)
const KN_CONTYPE_GENERAL = Int32(0)
const KN_CONTYPE_LINEAR = Int32(1)
const KN_CONTYPE_QUADRATIC = Int32(2)
const KN_CONTYPE_CONIC = Int32(3)
const KN_RSDTYPE_CONSTANT = Int32(-1)
const KN_RSDTYPE_GENERAL = Int32(0)
const KN_RSDTYPE_LINEAR = Int32(1)
const KN_CCTYPE_VARVAR = Int32(0)
const KN_CCTYPE_VARCON = Int32(1)
const KN_CCTYPE_CONCON = Int32(2)
const KN_VARTYPE_CONTINUOUS = Int32(0)
const KN_VARTYPE_INTEGER = Int32(1)
const KN_VARTYPE_BINARY = Int32(2)
const KN_VAR_LINEAR = Int32(1)
const KN_OBJ_CONVEX = Int32(1)
const KN_OBJ_CONCAVE = Int32(2)
const KN_OBJ_CONTINUOUS = Int32(4)
const KN_OBJ_DIFFERENTIABLE = Int32(8)
const KN_OBJ_TWICE_DIFFERENTIABLE = Int32(16)
const KN_OBJ_NOISY = Int32(32)
const KN_OBJ_NONDETERMINISTIC = Int32(64)
const KN_CON_CONVEX = Int32(1)
const KN_CON_CONCAVE = Int32(2)
const KN_CON_CONTINUOUS = Int32(4)
const KN_CON_DIFFERENTIABLE = Int32(8)
const KN_CON_TWICE_DIFFERENTIABLE = Int32(16)
const KN_CON_NOISY = Int32(32)
const KN_CON_NONDETERMINISTIC = Int32(64)
const KN_DENSE = Int32(-1)
const KN_DENSE_ROWMAJOR = Int32(-2)
const KN_DENSE_COLMAJOR = Int32(-3)
const KN_RC_EVALFC = Int32(1)
const KN_RC_EVALGA = Int32(2)
const KN_RC_EVALH = Int32(3)
const KN_RC_EVALHV = Int32(7)
const KN_RC_EVALH_NO_F = Int32(8)
const KN_RC_EVALHV_NO_F = Int32(9)
const KN_RC_EVALR = Int32(10)
const KN_RC_EVALRJ = Int32(11)
const KN_RC_EVALFCGA = Int32(12)
const KN_RC_OPTIMAL_OR_SATISFACTORY = Int32(0)
const KN_RC_OPTIMAL = Int32(0)
const KN_RC_NEAR_OPT = Int32(-100)
const KN_RC_FEAS_XTOL = Int32(-101)
const KN_RC_FEAS_NO_IMPROVE = Int32(-102)
const KN_RC_FEAS_FTOL = Int32(-103)
const KN_RC_INFEASIBLE = Int32(-200)
const KN_RC_INFEAS_XTOL = Int32(-201)
const KN_RC_INFEAS_NO_IMPROVE = Int32(-202)
const KN_RC_INFEAS_MULTISTART = Int32(-203)
const KN_RC_INFEAS_CON_BOUNDS = Int32(-204)
const KN_RC_INFEAS_VAR_BOUNDS = Int32(-205)
const KN_RC_UNBOUNDED = Int32(-300)
if KNITRO_VERSION > v"12.0"
    const KN_RC_UNBOUNDED_OR_INFEAS = Int32(-301)
end
const KN_RC_ITER_LIMIT_FEAS = Int32(-400)
const KN_RC_TIME_LIMIT_FEAS = Int32(-401)
const KN_RC_FEVAL_LIMIT_FEAS = Int32(-402)
const KN_RC_MIP_EXH_FEAS = Int32(-403)
const KN_RC_MIP_TERM_FEAS = Int32(-404)
const KN_RC_MIP_SOLVE_LIMIT_FEAS = Int32(-405)
const KN_RC_MIP_NODE_LIMIT_FEAS = Int32(-406)
const KN_RC_ITER_LIMIT_INFEAS = Int32(-410)
const KN_RC_TIME_LIMIT_INFEAS = Int32(-411)
const KN_RC_FEVAL_LIMIT_INFEAS = Int32(-412)
const KN_RC_MIP_EXH_INFEAS = Int32(-413)
const KN_RC_MIP_SOLVE_LIMIT_INFEAS = Int32(-415)
const KN_RC_MIP_NODE_LIMIT_INFEAS = Int32(-416)
const KN_RC_CALLBACK_ERR = Int32(-500)
const KN_RC_LP_SOLVER_ERR = Int32(-501)
const KN_RC_EVAL_ERR = Int32(-502)
const KN_RC_OUT_OF_MEMORY = Int32(-503)
const KN_RC_USER_TERMINATION = Int32(-504)
const KN_RC_OPEN_FILE_ERR = Int32(-505)
const KN_RC_BAD_N_OR_F = Int32(-506)
const KN_RC_BAD_CONSTRAINT = Int32(-507)
const KN_RC_BAD_JACOBIAN = Int32(-508)
const KN_RC_BAD_HESSIAN = Int32(-509)
const KN_RC_BAD_CON_INDEX = Int32(-510)
const KN_RC_BAD_JAC_INDEX = Int32(-511)
const KN_RC_BAD_HESS_INDEX = Int32(-512)
const KN_RC_BAD_CON_BOUNDS = Int32(-513)
const KN_RC_BAD_VAR_BOUNDS = Int32(-514)
const KN_RC_ILLEGAL_CALL = Int32(-515)
const KN_RC_BAD_KCPTR = Int32(-516)
const KN_RC_NULL_POINTER = Int32(-517)
const KN_RC_BAD_INIT_VALUE = Int32(-518)
const KN_RC_LICENSE_ERROR = Int32(-520)
const KN_RC_BAD_PARAMINPUT = Int32(-521)
const KN_RC_LINEAR_SOLVER_ERR = Int32(-522)
const KN_RC_DERIV_CHECK_FAILED = Int32(-523)
const KN_RC_DERIV_CHECK_TERMINATE = Int32(-524)
const KN_RC_OVERFLOW_ERR = Int32(-525)
const KN_RC_BAD_SIZE = Int32(-526)
const KN_RC_BAD_VARIABLE = Int32(-527)
const KN_RC_BAD_VAR_INDEX = Int32(-528)
const KN_RC_BAD_OBJECTIVE = Int32(-529)
const KN_RC_BAD_OBJ_INDEX = Int32(-530)
const KN_RC_BAD_RESIDUAL = Int32(-531)
const KN_RC_BAD_RSD_INDEX = Int32(-532)
const KN_RC_INTERNAL_ERROR = Int32(-600)
const KN_PARAM_NEWPOINT = Int32(1001)
const KN_PARAM_HONORBNDS = Int32(1002)
const KN_PARAM_ALGORITHM = Int32(1003)
const KN_PARAM_ALG = Int32(1003)
const KN_PARAM_BAR_MURULE = Int32(1004)
const KN_PARAM_BAR_FEASIBLE = Int32(1006)
const KN_PARAM_GRADOPT = Int32(1007)
const KN_PARAM_HESSOPT = Int32(1008)
const KN_PARAM_BAR_INITPT = Int32(1009)
const KN_PARAM_ACT_LPSOLVER = Int32(1012)
const KN_PARAM_CG_MAXIT = Int32(1013)
const KN_PARAM_MAXIT = Int32(1014)
const KN_PARAM_OUTLEV = Int32(1015)
const KN_PARAM_OUTMODE = Int32(1016)
const KN_PARAM_SCALE = Int32(1017)
const KN_PARAM_SOC = Int32(1019)
const KN_PARAM_DELTA = Int32(1020)
const KN_PARAM_BAR_FEASMODETOL = Int32(1021)
const KN_PARAM_FEASTOL = Int32(1022)
const KN_PARAM_FEASTOLABS = Int32(1023)
const KN_PARAM_MAXTIMECPU = Int32(1024)
const KN_PARAM_BAR_INITMU = Int32(1025)
const KN_PARAM_OBJRANGE = Int32(1026)
const KN_PARAM_OPTTOL = Int32(1027)
const KN_PARAM_OPTTOLABS = Int32(1028)
const KN_PARAM_LINSOLVER_PIVOTTOL = Int32(1029)
const KN_PARAM_XTOL = Int32(1030)
const KN_PARAM_DEBUG = Int32(1031)
const KN_PARAM_MULTISTART = Int32(1033)
const KN_PARAM_MSENABLE = Int32(1033)
const KN_PARAM_MSMAXSOLVES = Int32(1034)
const KN_PARAM_MSMAXBNDRANGE = Int32(1035)
const KN_PARAM_MSMAXTIMECPU = Int32(1036)
const KN_PARAM_MSMAXTIMEREAL = Int32(1037)
const KN_PARAM_LMSIZE = Int32(1038)
const KN_PARAM_BAR_MAXCROSSIT = Int32(1039)
const KN_PARAM_MAXTIMEREAL = Int32(1040)
const KN_PARAM_CG_PRECOND = Int32(1041)
const KN_PARAM_BLASOPTION = Int32(1042)
const KN_PARAM_BAR_MAXREFACTOR = Int32(1043)
const KN_PARAM_LINESEARCH_MAXTRIALS = Int32(1044)
const KN_PARAM_BLASOPTIONLIB = Int32(1045)
const KN_PARAM_OUTAPPEND = Int32(1046)
const KN_PARAM_OUTDIR = Int32(1047)
const KN_PARAM_CPLEXLIB = Int32(1048)
const KN_PARAM_BAR_PENRULE = Int32(1049)
const KN_PARAM_BAR_PENCONS = Int32(1050)
const KN_PARAM_MSNUMTOSAVE = Int32(1051)
const KN_PARAM_MSSAVETOL = Int32(1052)
const KN_PARAM_PRESOLVEDEBUG = Int32(1053)
const KN_PARAM_MSTERMINATE = Int32(1054)
const KN_PARAM_MSSTARTPTRANGE = Int32(1055)
const KN_PARAM_INFEASTOL = Int32(1056)
const KN_PARAM_LINSOLVER = Int32(1057)
const KN_PARAM_BAR_DIRECTINTERVAL = Int32(1058)
const KN_PARAM_PRESOLVE = Int32(1059)
const KN_PARAM_PRESOLVE_TOL = Int32(1060)
const KN_PARAM_BAR_SWITCHRULE = Int32(1061)
const KN_PARAM_HESSIAN_NO_F = Int32(1062)
const KN_PARAM_MA_TERMINATE = Int32(1063)
const KN_PARAM_MA_MAXTIMECPU = Int32(1064)
const KN_PARAM_MA_MAXTIMEREAL = Int32(1065)
const KN_PARAM_MSSEED = Int32(1066)
const KN_PARAM_MA_OUTSUB = Int32(1067)
const KN_PARAM_MS_OUTSUB = Int32(1068)
const KN_PARAM_XPRESSLIB = Int32(1069)
const KN_PARAM_TUNER = Int32(1070)
const KN_PARAM_TUNER_OPTIONSFILE = Int32(1071)
const KN_PARAM_TUNER_MAXTIMECPU = Int32(1072)
const KN_PARAM_TUNER_MAXTIMEREAL = Int32(1073)
const KN_PARAM_TUNER_OUTSUB = Int32(1074)
const KN_PARAM_TUNER_TERMINATE = Int32(1075)
const KN_PARAM_LINSOLVER_OOC = Int32(1076)
const KN_PARAM_BAR_RELAXCONS = Int32(1077)
const KN_PARAM_MSDETERMINISTIC = Int32(1078)
const KN_PARAM_BAR_REFINEMENT = Int32(1079)
const KN_PARAM_DERIVCHECK = Int32(1080)
const KN_PARAM_DERIVCHECK_TYPE = Int32(1081)
const KN_PARAM_DERIVCHECK_TOL = Int32(1082)
const KN_PARAM_LINSOLVER_INEXACT = Int32(1083)
const KN_PARAM_LINSOLVER_INEXACTTOL = Int32(1084)
const KN_PARAM_MAXFEVALS = Int32(1085)
const KN_PARAM_FSTOPVAL = Int32(1086)
const KN_PARAM_DATACHECK = Int32(1087)
const KN_PARAM_DERIVCHECK_TERMINATE = Int32(1088)
const KN_PARAM_BAR_WATCHDOG = Int32(1089)
const KN_PARAM_FTOL = Int32(1090)
const KN_PARAM_FTOL_ITERS = Int32(1091)
const KN_PARAM_ACT_QPALG = Int32(1092)
const KN_PARAM_BAR_INITPI_MPEC = Int32(1093)
const KN_PARAM_XTOL_ITERS = Int32(1094)
const KN_PARAM_LINESEARCH = Int32(1095)
const KN_PARAM_OUT_CSVINFO = Int32(1096)
const KN_PARAM_INITPENALTY = Int32(1097)
const KN_PARAM_ACT_LPFEASTOL = Int32(1098)
const KN_PARAM_CG_STOPTOL = Int32(1099)
const KN_PARAM_RESTARTS = Int32(1100)
const KN_PARAM_RESTARTS_MAXIT = Int32(1101)
const KN_PARAM_BAR_SLACKBOUNDPUSH = Int32(1102)
const KN_PARAM_CG_PMEM = Int32(1103)
const KN_PARAM_BAR_SWITCHOBJ = Int32(1104)
const KN_PARAM_OUTNAME = Int32(1105)
const KN_PARAM_OUT_CSVNAME = Int32(1106)
const KN_PARAM_ACT_PARAMETRIC = Int32(1107)
const KN_PARAM_ACT_LPDUMPMPS = Int32(1108)
const KN_PARAM_ACT_LPALG = Int32(1109)
const KN_PARAM_ACT_LPPRESOLVE = Int32(1110)
const KN_PARAM_ACT_LPPENALTY = Int32(1111)
const KN_PARAM_BNDRANGE = Int32(1112)
const KN_PARAM_BAR_CONIC_ENABLE = Int32(1113)
const KN_PARAM_CONVEX = Int32(1114)
const KN_PARAM_OUT_HINTS = Int32(1115)
const KN_PARAM_EVAL_FCGA = Int32(1116)
const KN_PARAM_MIP_METHOD = Int32(2001)
const KN_PARAM_MIP_BRANCHRULE = Int32(2002)
const KN_PARAM_MIP_SELECTRULE = Int32(2003)
const KN_PARAM_MIP_INTGAPABS = Int32(2004)
const KN_PARAM_MIP_INTGAPREL = Int32(2005)
const KN_PARAM_MIP_MAXTIMECPU = Int32(2006)
const KN_PARAM_MIP_MAXTIMEREAL = Int32(2007)
const KN_PARAM_MIP_MAXSOLVES = Int32(2008)
const KN_PARAM_MIP_INTEGERTOL = Int32(2009)
const KN_PARAM_MIP_OUTLEVEL = Int32(2010)
const KN_PARAM_MIP_OUTINTERVAL = Int32(2011)
const KN_PARAM_MIP_OUTSUB = Int32(2012)
const KN_PARAM_MIP_DEBUG = Int32(2013)
const KN_PARAM_MIP_IMPLICATNS = Int32(2014)
const KN_PARAM_MIP_GUB_BRANCH = Int32(2015)
const KN_PARAM_MIP_KNAPSACK = Int32(2016)
const KN_PARAM_MIP_ROUNDING = Int32(2017)
const KN_PARAM_MIP_ROOTALG = Int32(2018)
const KN_PARAM_MIP_LPALG = Int32(2019)
const KN_PARAM_MIP_TERMINATE = Int32(2020)
const KN_PARAM_MIP_MAXNODES = Int32(2021)
const KN_PARAM_MIP_HEURISTIC = Int32(2022)
const KN_PARAM_MIP_HEUR_MAXIT = Int32(2023)
const KN_PARAM_MIP_HEUR_MAXTIMECPU = Int32(2024)
const KN_PARAM_MIP_HEUR_MAXTIMEREAL = Int32(2025)
const KN_PARAM_MIP_PSEUDOINIT = Int32(2026)
const KN_PARAM_MIP_STRONG_MAXIT = Int32(2027)
const KN_PARAM_MIP_STRONG_CANDLIM = Int32(2028)
const KN_PARAM_MIP_STRONG_LEVEL = Int32(2029)
const KN_PARAM_MIP_INTVAR_STRATEGY = Int32(2030)
const KN_PARAM_MIP_RELAXABLE = Int32(2031)
const KN_PARAM_MIP_NODEALG = Int32(2032)
const KN_PARAM_MIP_HEUR_TERMINATE = Int32(2033)
const KN_PARAM_MIP_SELECTDIR = Int32(2034)
const KN_PARAM_PAR_NUMTHREADS = Int32(3001)
const KN_PARAM_PAR_CONCURRENT_EVALS = Int32(3002)
const KN_PARAM_PAR_BLASNUMTHREADS = Int32(3003)
const KN_PARAM_PAR_LSNUMTHREADS = Int32(3004)
const KN_PARAM_PAR_MSNUMTHREADS = Int32(3005)
const KN_NEWPOINT_NONE = Int32(0)
const KN_NEWPOINT_SAVEONE = Int32(1)
const KN_NEWPOINT_SAVEALL = Int32(2)
const KN_HONORBNDS_NO = Int32(0)
const KN_HONORBNDS_ALWAYS = Int32(1)
const KN_HONORBNDS_INITPT = Int32(2)
const KN_ALG_AUTOMATIC = Int32(0)
const KN_ALG_AUTO = Int32(0)
const KN_ALG_BAR_DIRECT = Int32(1)
const KN_ALG_BAR_CG = Int32(2)
const KN_ALG_ACT_CG = Int32(3)
const KN_ALG_ACT_SQP = Int32(4)
const KN_ALG_MULTI = Int32(5)
const KN_BAR_MURULE_AUTOMATIC = Int32(0)
const KN_BAR_MURULE_AUTO = Int32(0)
const KN_BAR_MURULE_MONOTONE = Int32(1)
const KN_BAR_MURULE_ADAPTIVE = Int32(2)
const KN_BAR_MURULE_PROBING = Int32(3)
const KN_BAR_MURULE_DAMPMPC = Int32(4)
const KN_BAR_MURULE_FULLMPC = Int32(5)
const KN_BAR_MURULE_QUALITY = Int32(6)
const KN_BAR_FEASIBLE_NO = Int32(0)
const KN_BAR_FEASIBLE_STAY = Int32(1)
const KN_BAR_FEASIBLE_GET = Int32(2)
const KN_BAR_FEASIBLE_GET_STAY = Int32(3)
const KN_GRADOPT_EXACT = Int32(1)
const KN_GRADOPT_FORWARD = Int32(2)
const KN_GRADOPT_CENTRAL = Int32(3)
const KN_GRADOPT_USER_FORWARD = Int32(4)
const KN_GRADOPT_USER_CENTRAL = Int32(5)
const KN_HESSOPT_AUTO = Int32(0)
const KN_HESSOPT_EXACT = Int32(1)
const KN_HESSOPT_BFGS = Int32(2)
const KN_HESSOPT_SR1 = Int32(3)
const KN_HESSOPT_PRODUCT_FINDIFF = Int32(4)
const KN_HESSOPT_PRODUCT = Int32(5)
const KN_HESSOPT_LBFGS = Int32(6)
const KN_HESSOPT_GAUSS_NEWTON = Int32(7)
const KN_BAR_INITPT_AUTO = Int32(0)
const KN_BAR_INITPT_CONVEX = Int32(1)
const KN_BAR_INITPT_NEARBND = Int32(2)
const KN_BAR_INITPT_CENTRAL = Int32(3)
const KN_ACT_LPSOLVER_INTERNAL = Int32(1)
const KN_ACT_LPSOLVER_CPLEX = Int32(2)
const KN_ACT_LPSOLVER_XPRESS = Int32(3)
const KN_OUTLEV_NONE = Int32(0)
const KN_OUTLEV_SUMMARY = Int32(1)
const KN_OUTLEV_ITER_10 = Int32(2)
const KN_OUTLEV_ITER = Int32(3)
const KN_OUTLEV_ITER_VERBOSE = Int32(4)
const KN_OUTLEV_ITER_X = Int32(5)
const KN_OUTLEV_ALL = Int32(6)
const KN_OUTMODE_SCREEN = Int32(0)
const KN_OUTMODE_FILE = Int32(1)
const KN_OUTMODE_BOTH = Int32(2)
const KN_SCALE_NEVER = Int32(0)
const KN_SCALE_NO = Int32(0)
const KN_SCALE_USER_INTERNAL = Int32(1)
const KN_SCALE_USER_NONE = Int32(2)
const KN_SCALE_INTERNAL = Int32(3)
const KN_SOC_NO = Int32(0)
const KN_SOC_MAYBE = Int32(1)
const KN_SOC_YES = Int32(2)
const KN_DEBUG_NONE = Int32(0)
const KN_DEBUG_PROBLEM = Int32(1)
const KN_DEBUG_EXECUTION = Int32(2)
const KN_MULTISTART_NO = Int32(0)
const KN_MULTISTART_YES = Int32(1)
const KN_CG_PRECOND_NONE = Int32(0)
const KN_CG_PRECOND_CHOL = Int32(1)
const KN_BLASOPTION_KNITRO = Int32(0)
const KN_BLASOPTION_INTEL = Int32(1)
const KN_BLASOPTION_DYNAMIC = Int32(2)
const KN_OUTAPPEND_NO = Int32(0)
const KN_OUTAPPEND_YES = Int32(1)
const KN_BAR_PENRULE_AUTO = Int32(0)
const KN_BAR_PENRULE_SINGLE = Int32(1)
const KN_BAR_PENRULE_FLEX = Int32(2)
const KN_BAR_PENCONS_AUTO = Int32(0)
const KN_BAR_PENCONS_NONE = Int32(1)
const KN_BAR_PENCONS_ALL = Int32(2)
const KN_BAR_PENCONS_EQUALITIES = Int32(3)
const KN_BAR_PENCONS_INFEAS = Int32(4)
const KN_PRESOLVEDBG_NONE = Int32(0)
const KN_PRESOLVEDBG_BASIC = Int32(1)
const KN_PRESOLVEDBG_VERBOSE = Int32(2)
const KN_MSTERMINATE_MAXSOLVES = Int32(0)
const KN_MSTERMINATE_OPTIMAL = Int32(1)
const KN_MSTERMINATE_FEASIBLE = Int32(2)
const KN_MSTERMINATE_ANY = Int32(3)
const KN_LINSOLVER_AUTO = Int32(0)
const KN_LINSOLVER_INTERNAL = Int32(1)
const KN_LINSOLVER_HYBRID = Int32(2)
const KN_LINSOLVER_DENSEQR = Int32(3)
const KN_LINSOLVER_MA27 = Int32(4)
const KN_LINSOLVER_MA57 = Int32(5)
const KN_LINSOLVER_MKLPARDISO = Int32(6)
const KN_LINSOLVER_MA97 = Int32(7)
const KN_LINSOLVER_MA86 = Int32(8)
const KN_PRESOLVE_NONE = Int32(0)
const KN_PRESOLVE_BASIC = Int32(1)
const KN_PRESOLVE_ADVANCED = Int32(2)
if KNITRO_VERSION >= v"12.0"
    const KN_PRESOLVE_NO = Int32(0)
    const KN_PRESOLVE_YES = Int32(1)
    const KN_PRESOLVE_ADVANCED = Int32(2)
end
const KN_BAR_SWITCHRULE_AUTO = Int32(0)
const KN_BAR_SWITCHRULE_NEVER = Int32(1)
const KN_BAR_SWITCHRULE_MODERATE = Int32(2)
const KN_BAR_SWITCHRULE_AGGRESSIVE = Int32(3)
const KN_HESSIAN_NO_F_FORBID = Int32(0)
const KN_HESSIAN_NO_F_ALLOW = Int32(1)
const KN_MA_TERMINATE_ALL = Int32(0)
const KN_MA_TERMINATE_OPTIMAL = Int32(1)
const KN_MA_TERMINATE_FEASIBLE = Int32(2)
const KN_MA_TERMINATE_ANY = Int32(3)
const KN_MA_OUTSUB_NONE = Int32(0)
const KN_MA_OUTSUB_YES = Int32(1)
const KN_MS_OUTSUB_NONE = Int32(0)
const KN_MS_OUTSUB_YES = Int32(1)
const KN_TUNER_OFF = Int32(0)
const KN_TUNER_ON = Int32(1)
const KN_TUNER_OUTSUB_NONE = Int32(0)
const KN_TUNER_OUTSUB_SUMMARY = Int32(1)
const KN_TUNER_OUTSUB_ALL = Int32(2)
const KN_TUNER_TERMINATE_ALL = Int32(0)
const KN_TUNER_TERMINATE_OPTIMAL = Int32(1)
const KN_TUNER_TERMINATE_FEASIBLE = Int32(2)
const KN_TUNER_TERMINATE_ANY = Int32(3)
const KN_LINSOLVER_OOC_NO = Int32(0)
const KN_LINSOLVER_OOC_MAYBE = Int32(1)
const KN_LINSOLVER_OOC_YES = Int32(2)
const KN_BAR_RELAXCONS_NONE = Int32(0)
const KN_BAR_RELAXCONS_EQS = Int32(1)
const KN_BAR_RELAXCONS_INEQS = Int32(2)
const KN_BAR_RELAXCONS_ALL = Int32(3)
const KN_MSDETERMINISTIC_NO = Int32(0)
const KN_MSDETERMINISTIC_YES = Int32(1)
const KN_BAR_REFINEMENT_NO = Int32(0)
const KN_BAR_REFINEMENT_YES = Int32(1)
const KN_DERIVCHECK_NONE = Int32(0)
const KN_DERIVCHECK_FIRST = Int32(1)
const KN_DERIVCHECK_SECOND = Int32(2)
const KN_DERIVCHECK_ALL = Int32(3)
const KN_DERIVCHECK_FORWARD = Int32(1)
const KN_DERIVCHECK_CENTRAL = Int32(2)
const KN_LINSOLVER_INEXACT_NO = Int32(0)
const KN_LINSOLVER_INEXACT_YES = Int32(1)
const KN_DATACHECK_NO = Int32(0)
const KN_DATACHECK_YES = Int32(1)
const KN_DERIVCHECK_STOPERROR = Int32(1)
const KN_DERIVCHECK_STOPALWAYS = Int32(2)
const KN_BAR_WATCHDOG_NO = Int32(0)
const KN_BAR_WATCHDOG_YES = Int32(1)
const KN_ACT_QPALG_AUTO = Int32(0)
const KN_ACT_QPALG_BAR_DIRECT = Int32(1)
const KN_ACT_QPALG_BAR_CG = Int32(2)
const KN_ACT_QPALG_ACT_CG = Int32(3)
const KN_LINESEARCH_AUTO = Int32(0)
const KN_LINESEARCH_BACKTRACK = Int32(1)
const KN_LINESEARCH_INTERPOLATE = Int32(2)
const KN_OUT_CSVINFO_NO = Int32(0)
const KN_OUT_CSVINFO_YES = Int32(1)
const KN_BAR_SWITCHOBJ_NONE = Int32(0)
const KN_BAR_SWITCHOBJ_SCALARPROX = Int32(1)
const KN_BAR_SWITCHOBJ_DIAGPROX = Int32(2)
const KN_ACT_PARAMETRIC_NO = Int32(0)
const KN_ACT_PARAMETRIC_MAYBE = Int32(1)
const KN_ACT_PARAMETRIC_YES = Int32(2)
const KN_ACT_LPDUMPMPS_NO = Int32(0)
const KN_ACT_LPDUMPMPS_YES = Int32(1)
const KN_ACT_LPALG_DEFAULT = Int32(0)
const KN_ACT_LPALG_PRIMAL = Int32(1)
const KN_ACT_LPALG_DUAL = Int32(2)
const KN_ACT_LPALG_BARRIER = Int32(3)
const KN_ACT_LPPRESOLVE_OFF = Int32(0)
const KN_ACT_LPPRESOLVE_ON = Int32(1)
const KN_ACT_LPPENALTY_ALL = Int32(1)
const KN_ACT_LPPENALTY_NONLINEAR = Int32(2)
const KN_ACT_LPPENALTY_DYNAMIC = Int32(3)
const KN_BAR_CONIC_ENABLE_NONE = Int32(0)
const KN_BAR_CONIC_ENABLE_SOC = Int32(1)
const KN_CONVEX_NO = Int32(0)
const KN_CONVEX_YES = Int32(1)
const KN_OUT_HINTS_NO = Int32(0)
const KN_OUT_HINTS_YES = Int32(1)
const KN_EVAL_FCGA_NO = Int32(0)
const KN_EVAL_FCGA_YES = Int32(1)
const KN_MIP_METHOD_AUTO = Int32(0)
const KN_MIP_METHOD_BB = Int32(1)
const KN_MIP_METHOD_HQG = Int32(2)
const KN_MIP_METHOD_MISQP = Int32(3)
const KN_MIP_BRANCH_AUTO = Int32(0)
const KN_MIP_BRANCH_MOSTFRAC = Int32(1)
const KN_MIP_BRANCH_PSEUDOCOST = Int32(2)
const KN_MIP_BRANCH_STRONG = Int32(3)
const KN_MIP_SEL_AUTO = Int32(0)
const KN_MIP_SEL_DEPTHFIRST = Int32(1)
const KN_MIP_SEL_BESTBOUND = Int32(2)
const KN_MIP_SEL_COMBO_1 = Int32(3)
const KN_MIP_OUTLEVEL_NONE = Int32(0)
const KN_MIP_OUTLEVEL_ITERS = Int32(1)
const KN_MIP_OUTLEVEL_ITERSTIME = Int32(2)
const KN_MIP_OUTLEVEL_ROOT = Int32(3)
const KN_MIP_OUTSUB_NONE = Int32(0)
const KN_MIP_OUTSUB_YES = Int32(1)
const KN_MIP_OUTSUB_YESPROB = Int32(2)
const KN_MIP_DEBUG_NONE = Int32(0)
const KN_MIP_DEBUG_ALL = Int32(1)
const KN_MIP_IMPLICATNS_NO = Int32(0)
const KN_MIP_IMPLICATNS_YES = Int32(1)
const KN_MIP_GUB_BRANCH_NO = Int32(0)
const KN_MIP_GUB_BRANCH_YES = Int32(1)
const KN_MIP_KNAPSACK_NO = Int32(0)
const KN_MIP_KNAPSACK_NONE = Int32(0)
const KN_MIP_KNAPSACK_INEQ = Int32(1)
const KN_MIP_KNAPSACK_INEQ_EQ = Int32(2)
const KN_MIP_ROUND_AUTO = Int32(0)
const KN_MIP_ROUND_NONE = Int32(1)
const KN_MIP_ROUND_HEURISTIC = Int32(2)
const KN_MIP_ROUND_NLP_SOME = Int32(3)
const KN_MIP_ROUND_NLP_ALWAYS = Int32(4)
const KN_MIP_ROOTALG_AUTO = Int32(0)
const KN_MIP_ROOTALG_BAR_DIRECT = Int32(1)
const KN_MIP_ROOTALG_BAR_CG = Int32(2)
const KN_MIP_ROOTALG_ACT_CG = Int32(3)
const KN_MIP_ROOTALG_ACT_SQP = Int32(4)
const KN_MIP_ROOTALG_MULTI = Int32(5)
const KN_MIP_LPALG_AUTO = Int32(0)
const KN_MIP_LPALG_BAR_DIRECT = Int32(1)
const KN_MIP_LPALG_BAR_CG = Int32(2)
const KN_MIP_LPALG_ACT_CG = Int32(3)
const KN_MIP_TERMINATE_OPTIMAL = Int32(0)
const KN_MIP_TERMINATE_FEASIBLE = Int32(1)
const KN_MIP_HEURISTIC_AUTO = Int32(0)
const KN_MIP_HEURISTIC_NONE = Int32(1)
const KN_MIP_HEURISTIC_FEASPUMP = Int32(2)
const KN_MIP_HEURISTIC_MPEC = Int32(3)
const KN_MIP_PSEUDOINIT_AUTO = Int32(0)
const KN_MIP_PSEUDOINIT_AVE = Int32(1)
const KN_MIP_PSEUDOINIT_STRONG = Int32(2)
const KN_MIP_INTVAR_STRATEGY_NONE = Int32(0)
const KN_MIP_INTVAR_STRATEGY_RELAX = Int32(1)
const KN_MIP_INTVAR_STRATEGY_MPEC = Int32(2)
const KN_MIP_RELAXABLE_NONE = Int32(0)
const KN_MIP_RELAXABLE_ALL = Int32(1)
const KN_MIP_NODEALG_AUTO = Int32(0)
const KN_MIP_NODEALG_BAR_DIRECT = Int32(1)
const KN_MIP_NODEALG_BAR_CG = Int32(2)
const KN_MIP_NODEALG_ACT_CG = Int32(3)
const KN_MIP_NODEALG_ACT_SQP = Int32(4)
const KN_MIP_NODEALG_MULTI = Int32(5)
const KN_MIP_HEUR_TERMINATE_FEASIBLE = Int32(1)
const KN_MIP_HEUR_TERMINATE_LIMIT = Int32(2)
const KN_MIP_SELECTDIR_DOWN = Int32(0)
const KN_MIP_SELECTDIR_UP = Int32(1)
const KN_PAR_CONCURRENT_EVALS_NO = Int32(0)
const KN_PAR_CONCURRENT_EVALS_YES = Int32(1)
const KN_PAR_MSNUMTHREADS_AUTO = Int32(0)
if KNITRO_VERSION >= v"12.0"
    # Cuts parameters
    const KN_PARAM_MIP_ZEROHALF = Int32(2036)
    const KN_PARAM_MIP_ZEROHALF_NONE = Int32(0)
    const KN_PARAM_MIP_ZEROHALF_ROOT = Int32(1)
    const KN_PARAM_MIP_ZEROHALF_TREE = Int32(2)
    const KN_PARAM_MIP_ZEROHALF_ALL = Int32(3)
    const KN_PARAM_MIP_MIR = Int32(2037)
    const KN_PARAM_MIP_MIR_NONE = Int32(0)
    const KN_PARAM_MIP_MIR_TREE = Int32(1)
    const KN_PARAM_MIP_CLIQUE = Int32(2038)
    const KN_PARAM_MIP_CLIQUE_NONE = Int32(0)
    const KN_PARAM_MIP_CLIQUE_ROOT = Int32(1)
    const KN_PARAM_MIP_CLIQUE_TREE = Int32(2)
    const KN_PARAM_MIP_CLIQUE_ALL = Int32(3)
    # SOCP parameters
    const KN_PARAM_PAR_CONICNUMTHREADS = Int32(3006)
    # Presolve parameters
    const KN_PARAM_PRESOLVE_PASSES = Int32(1121)
    const KN_PARAM_PRESOLVE_LEVEL = Int32(1122)
    const KN_PARAM_PRESOLVE_LEVEL_AUTO = Int32(-1)
    const KN_PARAM_PRESOLVE_LEVEL_1 = Int32(1)
    const KN_PARAM_PRESOLVE_LEVEL_2 = Int32(2)
end


text:
ENV["GKSwstype"] = "100"

using Printf
using Statistics
using Plots

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary

#####
##### Grid
#####

arch = GPU()
const Ly = 2000kilometers # north-south extent [m]

Ny = 400
Nz = 35

# s = 1.5 # stretching factor
# z_faces(k) = - Lz * (1 - tanh(s * (k - 1) / Nz) / tanh(s))

k_center = collect(1:Nz)
Δz_center = @. 10 * 1.104^(Nz - k_center)
const Lz = sum(Δz_center)
z_faces = vcat([-Lz], -Lz .+ cumsum(Δz_center))
z_faces[Nz+1] = 0

grid = VerticallyStretchedRectilinearGrid(architecture = arch,
                                          topology = (Periodic, Bounded, Bounded),
                                          size = (1, Ny, Nz),
                                          halo = (3, 3, 3),
                                          x = (0, Ly),
                                          y = (0, Ly),
                                          z_faces = z_faces)

@show grid

#=
# We visualize the cell interfaces by plotting the cell height
# as a function of depth,

p = plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
         marker = :circle,
         ylabel = "Depth (m)",
         xlabel = "Vertical spacing (m)",
         legend = nothing)

display(p)
=#

#####
##### Boundary conditions: buoyancy flux
#####

Qᵇ = 1e-8            # buoyancy flux magnitude [m² s⁻³]
y_shutoff = 5/6 * Ly # shutoff location for buoyancy flux [m]
τ = 1e-4             # surface kinematic wind stress [m² s⁻²]
μ = 1 / 100days      # bottom drag damping time-scale [s⁻¹]

@inline buoyancy_flux(x, y, t, p) = ifelse(y < p.y_shutoff, p.Qᵇ * cos(3π * y / p.Ly), 0)

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, parameters=(Ly=grid.Ly, y_shutoff=y_shutoff, Qᵇ=Qᵇ))

#####
##### Boundary conditions: wind stress
#####

@inline u_stress(x, y, t, p) = - p.τ * sin(π * y / p.Ly)

u_stress_bc = FluxBoundaryCondition(u_stress, parameters=(τ=τ, Ly=grid.Ly))

#####
##### Boundary conditions: linear bottom drag
#####

@inline u_drag(x, y, t, u, μ) = - μ * u
@inline v_drag(x, y, t, v, μ) = - μ * v

u_drag_bc = FluxBoundaryCondition(u_drag, field_dependencies=:u, parameters=μ)
v_drag_bc = FluxBoundaryCondition(v_drag, field_dependencies=:v, parameters=μ)

b_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc)
u_bcs = FieldBoundaryConditions(top = u_stress_bc, bottom = u_drag_bc)
v_bcs = FieldBoundaryConditions(bottom = v_drag_bc)

#####
##### Coriolis forces
#####

coriolis = BetaPlane(latitude=-45)

const N²₀ = 1.6e-5   # surface vertical buoyancy gradient [s⁻²]
const h = 1kilometer # decay scale of stable stratification [m]

#####
##### Initial stratification and restoring
#####

@inline b_stratification(z) = N²₀ * h * exp(z / h) / (1 - exp(-Lz / h))

const y_sponge = 19/20 * Ly # southern boundary of sponge layer [m]

## Mask that limits sponge layer to a thin region near the northern boundary
@inline northern_mask(x, y, z) = max(0, y - y_sponge) / (Ly - y_sponge)

## Target and initial buoyancy profile
@inline b_target(x, y, z, t) = b_stratification(z)

b_forcing = Relaxation(target=b_target, mask=northern_mask, rate=1/7days)

using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization
using Oceananigans.Fields: TracerFields, FunctionField

tracers = TracerFields(tuple(:b), arch, grid)

b = tracers.b

const f₀ = coriolis.f₀
const β = coriolis.β
const K₀ = 100

K_func(x, y, z) = K₀ + 3000 * sin(π * y / Ly)
f²_func(x, y, z) = (f₀ + β * y)^2

K = FunctionField{Center, Center, Center}(K_func, grid)
f² = FunctionField{Center, Center, Center}(f²_func, grid)

ν_op = @at (Center, Center, Center) K * f² / ∂z(b)
ν = ComputedField(ν_op)

closure = AnisotropicDiffusivity(νh = 100, νz = 10, κh = 10, κz = 10,
                                 time_discretization = VerticallyImplicitTimeDiscretization())

model = IncompressibleModel(architecture = arch,
                            grid = grid,
                            advection = UpwindBiasedFifthOrder(),
                            buoyancy = BuoyancyTracer(),
                            coriolis = coriolis,
                            closure = closure,
                            tracers = tracers,
                            boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                            auxiliary_fields = (; ν=ν),
                            forcing = (b=b_forcing,))

@show grid === model.grid

bᵢ(x, y, z) = b_stratification(z)
uᵢ(x, y, z) = 0.0 #0.1 * (1 - z / Lz)

set!(model, u=uᵢ, b=bᵢ)

using Oceananigans.Fields: similar_cpu_field
using Oceananigans.Grids: topology, halo_size

cpu_grid(grid::RegularRectilinearGrid) = grid

cpu_grid(grid::VerticallyStretchedRectilinearGrid) =
    VerticallyStretchedRectilinearGrid(architecture = CPU(),
                                       topology = topology(grid),
                                       size = size(grid),
                                       halo = halo_size(grid),
                                       x = (0, grid.Ly),
                                       y = (0, grid.Ly),
                                       z_faces = grid.zᵃᵃᶠ)

function channel_plot(u_device, b_device)

    grid = cpu_grid(u_device.grid)

    u = XFaceField(CPU(), grid)
    b = CenterField(CPU(), grid)

    copyto!(parent(u), parent(u_device))
    copyto!(parent(b), parent(b_device))

    grid = u.grid

    xb, yb, zb = nodes(b)
    xu, yu, zu = nodes(u)

    ui = interior(u)[1, :, :]
    bi = interior(b)[1, :, :]

    umax = maximum(abs, ui)
    ulim = max(0.001, 0.5 * umax)
    ulevels = range(-ulim, ulim, length=31)
    umax > ulim && (ulevels = vcat([-umax], ulevels, [umax]))

    bmin = minimum(bi)
    bmax = maximum(bi)
    blevels = range(bmin, bmax, length=10)

    pl = contourf(yu * 1e-3, zu, ui',
                  xlabel = "y (km)",
                  ylabel = "z (m)",
                  aspectratio = 0.2,
                  linewidth = 0,
                  levels = ulevels,
                  clims = (-ulim, ulim),
                  xlims = (0, grid.Ly) .* 1e-3,
                  ylims = (-grid.Lz, 0),
                  color = :balance)

    contour!(pl, yb * 1e-3, zb, bi', linewidth=1, seriescolor=:gray, levels=blevels, clims=(-ulim, ulim))

    return pl
end

u, v, w = model.velocities
b = model.tracers.b

p = channel_plot(u, b)
display(p)

# # Simulation setup
#
# We set up a simulation with adaptive time-stepping and a simple progress message.

#wizard = TimeStepWizard(cfl=0.2, Δt=10minutes, max_change=1.1, max_Δt=2hours)

wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))
            #prettytime(sim.Δt.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

diffusion_Δt = grid.Δy^2 / closure.νy
@show Δt = min(10minutes, diffusion_Δt)

simulation = Simulation(model, Δt=Δt, stop_time=10years, progress=print_progress, iteration_interval=1000)

u, v, w = model.velocities
b = model.tracers.b

outputs = merge(model.velocities, model.tracers)

#=
simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(100days),
                                                        prefix = "eddying_channel",
                                                        force = true)

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                      schedule = TimeInterval(10day),
                                                      prefix = "eddying_channel",
                                                      field_slicer = nothing,
                                                      force = true)

try
    run!(simulation, pickup=false)
catch err
    showerror(stdout, err)
end
=#

u_timeseries = FieldTimeSeries("eddying_channel.jld2", "u", grid=cpu_grid(grid))
b_timeseries = FieldTimeSeries("eddying_channel.jld2", "b", grid=cpu_grid(grid))

@show b_timeseries

anim = @animate for i in 5:120
    @info "Animating frame $i of $(length(b_timeseries.times))..."
    u = b_timeseries[i]
    b = b_timeseries[i]
    channel_plot(u, b)
end

mp4(anim, "two_dimensional_channel.mp4", fps = 8) # hide


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct ThreePartCost <: OperationalCost
        variable::VariableCost
        fixed::Float64
        start_up::Float64
        shut_down::Float64
    end

Data Structure Operational Cost Data in Three parts fixed, variable cost and start - stop costs.

# Arguments
- `variable::VariableCost`: variable cost
- `fixed::Float64`: fixed cost
- `start_up::Float64`: start-up cost, validation range: `(0, nothing)`, action if invalid: `warn`
- `shut_down::Float64`: shut-down cost, validation range: `(0, nothing)`, action if invalid: `warn`
"""
mutable struct ThreePartCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "fixed cost"
    fixed::Float64
    "start-up cost"
    start_up::Float64
    "shut-down cost"
    shut_down::Float64
end

function ThreePartCost(; variable, fixed, start_up, shut_down, )
    ThreePartCost(variable, fixed, start_up, shut_down, )
end

# Constructor for demo purposes; non-functional.
function ThreePartCost(::Nothing)
    ThreePartCost(;
        variable=VariableCost((0.0, 0.0)),
        fixed=0.0,
        start_up=0.0,
        shut_down=0.0,
    )
end

"""Get [`ThreePartCost`](@ref) `variable`."""
get_variable(value::ThreePartCost) = value.variable
"""Get [`ThreePartCost`](@ref) `fixed`."""
get_fixed(value::ThreePartCost) = value.fixed
"""Get [`ThreePartCost`](@ref) `start_up`."""
get_start_up(value::ThreePartCost) = value.start_up
"""Get [`ThreePartCost`](@ref) `shut_down`."""
get_shut_down(value::ThreePartCost) = value.shut_down

"""Set [`ThreePartCost`](@ref) `variable`."""
set_variable!(value::ThreePartCost, val) = value.variable = val
"""Set [`ThreePartCost`](@ref) `fixed`."""
set_fixed!(value::ThreePartCost, val) = value.fixed = val
"""Set [`ThreePartCost`](@ref) `start_up`."""
set_start_up!(value::ThreePartCost, val) = value.start_up = val
"""Set [`ThreePartCost`](@ref) `shut_down`."""
set_shut_down!(value::ThreePartCost, val) = value.shut_down = val


text:
module MethodsFor

export methodsfor

using DataStructures: DefaultDict
using InteractiveUtils: methodswith
using Crayons: Crayon

struct MethodInfo
    name::Symbol
    def_mod::Module
    fx_mod::Module
end

struct GroupedMethods
    method_dct
end

const method_color = Crayon(foreground=:light_yellow)
const module_color = Crayon(foreground=:light_cyan)
const title_color = Crayon(foreground=:light_green)

function _groupmethods(methods::Array{Method}, q)
    exported = Dict{MethodInfo, Int}()
    method_dct = DefaultDict(() -> DefaultDict(() -> 0))

    for method in methods
        if !isnothing(q) && !occursin(q, string(method.name))
            continue
        end

        sig = Base.unwrap_unionall(method.sig)
        generic_fx_mod = sig.types[1].name.module

        mi = MethodInfo(method.name, method.module, generic_fx_mod)

        is_exported = mi.name in names(mi.def_mod) || mi.fx_mod == Base

        if is_exported
            method_dct[mi.fx_mod][mi.name] += 1
        end
    end

    return GroupedMethods(method_dct)
end

function Base.show(io::IO, ::MIME"text/plain", g::GroupedMethods)
    method_dct = g.method_dct
    colwidth = 45

    group_lengths = Dict((k, length(method_dct[k])) for k in keys(method_dct))
    ordered_keys = sort(collect(keys(group_lengths)), by=x->group_lengths[x]; rev=true)

    cols = [[], []]
    col_len(i) = isempty(cols[i]) ? 0 : sum(group_lengths[x] for x in cols[i])

    for k in ordered_keys
        next_length = group_lengths[k]
        if col_len(1) > col_len(2) + next_length + 2
            push!(cols[2], k)
        else
            push!(cols[1], k)
        end
    end
    coltexts = [[], []]
    for i in 1:2
        for k in cols[i]
            push!(coltexts[i], string(module_color(string(k))) * ":")

            for f in sort(collect(keys(method_dct[k])))
                count = method_dct[k][f]
                push!(coltexts[i], "  " * string(method_color(string(f))) * " with $count methods")
            end
        end
    end

    col_text(i, row) = row > length(coltexts[i]) ? "" : coltexts[i][row]

    for row in 1:max(length(coltexts[1]), length(coltexts[2]))
        one = rpad(col_text(1, row), colwidth)
        two = col_text(2, row)
        println(io, one, two)
    end
end

function _printmodule(types, modules)
    colwidth = 30

    coltexts = [
        [string(module_color("Types:"))],
        [string(module_color("Modules:"))],
    ]
    for typ in types
        push!(coltexts[1], "  " * string(typ))
    end
    for mod in modules
        push!(coltexts[2], "  " * string(mod))
    end

    col_text(i, row) = row > length(coltexts[i]) ? "" : coltexts[i][row]

    for row in 1:max(length(coltexts[1]), length(coltexts[2]))
        extra_pad = row == 1 ? 10 : 0
        one = rpad(col_text(1, row), colwidth + extra_pad)
        two = col_text(2, row)
        println(one, two)
    end
end

"""
    methodsfor(obj; q=nothing, supertypes=true, ret=false)

Pretty print list of exported methods related to `obj`
 - `obj` is a value -> functions that take a  typeof(obj)
 - `obj` is a type -> function that take a `obj`
 - `obj` is a Module -> all public types, modules and functions of `obj`
"""
function methodsfor(::Type{T}; q=nothing, supertypes=true) where {T}
    methods = methodswith(T; supertypes=supertypes)
    grouped = _groupmethods(methods, q)
    return grouped
end
methodsfor(val; kwargs...) = methodsfor(typeof(val); kwargs...)

function methodsfor(m::Module; functions=false, q=nothing)
    types = []
    fs = Method[]
    modules = []

    for name in names(m)
        obj = nothing
        try
            obj = getfield(m, name)
        catch
            continue
        end

        typ = typeof(obj)

        if typ == UnionAll || typ == DataType
            push!(types, name)
        elseif typ == Module
            push!(modules, name)
        else
            append!(fs, methods(obj))
        end
    end

    method_dct = _groupmethods(fs, q)

    _printmodule(types, modules)

    if functions
        println("\nExported Functions:\n")
        _printmethods(method_dct)
    end
end

end


text:
module fixed_parameters
    tick_second_rate = 1
    coord_meter_rate = 1

    # Basic
    window_size = (16, 7)

    radius_with_masks = 0.4*coord_meter_rate
    radius_without_masks = 1.5*coord_meter_rate

    move_random_rate = 0.2
    move_speed = (coord_meter_rate*1.5)/tick_second_rate
    indoor_speed = (coord_meter_rate*0.5)/tick_second_rate
    quarantine_speed = (coord_meter_rate*0.2)/tick_second_rate

    day_tick = 86400*tick_second_rate÷60
    sim_time = 10

    end_infection_time = 14
    death_rate = 0.40
    fake_negative_test = 0.04

    # Advanced
    infection_period = 60÷60
    d = radius_without_masks/2

    lund_size = (400, 90)
    cwh, hwh = 10, 10
    lundbounds = [([0,lund_size[1]], [lund_size[2],lund_size[2]]), ([0,lund_size[1]], [0,0]), ([0,0], [0,lund_size[2]]), ([lund_size[1],lund_size[1]], [0,lund_size[2]])]

    num_in = [16,16,20,12,21,14,2,3,23,7,6,7,11,16,12,13,6,14,11,17,7,14,21,4,4,4]
    num_intern_students = 301
    num_houses = length(num_in)
    num_classrooms = 30
    lundwalls = Array{Tuple}(undef, num_houses+num_classrooms-2)
    start = 0.0
    for i in 1:num_houses-1
        global start
        start += num_in[i]/num_intern_students*lund_size[1]
        lundwalls[i] = ([start,start], [lund_size[2]-hwh,lund_size[2]])
    end
    start = 0.0
    for i in num_houses:num_houses+num_classrooms-2
        global start
        start += 1/num_classrooms*lund_size[1]
        lundwalls[i] = ([start,start], [0,cwh])
    end

    function get_infection_distance(radius_1::Float64, radius_2::Float64)
        return max(radius_1, radius_2)-abs((radius_1-radius_2)/2)
    end

    function get_infection_rate(distance::Float64, infection_distance::Float64)
        return (0.95/infection_distance)*((infection_distance-distance)/(distance+1))
    end

    function get_show_infection_sign_rate(days)
        return 0.00156566*days^3-0.0493434*days^2+0.383939*days
    end

    function get_infection_active_rate(days)
        return exp((days-5.7)/1.1)/(exp((days-5.7)/1.1)+1)
    end

    function get_cure_rate(days)
        return log(end_infection_time-days+1)/log(end_infection_time+1)
    end

    export tick_second_rate, coord_meter_rate, window_size, radius_with_masks, radius_without_masks, move_random_rate, move_speed, indoor_speed, quarantine_speed, day_tick, sim_time, end_infection_time, death_rate, fake_negative_test, infection_period, d, lund_size, lundbounds, lundwalls, hwh, cwh, num_houses, num_classrooms, get_infection_distance, get_infection_rate, get_infection_active_rate, get_show_infection_sign_rate, get_cure_rate
end


text:
function compute_rho_cic!( rho, mesh, x, w)

dx = mesh.dx
dy = mesh.dy
dz = mesh.dz

fill!(rho, 0.0)

vol = w / (dx*dy*dz)

nbpart = size(x)[2]

for m in 1:nbpart

   xp = x[1, m]/dx
   yp = x[2, m]/dy
   zp = x[3, m]/dz

   ip = floor(Int, xp)
   jp = floor(Int, yp)
   kP = floor(Int, zp)

   dxp = xp - ip
   dyp = yp - jp
   dzp = zp - kp

   a1 = (1 - dxp) * (1 - dyp) * (1 - dzp)
   a2 = dxp * (1 - dyp) * (1 - dzp)
   a3 = (1 - dxp) * dyp * (1 - dzp)
   a4 = dxp * dyp * (1 - dzp)
   a5 = (1 - dxp) * (1 - dyp) * dzp
   a6 = dxp * (1 - dyp) * dzp
   a7 = (1 - dxp) * dyp * dzp
   a8 = dxp * dyp * dzp

   ip = ip+1
   jp = jp+1
   kp = kp+1

   rho[ip  ,jp  ,kp  ] += a1 * vol
   rho[ip+1,jp  ,kp  ] += a2 * vol
   rho[ip  ,jp+1,kp  ] += a3 * vol
   rho[ip+1,jp+1,kp  ] += a4 * vol
   rho[ip  ,jp  ,kp+1] += a5 * vol
   rho[ip+1,jp  ,kp+1] += a6 * vol
   rho[ip  ,jp+1,kp+1] += a7 * vol
   rho[ip+1,jp+1,kp+1] += a8 * vol

end

nx = mesh.nx
ny = mesh.ny
nz = mesh.nz

rho[nx+1,:,:] = rho[1,:,:]
rho[:,ny+1,:] = rho[:,1,:]
rho[:,:,nz+1] = rho[:,:,1]

end


text:
using ForwardDiff
using DelimitedFiles
using MAT

nntype = length(ARGS)>=1 ? ARGS[1] : "piecewise"
idx = length(ARGS)>=2 ? parse(Int,ARGS[2]) : 3
nn_init_id = length(ARGS)>=3 ? parse(Int,ARGS[3]) : 0

config = [20,20,20,1]

if idx == 1
    config = [20,1]
elseif idx == 2
    config = [20,20,1]
elseif idx == 3
    config = [20,20,20,1]
elseif idx == 4
    config = [20,20,20,20,1]
elseif idx == 5
    config = [20,20,20,20,20,1]
else
    error("idx <= 5")
end

stress_scale = 1e2
strain_scale = 1e-3

E0 = 200.0e3/(stress_scale/strain_scale)

function nn(ε, ε0, σ0)
    # make sure all inputs are 2d matrix
    @show ε, ε0, σ0

    local y
    if nntype=="linear"
        y = E0 * ε
    elseif nntype=="ae_scaled"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        y = ae(x, config, nntype)*stress_scale
    elseif nntype=="ae_scaled2"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        y = ae(x, config, nntype)*stress_scale + σ0
    elseif nntype=="piecewise"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        H = ae(x, config, nntype)^2
        s = σ0^2
        i = sigmoid((s - 0.01e6)/1.0e3)

        y = ( H.* i + E0 * (1-i) ) .* (ε-ε0)*stress_scale/strain_scale + σ0

    elseif nntype=="piecewise2"
        x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
        H = ae(x, config, nntype)
        s = σ0^2
        i = sigmoid((s - 0.01e6))

        y = ( H.* i + E0 * (1-i) ) .* (ε-ε0)*stress_scale/strain_scale + σ0
    else
        error("nntype must be specified.")
    end

    return y

end

function sigmoid_(z)
    return 1.0 / (1.0 + exp(-z))
end

function post_nn(ε::Float64, ε0::Float64, σ0::Float64, Δt::Float64)
    # @show "Post NN"
    local f, df
    if nntype=="linear"
        f = x -> E0*x
        df = ForwardDiff.derivative(f, ε)
    elseif nntype=="piecewise"
        f = x -> begin
            H = nnpiecewise(reshape([x/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]^2

            @show H
            s = σ0^2
            i = sigmoid_((s - 0.01e6)/1.0e3)

            (H * i + E0 * (1.0-i) )*stress_scale/strain_scale * (x-ε0) + σ0

        end
        df = ForwardDiff.derivative(f, ε)
    elseif nntype=="piecewise2"
        f = x -> begin
            H = nnpiecewise2(reshape([x/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]

            @show H

            s = σ0^2
            i = sigmoid_((s - 0.01e6))
            (H * i + E0 * (1.0-i))*stress_scale/strain_scale * (x-ε0) + σ0
        end
        df = ForwardDiff.derivative(f, ε)
    elseif nntype=="ae_scaled"
        f = x -> nnae_scaled(reshape([x/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]*stress_scale
        df = ForwardDiff.derivative(f, ε)
    elseif nntype=="ae_scaled2"
        f = x -> nnae_scaled2(reshape([x/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]*stress_scale + σ0
        df = ForwardDiff.derivative(f, ε)
    else
        error("$(nntype) must be specified.")
    end
    return f(ε), df
end

# function nn_helper(ε, ε0, σ0)
#     if nntype=="linear"
#         E0 * ε
#     elseif nntype=="piecewise"
#         H = nnpiecewise(reshape([ε/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]^2*stress_scale
#         s = σ0^2
#         i = sigmoid_((s - 0.01e6))
#         ( H * i + E0 * (1.0-i) ) * (ε-ε0) + σ0
#     elseif nntype=="ae_scaled"
#         nnae_scaled(reshape([ε/strain_scale;ε0/strain_scale;σ0/stress_scale],1,3))[1,1]*stress_scale
#     else
#         error("$nntype does not exist")
#     end
# end

# function post_nn(ε, ε0, σ0, Δt)
#     f = x -> nn_helper(x, ε0, σ0)
#     df = ForwardDiff.jacobian(f, ε)
#     return f(ε), df
# end


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using Compat
using Compat.Test
using Compat: view, String
using LegacyStrings
using LegacyStrings: ASCIIString, UTF8String # override Compat's version
import LegacyStrings:
    checkstring,
    UnicodeError,
    UTF_ERR_SHORT

# types
nullstring16 = UInt16[]
badstring16  = UInt16[0x0065]
@test_throws UnicodeError UTF16String(nullstring16)
@test_throws UnicodeError UTF16String(badstring16)

nullstring32 = UInt32[]
badstring32  = UInt32['a']
@test_throws UnicodeError UTF32String(nullstring32)
@test_throws UnicodeError UTF32String(badstring32)

# Unicode errors
let io = IOBuffer()
    show(io, UnicodeError(UTF_ERR_SHORT, 1, 10))
    check = "UnicodeError: invalid UTF-8 sequence starting at index 1 (0xa missing one or more continuation bytes)"
    @test String(take!(io)) == check
end

## Test invalid sequences

# Continuation byte not after lead
for byt in 0x80:0xbf
    @test_throws UnicodeError checkstring(UInt8[byt])
end

# Test lead bytes
for byt in 0xc0:0xff
    # Single lead byte at end of string
    @test_throws UnicodeError checkstring(UInt8[byt])
    # Lead followed by non-continuation character < 0x80
    @test_throws UnicodeError checkstring(UInt8[byt,0])
    # Lead followed by non-continuation character > 0xbf
    @test_throws UnicodeError checkstring(UInt8[byt,0xc0])
end

# Test overlong 2-byte
for byt in 0x81:0xbf
    @test_throws UnicodeError checkstring(UInt8[0xc0,byt])
end
for byt in 0x80:0xbf
    @test_throws UnicodeError checkstring(UInt8[0xc1,byt])
end

# Test overlong 3-byte
for byt in 0x80:0x9f
    @test_throws UnicodeError checkstring(UInt8[0xe0,byt,0x80])
end

# Test overlong 4-byte
for byt in 0x80:0x8f
    @test_throws UnicodeError checkstring(UInt8[0xef,byt,0x80,0x80])
end

# Test 4-byte > 0x10ffff
for byt in 0x90:0xbf
    @test_throws UnicodeError checkstring(UInt8[0xf4,byt,0x80,0x80])
end
for byt in 0xf5:0xf7
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80,0x80])
end

# Test 5-byte
for byt in 0xf8:0xfb
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80,0x80,0x80])
end

# Test 6-byte
for byt in 0xfc:0xfd
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80,0x80,0x80,0x80])
end

# Test 7-byte
@test_throws UnicodeError checkstring(UInt8[0xfe,0x80,0x80,0x80,0x80,0x80,0x80])

# Three and above byte sequences
for byt in 0xe0:0xef
    # Lead followed by only 1 continuation byte
    @test_throws UnicodeError checkstring(UInt8[byt,0x80])
    # Lead ended by non-continuation character < 0x80
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0])
    # Lead ended by non-continuation character > 0xbf
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0xc0])
end

# 3-byte encoded surrogate character(s)
# Single surrogate
@test_throws UnicodeError checkstring(UInt8[0xed,0xa0,0x80])
# Not followed by surrogate
@test_throws UnicodeError checkstring(UInt8[0xed,0xa0,0x80,0xed,0x80,0x80])
# Trailing surrogate first
@test_throws UnicodeError checkstring(UInt8[0xed,0xb0,0x80,0xed,0xb0,0x80])
# Followed by lead surrogate
@test_throws UnicodeError checkstring(UInt8[0xed,0xa0,0x80,0xed,0xa0,0x80])

# Four byte sequences
for byt in 0xf0:0xf4
    # Lead followed by only 2 continuation bytes
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80])
    # Lead followed by non-continuation character < 0x80
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80,0])
    # Lead followed by non-continuation character > 0xbf
    @test_throws UnicodeError checkstring(UInt8[byt,0x80,0x80,0xc0])
end

# Long encoding of 0x01
@test_throws UnicodeError utf8(b"\xf0\x80\x80\x80")
# Test ends of long encoded surrogates
@test_throws UnicodeError utf8(b"\xf0\x8d\xa0\x80")
@test_throws UnicodeError utf8(b"\xf0\x8d\xbf\xbf")
@test_throws UnicodeError checkstring(b"\xf0\x80\x80\x80")
@test checkstring(b"\xc0\x81"; accept_long_char=true) == (1,0x1,0,0,0)
@test checkstring(b"\xf0\x80\x80\x80"; accept_long_char=true) == (1,0x1,0,0,0)

# Surrogates
@test_throws UnicodeError checkstring(UInt16[0xd800])
@test_throws UnicodeError checkstring(UInt16[0xdc00])
@test_throws UnicodeError checkstring(UInt16[0xdc00,0xd800])

# Surrogates in UTF-32
@test_throws UnicodeError checkstring(UInt32[0xd800])
@test_throws UnicodeError checkstring(UInt32[0xdc00])
@test_throws UnicodeError checkstring(UInt32[0xdc00,0xd800])

# Characters > 0x10ffff
@test_throws UnicodeError checkstring(UInt32[0x110000])

# Test starting and different position
@test checkstring(UInt32[0x110000, 0x1f596], 2) == (1,0x10,1,0,0)

# Test valid sequences
for (seq, res) in (
    (UInt8[0x0],                (1,0,0,0,0)),   # Nul byte, beginning of ASCII range
    (UInt8[0x7f],               (1,0,0,0,0)),   # End of ASCII range
    (UInt8[0xc0,0x80],          (1,1,0,0,0)),   # Long encoded Nul byte (Modified UTF-8, Java)
    (UInt8[0xc2,0x80],          (1,2,0,0,1)),   # \u80, beginning of Latin1 range
    (UInt8[0xc3,0xbf],          (1,2,0,0,1)),   # \uff, end of Latin1 range
    (UInt8[0xc4,0x80],          (1,4,0,0,1)),   # \u100, beginning of non-Latin1 2-byte range
    (UInt8[0xdf,0xbf],          (1,4,0,0,1)),   # \u7ff, end of non-Latin1 2-byte range
    (UInt8[0xe0,0xa0,0x80],     (1,8,0,1,0)),   # \u800, beginning of 3-byte range
    (UInt8[0xed,0x9f,0xbf],     (1,8,0,1,0)),   # \ud7ff, end of first part of 3-byte range
    (UInt8[0xee,0x80,0x80],     (1,8,0,1,0)),   # \ue000, beginning of second part of 3-byte range
    (UInt8[0xef,0xbf,0xbf],     (1,8,0,1,0)),   # \uffff, end of 3-byte range
    (UInt8[0xf0,0x90,0x80,0x80],(1,16,1,0,0)),  # \U10000, beginning of 4-byte range
    (UInt8[0xf4,0x8f,0xbf,0xbf],(1,16,1,0,0)),  # \U10ffff, end of 4-byte range
    (UInt8[0xed,0xa0,0x80,0xed,0xb0,0x80], (1,0x30,1,0,0)), # Overlong \U10000, (CESU-8)
    (UInt8[0xed,0xaf,0xbf,0xed,0xbf,0xbf], (1,0x30,1,0,0)), # Overlong \U10ffff, (CESU-8)
    (UInt16[0x0000],            (1,0,0,0,0)),   # Nul byte, beginning of ASCII range
    (UInt16[0x007f],            (1,0,0,0,0)),   # End of ASCII range
    (UInt16[0x0080],            (1,2,0,0,1)),   # Beginning of Latin1 range
    (UInt16[0x00ff],            (1,2,0,0,1)),   # End of Latin1 range
    (UInt16[0x0100],            (1,4,0,0,1)),   # Beginning of non-Latin1 2-byte range
    (UInt16[0x07ff],            (1,4,0,0,1)),   # End of non-Latin1 2-byte range
    (UInt16[0x0800],            (1,8,0,1,0)),   # Beginning of 3-byte range
    (UInt16[0xd7ff],            (1,8,0,1,0)),   # End of first part of 3-byte range
    (UInt16[0xe000],            (1,8,0,1,0)),   # Beginning of second part of 3-byte range
    (UInt16[0xffff],            (1,8,0,1,0)),   # End of 3-byte range
    (UInt16[0xd800,0xdc00],     (1,16,1,0,0)),  # \U10000, beginning of 4-byte range
    (UInt16[0xdbff,0xdfff],     (1,16,1,0,0)),  # \U10ffff, end of 4-byte range
    (UInt32[0x0000],            (1,0,0,0,0)),   # Nul byte, beginning of ASCII range
    (UInt32[0x007f],            (1,0,0,0,0)),   # End of ASCII range
    (UInt32[0x0080],            (1,2,0,0,1)),   # Beginning of Latin1 range
    (UInt32[0x00ff],            (1,2,0,0,1)),   # End of Latin1 range
    (UInt32[0x0100],            (1,4,0,0,1)),   # Beginning of non-Latin1 2-byte range
    (UInt32[0x07ff],            (1,4,0,0,1)),   # End of non-Latin1 2-byte range
    (UInt32[0x0800],            (1,8,0,1,0)),   # Beginning of 3-byte range
    (UInt32[0xd7ff],            (1,8,0,1,0)),   # End of first part of 3-byte range
    (UInt32[0xe000],            (1,8,0,1,0)),   # Beginning of second part of 3-byte range
    (UInt32[0xffff],            (1,8,0,1,0)),   # End of 3-byte range
    (UInt32[0x10000],           (1,16,1,0,0)),  # \U10000, beginning of 4-byte range
    (UInt32[0x10ffff],          (1,16,1,0,0)),  # \U10ffff, end of 4-byte range
    (UInt32[0xd800,0xdc00],     (1,0x30,1,0,0)),# Overlong \U10000, (CESU-8)
    (UInt32[0xdbff,0xdfff],     (1,0x30,1,0,0)))# Overlong \U10ffff, (CESU-8)
    @test checkstring(seq) == res
end

# Test bounds checking
@test_throws BoundsError checkstring(b"abcdef", -10)
@test_throws BoundsError checkstring(b"abcdef", 0)
@test_throws BoundsError checkstring(b"abcdef", 7)
@test_throws BoundsError checkstring(b"abcdef", 3, -10)
@test_throws BoundsError checkstring(b"abcdef", 3, 0)
@test_throws BoundsError checkstring(b"abcdef", 3, 7)
@test_throws ArgumentError checkstring(b"abcdef", 3, 1)

## UTF-8 tests

# Test for CESU-8 sequences
let ch = 0x10000
    for hichar = 0xd800:0xdbff
        for lochar = 0xdc00:0xdfff
            @test convert(UTF8String, utf8(Char[hichar, lochar]).data) == string(Char(ch))
            ch += 1
        end
    end
end

let str = UTF8String(b"this is a test\xed\x80")
    @static if VERSION < v"0.7-"
        @test next(str, 15) == ('\ufffd', 16)
    else
        @test iterate(str, 15) == ('\ufffd', 16)
    end
    @test_throws BoundsError getindex(str, 0:3)
    @test_throws BoundsError getindex(str, 17:18)
    @test_throws BoundsError getindex(str, 2:17)
    @test_throws UnicodeError getindex(str, 16:17)
    # @test string(Char(0x110000)) == "\ufffd"
    sa = SubString{ASCIIString}(LegacyStrings.ascii("This is a silly test"), 1, 14)
    s8 = convert(SubString{UTF8String}, sa)
    @test typeof(s8) == SubString{UTF8String}
    @test s8 == "This is a sill"
    @test convert(UTF8String, b"this is a test\xed\x80\x80") == "this is a test\ud000"
end

# Reverse of UTF8String
@test reverse(UTF8String("")) == ""
@test reverse(UTF8String("a")) == "a"
@test reverse(UTF8String("abc")) == "cba"
@test reverse(UTF8String("xyz\uff\u800\uffff\U10ffff")) == "\U10ffff\uffff\u800\uffzyx"
for str in (b"xyz\xc1", b"xyz\xd0", b"xyz\xe0", b"xyz\xed\x80", b"xyz\xf0", b"xyz\xf0\x80",  b"xyz\xf0\x80\x80")
    @test_throws UnicodeError reverse(UTF8String(str))
end

# Specifically check UTF-8 string whose lead byte is same as a surrogate
@test convert(UTF8String,b"\xed\x9f\xbf") == "\ud7ff"

# issue #8
@test !isempty(methods(string, Tuple{Char}))

## UTF-16 tests

let u8 = "\U10ffff\U1d565\U1d7f6\U00066\U2008a"
    u16 = utf16(u8)
    @test sizeof(u16) == 18
    @test length(u16.data) == 10 && u16.data[end] == 0
    @test length(u16) == 5
    @test utf8(u16) == u8
    @test collect(u8) == collect(u16)
    @test u8 == utf16(u16.data[1:end-1]) == utf16(copyto!(Vector{UInt8}(undef, 18), 1, reinterpret(UInt8, u16.data), 1, 18))
    @test u8 == utf16(pointer(u16)) == utf16(convert(Ptr{Int16}, pointer(u16)))
    @test_throws UnicodeError utf16(utf32(Char(0x120000)))
    @test_throws UnicodeError utf16(UInt8[1,2,3])

    @test convert(UTF16String, "test") == "test"
    @test convert(UTF16String, u16) == u16
    @test convert(UTF16String, UInt16[[0x65, 0x66] [0x67, 0x68]]) == "efgh"
    @test convert(UTF16String, Int16[[0x65, 0x66] [0x67, 0x68]]) == "efgh"
    @test map(lowercase, utf16("TEST\U1f596")) == "test\U1f596"
    @test typeof(Base.unsafe_convert(Ptr{UInt16}, utf16("test"))) == Ptr{UInt16}
end

## UTF-32 tests

let u8 = "\U10ffff\U1d565\U1d7f6\U00066\U2008a"
    u32 = utf32(u8)
    @test sizeof(u32) == 20
    @test length(u32.data) == 6 && u32.data[end] == 0
    @test length(u32) == 5
    @test utf8(u32) == u8
    @test collect(u8) == collect(u32)
    @test u8 == utf32(u32.data[1:end-1]) == utf32(copyto!(Vector{UInt8}(undef, 20), 1, reinterpret(UInt8, u32.data), 1, 20))
    @test u8 == utf32(pointer(u32)) == utf32(convert(Ptr{Int32}, pointer(u32)))
    @test_throws UnicodeError utf32(UInt8[1,2,3])
end

# issue #11551 (#11004,#10959)
function tstcvt(strUTF8::UTF8String, strUTF16::UTF16String, strUTF32::UTF32String)
    @test utf16(strUTF8) == strUTF16
    @test utf32(strUTF8) == strUTF32
    @test utf8(strUTF16) == strUTF8
    @test utf32(strUTF16) == strUTF32
    @test utf8(strUTF32)  == strUTF8
    @test utf16(strUTF32) == strUTF16
end

# Create some ASCII, UTF8, UTF16, and UTF32 strings
strAscii = LegacyStrings.ascii("abcdefgh")
strA_UTF8 = utf8(("abcdefgh\uff")[1:8])
strL_UTF8 = utf8("abcdef\uff\uff")
str2_UTF8 = utf8("abcd\uff\uff\u7ff\u7ff")
str3_UTF8 = utf8("abcd\uff\uff\u7fff\u7fff")
str4_UTF8 = utf8("abcd\uff\u7ff\u7fff\U7ffff")
strS_UTF8 = UTF8String(b"abcd\xc3\xbf\xdf\xbf\xe7\xbf\xbf\xed\xa0\x80\xed\xb0\x80")
strC_UTF8 = UTF8String(b"abcd\xc3\xbf\xdf\xbf\xe7\xbf\xbf\U10000")
strz_UTF8 = UTF8String(b"abcd\xc3\xbf\xdf\xbf\xe7\xbf\xbf\0")
strZ      = b"abcd\xc3\xbf\xdf\xbf\xe7\xbf\xbf\xc0\x80"

strA_UTF16 = utf16(strA_UTF8)
strL_UTF16 = utf16(strL_UTF8)
str2_UTF16 = utf16(str2_UTF8)
str3_UTF16 = utf16(str3_UTF8)
str4_UTF16 = utf16(str4_UTF8)
strS_UTF16 = utf16(strS_UTF8)

strA_UTF32 = utf32(strA_UTF8)
strL_UTF32 = utf32(strL_UTF8)
str2_UTF32 = utf32(str2_UTF8)
str3_UTF32 = utf32(str3_UTF8)
str4_UTF32 = utf32(str4_UTF8)
strS_UTF32 = utf32(strS_UTF8)

@test utf8(strAscii) == strAscii
@test utf16(strAscii) == strAscii
@test utf32(strAscii) == strAscii

tstcvt(strA_UTF8,strA_UTF16,strA_UTF32)
tstcvt(strL_UTF8,strL_UTF16,strL_UTF32)
tstcvt(str2_UTF8,str2_UTF16,str2_UTF32)
tstcvt(str3_UTF8,str3_UTF16,str3_UTF32)
tstcvt(str4_UTF8,str4_UTF16,str4_UTF32)

# Test converting surrogate pairs
@test utf16(strS_UTF8) == strC_UTF8
@test utf32(strS_UTF8) == strC_UTF8
@test utf8(strS_UTF16) == strC_UTF8
@test utf32(strS_UTF16) == strC_UTF8
@test utf8(strS_UTF32)  == strC_UTF8
@test utf16(strS_UTF32) == strC_UTF8

# Test converting overlong \0
@test utf8(strZ)  == strz_UTF8
@test utf16(UTF8String(strZ)) == strz_UTF8
@test utf32(UTF8String(strZ)) == strz_UTF8

# Test invalid sequences

strval(::Type{UTF8String}, dat) = dat
strval(::Union{Type{UTF16String},Type{UTF32String}}, dat) = UTF8String(dat)

for T in (UTF8String, UTF16String, UTF32String)
    # Continuation byte not after lead
    for byt in 0x80:0xbf
        @test_throws UnicodeError convert(T,  strval(T, UInt8[byt]))
    end

    # Test lead bytes
    for byt in 0xc0:0xff
        # Single lead byte at end of string
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt]))
        # Lead followed by non-continuation character < 0x80
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0]))
        # Lead followed by non-continuation character > 0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0xc0]))
    end

    # Test overlong 2-byte
    for byt in 0x81:0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[0xc0,byt]))
    end
    for byt in 0x80:0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[0xc1,byt]))
    end

    # Test overlong 3-byte
    for byt in 0x80:0x9f
        @test_throws UnicodeError convert(T, strval(T, UInt8[0xe0,byt,0x80]))
    end

    # Test overlong 4-byte
    for byt in 0x80:0x8f
        @test_throws UnicodeError convert(T, strval(T, UInt8[0xef,byt,0x80,0x80]))
    end

    # Test 4-byte > 0x10ffff
    for byt in 0x90:0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[0xf4,byt,0x80,0x80]))
    end
    for byt in 0xf5:0xf7
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80,0x80]))
    end

    # Test 5-byte
    for byt in 0xf8:0xfb
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80,0x80,0x80]))
    end

    # Test 6-byte
    for byt in 0xfc:0xfd
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80,0x80,0x80,0x80]))
    end

    # Test 7-byte
    @test_throws UnicodeError convert(T, strval(T, UInt8[0xfe,0x80,0x80,0x80,0x80,0x80,0x80]))

    # Three and above byte sequences
    for byt in 0xe0:0xef
        # Lead followed by only 1 continuation byte
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80]))
        # Lead ended by non-continuation character < 0x80
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0]))
        # Lead ended by non-continuation character > 0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0xc0]))
    end

    # 3-byte encoded surrogate character(s)
    # Single surrogate
    @test_throws UnicodeError convert(T, strval(T, UInt8[0xed,0xa0,0x80]))
    # Not followed by surrogate
    @test_throws UnicodeError convert(T, strval(T, UInt8[0xed,0xa0,0x80,0xed,0x80,0x80]))
    # Trailing surrogate first
    @test_throws UnicodeError convert(T, strval(T, UInt8[0xed,0xb0,0x80,0xed,0xb0,0x80]))
    # Followed by lead surrogate
    @test_throws UnicodeError convert(T, strval(T, UInt8[0xed,0xa0,0x80,0xed,0xa0,0x80]))

    # Four byte sequences
    for byt in 0xf0:0xf4
        # Lead followed by only 2 continuation bytes
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80]))
        # Lead followed by non-continuation character < 0x80
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80,0]))
        # Lead followed by non-continuation character > 0xbf
        @test_throws UnicodeError convert(T, strval(T, UInt8[byt,0x80,0x80,0xc0]))
    end
end

# Wstring
let u8 = "\U10ffff\U1d565\U1d7f6\U00066\U2008a"
    w = wstring(u8)
    @test length(w) == 5 && utf8(w) == u8 && collect(u8) == collect(w)
    @test u8 == WString(w.data)
end

# 12268
for (fun, S, T) in ((utf16, UInt16, UTF16String), (utf32, UInt32, UTF32String))
    # AbstractString
    str = "abcd\0\uff\u7ff\u7fff\U7ffff"
    tst = SubString(convert(T,str),4)
    cmp = Char['d','\0','\uff','\u7ff','\u7fff','\U7ffff']
    cmp32 = UInt32['d','\0','\uff','\u7ff','\u7fff','\U7ffff','\0']
    cmp16 = UInt16[0x0064,0x0000,0x00ff,0x07ff,0x7fff,0xd9bf,0xdfff,0x0000]
    x = fun(tst)
    cmpx = (S == UInt16 ? cmp16 : cmp32)
    @test typeof(tst) == SubString{T}
    @test convert(T, tst) == str[4:end]
    @test Vector{Char}(x) == cmp
    # Vector{T} / Array{T}
    @test convert(Vector{S}, x) == cmpx
    @test convert(Array{S}, x) == cmpx
    # Embedded nul checking
    @test Base.containsnul(x)
    @test Base.containsnul(tst)
    # map
    @test_throws UnicodeError map(islowercase, x)
    @test_throws ArgumentError map(islowercase, tst)
    # SubArray conversion
    subarr = view(cmp, 1:6)
    @test convert(T, subarr) == str[4:end]
end

# Char to UTF32String
@test utf32('\U7ffff') == utf32("\U7ffff")
@test convert(UTF32String, '\U7ffff') == utf32("\U7ffff")

@test isvalid(UTF32String, Char['d','\uff','\u7ff','\u7fff','\U7ffff'])
@test reverse(utf32("abcd \uff\u7ff\u7fff\U7ffff")) == utf32("\U7ffff\u7fff\u7ff\uff dcba")

# Test pointer() functions
let str = LegacyStrings.ascii("this ")
    u8  = utf8(str)
    u16 = utf16(str)
    u32 = utf32(str)
    pa  = pointer(str)
    p8  = pointer(u8)
    p16 = pointer(u16)
    p32 = pointer(u32)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x74
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x74
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x74
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x74
    pa  = pointer(str, 2)
    p8  = pointer(u8,  2)
    p16 = pointer(u16, 2)
    p32 = pointer(u32, 2)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x68
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x68
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x68
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x68
    sa  = SubString{ASCIIString}(str, 3, 5)
    s8  = SubString{UTF8String}(u8,   3, 5)
    s16 = SubString{UTF16String}(u16, 3, 5)
    s32 = SubString{UTF32String}(u32, 3, 5)
    pa  = pointer(sa)
    p8  = pointer(s8)
    p16 = pointer(s16)
    p32 = pointer(s32)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x69
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x69
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x69
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x69
    pa  = pointer(sa, 2)
    p8  = pointer(s8,  2)
    p16 = pointer(s16, 2)
    p32 = pointer(s32, 2)
    @test typeof(pa) == Ptr{UInt8}
    @test unsafe_load(pa,1) == 0x73
    @test typeof(p8) == Ptr{UInt8}
    @test unsafe_load(p8,1) == 0x73
    @test typeof(p16) == Ptr{UInt16}
    @test unsafe_load(p16,1) == 0x73
    @test typeof(p32) == Ptr{UInt32}
    @test unsafe_load(p32,1) == 0x73
end

@test isvalid(Char['f','o','o','b','a','r'])

## Repeat strings ##

# Base Julia issue #7764
let
    rs = RepString("foo", 2)
    @test length(rs) == 6
    @test sizeof(rs) == 6
    @test isascii(rs)
    @test convert(RepString, "foobar") == "foobar"
    @test typeof(convert(RepString, "foobar")) == RepString

    srep = RepString("Σβ",2)
    s="Σβ"
    ss=SubString(s,1,lastindex(s))

    @test ss^2 == "ΣβΣβ"
    @test RepString(ss,2) == "ΣβΣβ"

    @test lastindex(srep) == 7

    @static if VERSION < v"0.7-"
        @test next(srep, 3) == ('β',5)
        @test next(srep, 7) == ('β',9)
    else
        @test iterate(srep, 3) == ('β',5)
        @test iterate(srep, 7) == ('β',9)
    end

    @test srep[7] == 'β'
    @static if VERSION < v"0.7.0-DEV.2924"
        @test_throws BoundsError srep[8]
    else
        @test_throws StringIndexError srep[8]
    end
end

## Reverse strings ##

let
    rs = RevString("foobar")
    @test length(rs) == 6
    @test sizeof(rs) == 6
    @test isascii(rs)

    # Base issue #4586
    @test rsplit(RevString("ailuj"),'l') == ["ju","ia"]
    @test parse(Float64,RevString("64")) === 46.0

    # reverseind
    for prefix in ("", "abcd", "\U0001d6a4\U0001d4c1", "\U0001d6a4\U0001d4c1c", " \U0001d6a4\U0001d4c1")
        for suffix in ("", "abcde", "\U0001d4c1β\U0001d6a4", "\U0001d4c1β\U0001d6a4c", " \U0001d4c1β\U0001d6a4")
            for c in ('X', 'δ', '\U0001d6a5')
                s = convert(String, string(prefix, c, suffix))
                rs = RevString(s)
                r = reverse(s)
                @test r == rs
                ri = something(findfirst(isequal(c), r), 0)
                @test c == s[reverseind(s, ri)] == r[ri]
            end
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.05"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	n1, n2 = 10, 15
	N = 10^6
	mu, sigma = 10, 4
	normDist = Normal(mu,sigma)
end;

# ╔═╡ 67a42f40-3a99-11eb-2199-cb72cdee58fa
begin
	fValues = Array{Float64}(undef, N)

	for i in 1:N
		data1 = rand(normDist,n1)
		data2 = rand(normDist,n2)
		fValues[i] = var(data1)/var(data2)
	end
end

# ╔═╡ 67a46e56-3a99-11eb-1254-3d0d490364f5
begin
	fRange = 0:0.1:5
	stephist(fValues, bins=400, c=:blue, label="Simulated", normed=true)
	plot!(fRange, pdf.(FDist(n1-1, n2-1), fRange),
		c=:red, label="Analytic", xlims=(0,5), ylims=(0,0.8),
			xlabel = "F", ylabel = "Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.05"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═67a42f40-3a99-11eb-2199-cb72cdee58fa
# ╠═67a46e56-3a99-11eb-1254-3d0d490364f5
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
@testset "Find" begin
    kmer = DNAMer("ACGAG")
    bigkmer = BigDNAMer("ACGAG")

    @test findnext(DNA_A, kmer, 1) == 1
    @test findnext(DNA_C, kmer, 1) == 2
    @test findnext(DNA_G, kmer, 1) == 3
    @test findnext(DNA_T, kmer, 1) == nothing
    @test findnext(DNA_A, kmer, 2) == 4

    @test findnext(DNA_A, bigkmer, 1) == 1
    @test findnext(DNA_C, bigkmer, 1) == 2
    @test findnext(DNA_G, bigkmer, 1) == 3
    @test findnext(DNA_T, bigkmer, 1) == nothing
    @test findnext(DNA_A, bigkmer, 2) == 4

    @test_throws BoundsError findnext(DNA_A, kmer, 0)
    @test_throws BoundsError findnext(DNA_A, kmer, 6)

    @test_throws BoundsError findnext(DNA_A, bigkmer, 0)
    @test_throws BoundsError findnext(DNA_A, bigkmer, 6)

    @test findprev(DNA_A, kmer, 5) == 4
    @test findprev(DNA_C, kmer, 5) == 2
    @test findprev(DNA_G, kmer, 5) == 5
    @test findprev(DNA_T, kmer, 5) == nothing
    @test findprev(DNA_G, kmer, 4) == 3

    @test findprev(DNA_A, bigkmer, 5) == 4
    @test findprev(DNA_C, bigkmer, 5) == 2
    @test findprev(DNA_G, bigkmer, 5) == 5
    @test findprev(DNA_T, bigkmer, 5) == nothing
    @test findprev(DNA_G, bigkmer, 4) == 3

    @test_throws BoundsError findprev(DNA_A, kmer, 0)
    @test_throws BoundsError findprev(DNA_A, kmer, 6)

    @test_throws BoundsError findprev(DNA_A, bigkmer, 0)
    @test_throws BoundsError findprev(DNA_A, bigkmer, 6)

    @test findfirst(DNA_A, kmer) == 1
    @test findfirst(DNA_G, kmer) == 3
    @test findlast(DNA_A, kmer) == 4
    @test findlast(DNA_G, kmer) == 5

    @test findfirst(DNA_A, bigkmer) == 1
    @test findfirst(DNA_G, bigkmer) == 3
    @test findlast(DNA_A, bigkmer) == 4
    @test findlast(DNA_G, bigkmer) == 5
end


text:
# Tests for validated_integ

using TaylorModels
# using LinearAlgebra: norm
using Test
using Random

const _num_tests = 1_000

setformat(:full)

# NOTE: IntervalArithmetic v0.16.0 includes this function; but
# IntervalRootFinding is bounded to use v0.15.x
interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)
interval_rand(X::IntervalBox) = interval_rand.(X)

# Function to check, against an exact solution of the ODE, the computed
# validted solution
function test_integ(fexact, t0, qTM, q0, δq0)
    normalized_box = symmetric_box(length(q0), Float64)
    # Time domain
    domt = domain(qTM[1])
    # Random time (within time domain) and random initial condition
    δt = rand(domt)
    δtI = (δt .. δt) ∩ domt
    q0ξ = interval_rand(δq0)
    q0ξB = IntervalBox([(q0ξ[i] .. q0ξ[i]) ∩ δq0[i] for i in eachindex(q0ξ)])
    # Box computed to overapproximate the solution at time δt
    q = evaluate.(evaluate.(qTM, δtI), (normalized_box,))
    # Box computed from the exact solution must be within q
    bb = all(fexact(t0+δtI, q0 .+ q0ξB) .⊆ q)
    # Display details if bb is false
    bb || @show(t0, domt, remainder.(qTM),
            δt, δtI, q0ξ, q0ξB, q,
            fexact(t0+δtI, q0 .+ q0ξB))
    return bb
end

@testset "Tests for `validated_integ`" begin
    @testset "falling_ball!" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) = (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 0.0, 10.0
        normalized_box = symmetric_box(2, Float64)
        q0 = [10.0, 0.0]
        δq0 = 0.25 * normalized_box
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))

        @testset "Forward integration 1" begin
            sol = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
            tTM = expansion_point(sol)
            qv  = flowpipe(sol)
            qTM = get_xTM(sol)
            @test length(qv) == length(qTM[1,:]) == length(sol)
            @test firstindex(sol) == 1
            @test sol[2] == get_xTM(sol,2)
            @test domain(sol,1) == 0..0

            end_idx = lastindex(sol)
            Random.seed!(1)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
            end

            solf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            qvf, qTMf = getfield.((solf,), 2:3)

            @test length(qvf) == length(qv)
            @test qTM == qTMf

            # initializaton with a Taylor model
            X0tm = sol[1]
            sol2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            qTM2 = getfield(sol2, 3)
            @test qTM == qTM2
            @test sol2[1,2] == X0tm[2]

            # Tests for TaylorModels.mince_in_time
            domT = TaylorModels.mince_in_time(sol)
            @test domT == expansion_point(sol) .+ domain(sol)
            timesdiv = TaylorModels.mince_in_time(sol, var=0, timediv=2)
            fpdiv = TaylorModels.mince_in_time(sol, var=1, timediv=2)
            @test timesdiv[3] ⊂ domT[2]
            @test hull(timesdiv[1],timesdiv[2]) == domT[1]
            @test fpdiv[3] ⊂ qv[2][1]
            @test hull(fpdiv[3],fpdiv[4]) ⊂ qv[2][1]
        end

        @testset "Forward integration 2" begin
            sol = validated_integ2(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
            tTM = expansion_point(sol)
            qv  = flowpipe(sol)
            qTM = get_xTM(sol)
            @test length(qv) == length(qTM[1,:]) == length(sol)
            @test firstindex(sol) == 1
            @test sol[2] == get_xTM(sol,2)
            @test domain(sol,1) == 0..0

            Random.seed!(1)
            end_idx = lastindex(sol)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), tTM[n], sol[n], q0, δq0)
            end

            # initializaton with a Taylor model
            X0tm = get_xTM(sol,1)
            sol2 = validated_integ2(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            qTM2 = getfield(sol2, 3)
            @test qTM == qTM2
            @test sol2[1,2] == X0tm[2]
        end

        # Initial conditions
        tini, tend = 10.0, 0.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        @testset "Backward integration 1" begin
            sol = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
            tTM, qv, qTM = getfield.((sol,), 1:3)
            @test length(qv) == length(qTM[1,:]) == length(sol)
            @test firstindex(sol) == 1
            @test sol[2] == get_xTM(sol,2)
            @test domain(sol,1) == 0..0

            Random.seed!(1)
            end_idx = lastindex(sol)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), expansion_point(sol,n), sol[n], q0, δq0)
            end

            solf = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            tTMf, qvf, qTMf = getfield.((solf,), 1:3)

            @test length(qvf) == length(qv)
            @test all(qTM .== qTMf)

            # initializaton with a Taylor model
            X0tm = sol[1]
            sol2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
            qv2, qTM2 = getfield.((sol2,), 2:3)
            @test qTM == qTM2

            sol2f = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            qv2f, qTM2f = getfield.((sol2f,), 2:3)
            @test length(qv2f) == length(qv2)
            @test all(qTM .== qTM2f)
        end

        @testset "Backward integration 2" begin
            sol = validated_integ2(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)
            tTM, qv, qTM = getfield.((sol,), 1:3)
            @test length(qv) == length(qTM[1,:]) == length(sol)
            @test firstindex(sol) == 1
            @test sol[2] == get_xTM(sol,2)
            @test domain(sol,1) == 0..0

            Random.seed!(1)
            end_idx = lastindex(sol)
            for it = 1:_num_tests
                n = rand(2:end_idx)
                @test test_integ((t,x)->exactsol(t,tini,x), expansion_point(sol,n), sol[n], q0, δq0)
            end
        end
    end

    @testset "x_square!" begin
        @taylorize function x_square!(dx, x, p, t)
            dx[1] = x[1]^2
            nothing
        end

        exactsol(t, x0) = 1 / (1/x0[1] - t)

        tini, tend = 0., 0.45
        normalized_box = symmetric_box(1, Float64)
        abstol = 1e-15
        orderQ = 5
        orderT = 20
        q0 = [2.]
        δq0 = 0.0625 * normalized_box
        X0 = IntervalBox(q0 .+ δq0)
        ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)

        @testset "Forward integration 1" begin
            sol = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol)
            tTM, qv, qTM = getfield.((sol,), 1:3)
            @test length(qv) == length(qTM[1, :]) == length(tTM)
            @test domain(sol,1) == 0..0

            Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(1:end_idx)
                @test test_integ((t,x)->exactsol(t,x), tTM[n], sol[n], q0, δq0)
            end

            solf = validated_integ(x_square!, X0, tini, tend, orderQ, orderT, abstol,
                adaptive=false)
            tTMf, qvf, qTMf = getfield.((solf,), 1:3)

            @test length(qvf) == length(qv)
            @test all(qTMf .== qTM)

            # initializaton with a Taylor model
            X0tm = copy(qTM[:, 1])
            sol2 = validated_integ(x_square!, X0tm, tini, tend, orderQ, orderT, abstol)
            tTM2, qv2, qTM2 = getfield.((sol2,), 1:3)
            @test qTM == qTM2
        end

        @testset "Forward integration 2" begin
            sol = validated_integ2(x_square!, X0, tini, tend, orderQ, orderT, abstol)
            tTM, qv, qTM = getfield.((sol,), 1:3)
            @test domain(sol,1) == 0..0

            @test length(qv) == length(qTM[1, :]) == length(tTM)

            Random.seed!(1)
            end_idx = lastindex(tTM)
            for it = 1:_num_tests
                n = rand(1:end_idx)
                @test test_integ((t,x)->exactsol(t,x), tTM[n], sol[n], q0, δq0)
            end
        end
    end

    @testset "Pendulum with constant torque" begin
        @taylorize function pendulum!(dx, x, p, t)
            si = sin(x[1])
            aux = 2 *  si
            dx[1] = x[2]
            dx[2] = aux + 8*x[3]
            dx[3] = zero(x[1])
            nothing
        end
        # Conserved quantity
        ene_pendulum(x) = x[2]^2/2 + 2 * cos(x[1]) - 8 * x[3]

        # Initial conditions
        tini, tend = 0.0, 12.0
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, 0..0)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        # Parameters
        abstol = 1e-10
        orderQ = 3
        orderT = 10
        ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))

        sol = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))

        sol = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            validatesteps=32);
        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))

        # Initial conditions 2
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, -0.01 .. 0.01)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        sol = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol);
        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))

        sol = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            validatesteps=32);
        @test all(ene0 .⊆ ene_pendulum.(flowpipe(sol)))
    end
end