text:
## Common Interface Solve Functions

function DiffEqBase.__solve(
  prob::Union{DiffEqBase.AbstractODEProblem,DiffEqBase.AbstractDAEProblem},
  alg::algType,timeseries=[],ts=[],ks=[],
  recompile::Type{Val{recompile_flag}}=Val{true};
  kwargs...) where {algType<:Union{SundialsODEAlgorithm,SundialsDAEAlgorithm},
                    recompile_flag}

  integrator = DiffEqBase.init(prob,alg,timeseries,ts,ks;kwargs...)
  if integrator.sol.retcode == :Default
    solve!(integrator)
  end
  integrator.sol
end

function DiffEqBase.__init(
    prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace},
    alg::SundialsODEAlgorithm{Method,LinearSolver},
    timeseries=[], ts=[], ks=[];

    verbose=true,
    callback=nothing, abstol=1/10^6, reltol=1/10^3,
    saveat=Float64[], tstops=Float64[],
    maxiters=Int(1e5),
    dt = nothing, dtmin = 0.0, dtmax = 0.0,
    timeseries_errors=true,
    dense_errors = false,
    save_everystep=isempty(saveat),
    save_on = true,
    save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[1] in saveat,
    save_end = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[2] in saveat,
    dense = save_everystep && isempty(saveat),
    progress=false,progress_name="ODE",
    progress_message = DiffEqBase.ODE_DEFAULT_PROG_MESSAGE,
    save_timeseries = nothing,
    advance_to_tstop = false,stop_at_next_tstop=false,
    userdata=nothing,
    alias_u0=false,
    kwargs...) where {uType, tupType, isinplace, Method, LinearSolver}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction) && typeof(alg) <: CVODE_BDF
            if DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    if prob.f.mass_matrix != I
        error("This solver is not able to use mass matrices.")
    end

    if typeof(reltol) <: AbstractArray
        error("Sundials only allows scalar reltol.")
    end

    progress && @logmsg(-1,progress_name,_id=_id = :Sundials,progress=0)

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
    if max_len_cb isa VectorContinuousCallback
      callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltype,uBottomEltype)
    else
      callback_cache = nothing
    end

    tspan = prob.tspan
    t0 = tspan[1]

    tdir = sign(tspan[2]-tspan[1])

    tstops_internal, saveat_internal =
      tstop_saveat_disc_handling(tstops,saveat,tdir,tspan,tType)

    if typeof(prob.u0) <: Number
        u0 = [prob.u0]
    else
        if alias_u0
            u0 = vec(prob.u0)
        else
            u0 = vec(copy(prob.u0))
        end
    end

    sizeu = size(prob.u0)

    ### Fix the more general function to Sundials allowed style
    if !isinplace && typeof(prob.u0)<:Number
        f! = (du, u, p, t) -> (du .= prob.f(first(u), p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:Vector{Float64}
        f! = (du, u, p, t) -> (du .= prob.f(u, p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:AbstractArray
        f! = (du, u, p, t) -> (du .= vec(prob.f(reshape(u, sizeu), p, t)); Cint(0))
    elseif typeof(prob.u0)<:Vector{Float64}
        f! = prob.f
    else # Then it's an in-place function on an abstract array
        f! = (du, u, p, t) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p, t);
                               du=vec(du); 0)
    end

    if typeof(alg) <: CVODE_BDF
        alg_code = CV_BDF
    elseif typeof(alg) <:  CVODE_Adams
        alg_code = CV_ADAMS
    end

    if Method == :Newton
        method_code = CV_NEWTON
    elseif Method ==  :Functional
        method_code = CV_FUNCTIONAL
    end

    mem_ptr = CVodeCreate(alg_code, method_code)
    (mem_ptr == C_NULL) && error("Failed to allocate CVODE solver object")
    mem = Handle(mem_ptr)

    !verbose && CVodeSetErrHandlerFn(mem,@cfunction(null_error_handler, Nothing,
                                    (Cint, Char,
                                    Char, Ptr{Cvoid})),C_NULL)

    ures  = Vector{uType}()
    dures = Vector{uType}()
    save_start ? ts = [t0] : ts = Float64[]

    u0nv = NVector(u0)
    _u0 = copy(u0)
    utmp = NVector(_u0)

    userfun = FunJac(f!,prob.f.jac,prob.p,nothing,prob.f.jac_prototype,alg.prec,alg.psetup,u0,_u0)

    function getcfunf(userfun::T) where T
        @cfunction(cvodefunjac, Cint, (realtype, N_Vector, N_Vector, Ref{T}))
    end

    flag = CVodeInit(mem,getcfunf(userfun),
                     t0, convert(N_Vector, utmp))

    dt != nothing && (flag = CVodeSetInitStep(mem, dt))
    flag = CVodeSetMinStep(mem, dtmin)
    flag = CVodeSetMaxStep(mem, dtmax)
    flag = CVodeSetUserData(mem, userfun)
    if typeof(abstol) <: Array
        flag = CVodeSVtolerances(mem, reltol, abstol)
    else
        flag = CVodeSStolerances(mem, reltol, abstol)
    end
    flag = CVodeSetMaxNumSteps(mem, maxiters)
    flag = CVodeSetMaxOrd(mem, alg.max_order)
    flag = CVodeSetMaxHnilWarns(mem, alg.max_hnil_warns)
    flag = CVodeSetStabLimDet(mem, alg.stability_limit_detect)
    flag = CVodeSetMaxErrTestFails(mem, alg.max_error_test_failures)
    flag = CVodeSetMaxNonlinIters(mem, alg.max_nonlinear_iters)
    flag = CVodeSetMaxConvFails(mem, alg.max_convergence_failures)

    if Method == :Newton # Only use a linear solver if it's a Newton-based method
        if LinearSolver == :Dense
            A = SUNDenseMatrix(length(u0),length(u0))
            LS = SUNDenseLinearSolver(u0,A)
            flag = CVDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,DenseMatrix())
            _LS = LinSolHandle(LS,Dense())
        elseif LinearSolver == :Band
            A = SUNBandMatrix(length(u0), alg.jac_upper, alg.jac_lower,
                                       alg.stored_upper)
            LS = SUNBandLinearSolver(u0,A)
            flag = CVDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,BandMatrix())
            _LS = LinSolHandle(LS,Band())
        elseif LinearSolver == :Diagonal
            flag = CVDiag(mem)
            _A = nothing
            _LS = nothing
        elseif LinearSolver == :GMRES
            LS = SUNSPGMR(u0, alg.prec_side, alg.krylov_dim)
            flag = CVSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = Sundials.LinSolHandle(LS,Sundials.SPGMR())
        elseif LinearSolver == :FGMRES
            LS = SUNSPFGMR(u0, alg.prec_side, alg.krylov_dim)
            flag = CVSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,SPFGMR())
        elseif LinearSolver == :BCG
            LS = SUNSPBCGS(u0, alg.prec_side, alg.krylov_dim)
            flag = CVSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,SPBCGS())
        elseif LinearSolver == :PCG
            LS = SUNPCG(u0, alg.prec_side, alg.krylov_dim)
            flag = CVSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,PCG())
        elseif LinearSolver == :TFQMR
            LS = SUNSPTFQMR(u0, alg.prec_side, alg.krylov_dim)
            flag = CVSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,PTFQMR())
	    elseif LinearSolver == :KLU
	        nnz = length(nonzeros(prob.f.jac_prototype))
	        A = SUNSparseMatrix(length(u0),length(u0), nnz, CSC_MAT)
            LS = SUNKLU(u0, A)
            flag = CVDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,SparseMatrix())
            _LS = LinSolHandle(LS,KLU())
        end
    else
        _A = nothing
        _LS = nothing
    end

    if DiffEqBase.has_jac(prob.f) && Method == :Newton
      function getcfunjac(userfun::T) where T
          @cfunction(cvodejac,
                          Cint,
                           (realtype,
                           N_Vector,
                           N_Vector,
                           SUNMatrix,
                           Ref{T},
                           N_Vector,
                           N_Vector,
                           N_Vector))
      end
      jac = getcfunjac(userfun)
      flag = CVodeSetUserData(mem, userfun)
      flag = CVDlsSetJacFn(mem, jac)
    else
        jac = nothing
    end

    if typeof(prob.f.jac_prototype) <: DiffEqBase.AbstractDiffEqLinearOperator
        function getcfunjtimes(userfun::T) where T
            @cfunction(jactimes,
                            Cint,
                            (N_Vector,
                             N_Vector,
                             realtype,
                             N_Vector,
                             N_Vector,
                             Ref{T},
                             N_Vector))
        end
        jtimes = getcfunjtimes(userfun)
        CVSpilsSetJacTimes(mem, C_NULL, jtimes)
    end

    if alg.prec !== nothing
        function getpercfun(userfun::T) where T
            @cfunction(precsolve,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             N_Vector,
                             N_Vector,Float64,Float64,Int,
                             Ref{T}))
        end
        precfun = getpercfun(userfun)

        function getpsetupfun(userfun::T) where T
            @cfunction(precsetup,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             Int,
                             Ptr{Int},Float64,Ref{T}))
        end
        psetupfun = alg.psetup === nothing ? C_NULL : getpsetupfun(userfun)

        CVSpilsSetPreconditioner(mem, psetupfun, precfun)
    end

    callbacks_internal == nothing ? tmp = nothing : tmp = similar(u0)
    callbacks_internal == nothing ? uprev = nothing : uprev = similar(u0)
    tout = [tspan[1]]

    if save_start
      save_value!(ures,u0,uType,sizeu)
      if dense
        f!(_u0,u0,prob.p,tspan[1])
        save_value!(dures,utmp,uType,sizeu)
      end
    end

    sol = DiffEqBase.build_solution(prob, alg, ts, ures,
                   dense = dense,
                   interp = dense ? DiffEqBase.HermiteInterpolation(ts,ures,dures) :
                                    DiffEqBase.LinearInterpolation(ts,ures),
                   timeseries_errors = timeseries_errors,
                   destats = DiffEqBase.DEStats(0),
                   calculate_error = false)
    opts = DEOptions(saveat_internal,tstops_internal,save_everystep,dense,
                     timeseries_errors,dense_errors,save_on,save_end,
                     callbacks_internal,abstol,reltol,verbose,advance_to_tstop,stop_at_next_tstop,
                     progress,progress_name,progress_message)
    integrator = CVODEIntegrator(u0,prob.p,t0,t0,mem,_LS,_A,sol,alg,f!,userfun,jac,opts,
                       tout,tdir,sizeu,false,tmp,uprev,Cint(flag),false,0,1,callback_cache,0.)

    initialize_callbacks!(integrator)
    integrator
end # function solve

function DiffEqBase.__init(
    prob::DiffEqBase.AbstractODEProblem{uType, tupType, isinplace},
    alg::ARKODE{Method,LinearSolver,MassLinearSolver},
    timeseries=[], ts=[], ks=[];

    verbose=true,
    callback=nothing, abstol=1/10^6, reltol=1/10^3,
    saveat=Float64[], tstops=Float64[],
    maxiters=Int(1e5),
    dt = nothing, dtmin = 0.0, dtmax = 0.0,
    timeseries_errors=true,
    dense_errors = false,
    save_everystep=isempty(saveat), dense = save_everystep,
    save_on = true, save_start = true, save_end = true,
    save_timeseries = nothing,
    progress=false,progress_name="ODE",
    progress_message = DiffEqBase.ODE_DEFAULT_PROG_MESSAGE,
    advance_to_tstop = false,stop_at_next_tstop=false,
    userdata=nothing,
    alias_u0=false,
    kwargs...) where {uType, tupType, isinplace, Method, LinearSolver, MassLinearSolver}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction)
            if typeof(prob.f) <: SplitFunction ? DiffEqBase.has_tgrad(prob.f.f1) : DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    if typeof(reltol) <: AbstractArray
        error("Sundials only allows scalar reltol.")
    end

    progress && @logmsg(-1,progress_name,_id=_id = :Sundials,progress=0)

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
    if max_len_cb isa VectorContinuousCallback
      callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltype,uBottomEltype)
    else
      callback_cache = nothing
    end

    tspan = prob.tspan
    t0 = tspan[1]

    tdir = sign(tspan[2]-tspan[1])

    tstops_internal, saveat_internal =
      tstop_saveat_disc_handling(tstops,saveat,tdir,tspan,tType)

    if typeof(prob.u0) <: Number
        u0 = [prob.u0]
    else
        if alias_u0
            u0 = vec(prob.u0)
        else
            u0 = vec(copy(prob.u0))
        end
    end

    sizeu = size(prob.u0)

    mem_ptr = ARKodeCreate()
    (mem_ptr == C_NULL) && error("Failed to allocate ARKODE solver object")
    mem = Handle(mem_ptr)

    !verbose && ARKodeSetErrHandlerFn(mem,@cfunction(null_error_handler, Nothing,
                                    (Cint, Char,
                                    Char, Ptr{Cvoid})),C_NULL)

    ures  = Vector{uType}()
    dures = Vector{uType}()
    save_start ? ts = [t0] : ts = Float64[]
    u0nv = NVector(u0)
    _u0 = copy(u0)
    utmp = NVector(_u0)

    ### Fix the more general function to Sundials allowed style
    if !isinplace && typeof(prob.u0)<:Number
        f! = (du, u, p, t) -> (du .= prob.f(first(u), p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:Vector{Float64}
        f! = (du, u, p, t) -> (du .= prob.f(u, p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:AbstractArray
        f! = (du, u, p, t) -> (du .= vec(prob.f(reshape(u, sizeu), p, t)); Cint(0))
    elseif typeof(prob.u0)<:Vector{Float64}
        f! = prob.f
    else # Then it's an in-place function on an abstract array
        f! = (du, u, p, t) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p, t);
                               du=vec(du); Cint(0))
    end

    if typeof(prob.problem_type) <: SplitODEProblem

        ### Fix the more general function to Sundials allowed style
        if !isinplace && typeof(prob.u0)<:Number
            f1! = (du, u, p, t) -> (du .= prob.f.f1(first(u), p, t); Cint(0))
            f2! = (du, u, p, t) -> (du .= prob.f.f2(first(u), p, t); Cint(0))
        elseif !isinplace && typeof(prob.u0)<:Vector{Float64}
            f1! = (du, u, p, t) -> (du .= prob.f.f1(u, p, t); Cint(0))
            f2! = (du, u, p, t) -> (du .= prob.f.f2(u, p, t); Cint(0))
        elseif !isinplace && typeof(prob.u0)<:AbstractArray
            f1! = (du, u, p, t) -> (du .= vec(prob.f.f1(reshape(u, sizeu), p, t)); Cint(0))
            f2! = (du, u, p, t) -> (du .= vec(prob.f.f2(reshape(u, sizeu), p, t)); Cint(0))
        elseif typeof(prob.u0)<:Vector{Float64}
            f1! = prob.f.f1
            f2! = prob.f.f2
        else # Then it's an in-place function on an abstract array
            f1! = (du, u, p, t) -> (prob.f.f1(reshape(du, sizeu), reshape(u, sizeu), p, t);
                                   du=vec(du); Cint(0))
            f2! = (du, u, p, t) -> (prob.f.f2(reshape(du, sizeu), reshape(u, sizeu), p, t);
                                  du=vec(du); Cint(0))
        end

        userfun = FunJac(f1!,f2!,prob.f.f1.jac,prob.p,prob.f.mass_matrix,
                         prob.f.f1.jac_prototype,alg.prec,alg.psetup,u0,_u0,nothing)

        function getcfunjac(userfun::T) where T
            @cfunction(cvodefunjac, Cint,
                     (realtype, N_Vector,
                     N_Vector, Ref{T}))
        end
        function getcfunjac2(userfun::T) where T
            @cfunction(cvodefunjac2, Cint,
                     (realtype, N_Vector,
                     N_Vector, Ref{T}))
        end
        cfj1 = getcfunjac(userfun)
        cfj2 = getcfunjac2(userfun)

        flag = ARKodeInit(mem,
                    cfj1,cfj2,
                    t0, convert(N_Vector, u0nv))
    else
        userfun = FunJac(f!,prob.f.jac,prob.p,prob.f.mass_matrix,prob.f.jac_prototype,alg.prec,alg.psetup,u0,_u0)
        if alg.stiffness == Explicit()
            function getcfun1(userfun::T) where T
                @cfunction(cvodefunjac, Cint,
                         (realtype, N_Vector,
                         N_Vector, Ref{T}))
            end
            cfj1 = getcfun1(userfun)
            flag = ARKodeInit(mem,
                        cfj1,
                        C_NULL,
                        t0, convert(N_Vector, u0nv))
        elseif alg.stiffness == Implicit()
            function getcfun2(userfun::T) where T
                @cfunction(cvodefunjac, Cint,
                         (realtype, N_Vector,
                         N_Vector, Ref{T}))
            end
            cfj2 = getcfun2(userfun)
            flag = ARKodeInit(mem,
                        C_NULL,cfj2,
                        t0, convert(N_Vector, u0nv))
        end
    end

    dt != nothing && (flag = ARKodeSetInitStep(mem, dt))
    flag = ARKodeSetMinStep(mem, dtmin)
    flag = ARKodeSetMaxStep(mem, dtmax)
    flag = ARKodeSetUserData(mem, userfun)
    if typeof(abstol) <: Array
        flag = ARKodeSVtolerances(mem, reltol, abstol)
    else
        flag = ARKodeSStolerances(mem, reltol, abstol)
    end
    flag = ARKodeSetMaxNumSteps(mem, maxiters)
    flag = ARKodeSetMaxHnilWarns(mem, alg.max_hnil_warns)
    flag = ARKodeSetMaxErrTestFails(mem, alg.max_error_test_failures)
    flag = ARKodeSetMaxNonlinIters(mem, alg.max_nonlinear_iters)
    flag = ARKodeSetMaxConvFails(mem, alg.max_convergence_failures)
    flag = ARKodeSetPredictorMethod(mem, alg.predictor_method)
    flag = ARKodeSetNonlinConvCoef(mem, alg.nonlinear_convergence_coefficient)
    flag = ARKodeSetDenseOrder(mem,alg.dense_order)

    if alg.itable == nothing && alg.etable == nothing
        flag = ARKodeSetOrder(mem,alg.order)
    elseif alg.itable == nothing && alg.etable != nothing
        flag = ARKodeSetERKTableNum(mem, alg.etable)
    elseif alg.itable != nothing && alg.etable == nothing
        flag = ARKodeSetIRKTableNum(mem, alg.itable)
    else
        flag = ARKodeSetARKTableNum(mem, alg.itable, alg.etable)
    end

    flag = ARKodeSetNonlinCRDown(mem,alg.crdown)
    flag = ARKodeSetNonlinRDiv(mem, alg.rdiv)
    flag = ARKodeSetDeltaGammaMax(mem, alg.dgmax)
    flag = ARKodeSetMaxStepsBetweenLSet(mem, alg.msbp)
    #flag = ARKodeSetAdaptivityMethod(mem,alg.adaptivity_method,1,0)

    #flag = ARKodeSetFixedStep(mem,)
    alg.set_optimal_params && (flag = ARKodeSetOptimalParams(mem))

    if Method == :Newton # Only use a linear solver if it's a Newton-based method
        if LinearSolver == :Dense
            A = SUNDenseMatrix(length(u0),length(u0))
            LS = SUNDenseLinearSolver(u0,A)
            flag = ARKDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,DenseMatrix())
            _LS = LinSolHandle(LS,Dense())
        elseif LinearSolver == :Band
            A = SUNBandMatrix(length(u0), alg.jac_upper, alg.jac_lower,
                                       alg.stored_upper)
            LS = SUNBandLinearSolver(u0,A)
            flag = ARKDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,BandMatrix())
            _LS = LinSolHandle(LS,Band())
        elseif LinearSolver == :GMRES
            LS = SUNSPGMR(u0, alg.prec_side, alg.krylov_dim)
            flag = ARKSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = Sundials.LinSolHandle(LS,Sundials.SPGMR())
        elseif LinearSolver == :FGMRES
            LS = SUNSPFGMR(u0, alg.prec_side, alg.krylov_dim)
            flag = ARKSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,SPFGMR())
        elseif LinearSolver == :BCG
            LS = SUNSPBCGS(u0, alg.prec_side, alg.krylov_dim)
            flag = ARKSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,SPBCGS())
        elseif LinearSolver == :PCG
            LS = SUNPCG(u0, alg.prec_side, alg.krylov_dim)
            flag = ARKSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,PCG())
        elseif LinearSolver == :TFQMR
            LS = SUNSPTFQMR(u0, alg.prec_side, alg.krylov_dim)
            flag = ARKSpilsSetLinearSolver(mem, LS)
            _A = nothing
            _LS = LinSolHandle(LS,PTFQMR())
        elseif LinearSolver == :KLU
            nnz = length(nonzeros(prob.f.jac_prototype))
            A = SUNSparseMatrix(length(u0),length(u0), nnz, CSC_MAT)
            LS = SUNKLU(u0, A)
            flag = ARKDlsSetLinearSolver(mem, LS, A)
            _A = MatrixHandle(A,SparseMatrix())
            _LS = LinSolHandle(LS,KLU())
        end
    elseif Method == :Functional
        ARKodeSetFixedPoint(mem, Clong(alg.krylov_dim))
    else
        _A = nothing
        _LS = nothing
    end

    if (typeof(prob.problem_type) <: SplitODEProblem &&
       typeof(prob.f.f1.jac_prototype) <: DiffEqBase.AbstractDiffEqLinearOperator) ||
       (!(typeof(prob.problem_type) <: SplitODEProblem) &&
       typeof(prob.f.jac_prototype) <: DiffEqBase.AbstractDiffEqLinearOperator)
       function getcfunjtimes(userfun::T) where T
           @cfunction(jactimes,
                           Cint,
                           (N_Vector,
                            N_Vector,
                            realtype,
                            N_Vector,
                            N_Vector,
                            Ref{T},
                            N_Vector))
        end
        jtimes = getcfunjtimes(userfun)
        ARKodeSpilsSetJacTimes(mem, C_NULL, jtimes)
    end

    if prob.f.mass_matrix != I
        if MassLinearSolver == :Dense
            M = SUNDenseMatrix(length(u0),length(u0))
            MLS = SUNDenseLinearSolver(u0,M)
            ARKDlsSetMassLinearSolver(mem,MLS,M,false)
            _M = MatrixHandle(M,DenseMatrix())
            _MLS = LinSolHandle(MLS,Dense())
        elseif MassLinearSolver == :Band
            M = SUNBandMatrix(length(u0), alg.mass_upper, alg.mass_lower,
                                       alg.mass_stored_upper)
            MLS = SUNBandLinearSolver(u0,M)
            ARKDlsSetMassLinearSolver(mem,MLS,M,false)
            _M = MatrixHandle(M,BandMatrix())
            _MLS = LinSolHandle(MLS,Band())
        elseif MassLinearSolver == :GMRES
            MLS = SUNSPGMR(u0, alg.prec_side, alg.mass_krylov_dim)
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
            _M = nothing
            _MLS = LinSolHandle(MLS,SPGMR())
        elseif MassLinearSolver == :FGMRES
            MLS = SUNSPGMR(u0, alg.prec_side, alg.mass_krylov_dim)
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
            _M = nothing
            _MLS = LinSolHandle(MLS,SPFGMR())
        elseif MassLinearSolver == :BCG
            MLS = SUNSPGMR(u0, alg.prec_side, alg.mass_krylov_dim)
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
            _M = nothing
            _MLS = LinSolHandle(MLS,SPBCGS())
        elseif MassLinearSolver == :PCG
            MLS = SUNSPGMR(u0, alg.prec_side, alg.mass_krylov_dim)
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
            _M = nothing
            _MLS = LinSolHandle(MLS,PCG())
        elseif MassLinearSolver == :TFQMR
            MLS = SUNSPGMR(u0, alg.prec_side, alg.mass_krylov_dim)
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
            _M = nothing
            _MLS = LinSolHandle(MLS,PTFQMR())
        elseif MassLinearSolver == :KLU
            nnz = length(nonzeros(prob.f.mass_matrix))
            M = SUNSparseMatrix(length(u0),length(u0), nnz, CSC_MAT)
            MLS = SUNKLU(u0, M)
            ARKDlsSetMassLinearSolver(mem,MLS,M,false)
            _M = MatrixHandle(M,SparseMatrix())
            _MLS = LinSolHandle(MLS,KLU())
        else
            ARKSpilsSetMassLinearSolver(mem,MLS,false)
        end
        function getmatfun(userfun::T) where T
            @cfunction(massmat,
                            Cint,
                            (realtype,
                             SUNMatrix,
                             Ref{T},
                             N_Vector,
                             N_Vector,
                             N_Vector))
        end
        matfun = getmatfun(userfun)
        ARKDlsSetMassFn(mem,matfun)
    else
        _M = nothing
        _MLS = nothing
    end

    if DiffEqBase.has_jac(prob.f)
      function getfunjac(userfun::T) where T
          @cfunction(cvodejac,
                          Cint,
                          (realtype,
                           N_Vector,
                           N_Vector,
                           SUNMatrix,
                           Ref{T},
                           N_Vector,
                           N_Vector,
                           N_Vector))
      end
      jac = getfunjac(userfun)
      flag = ARKodeSetUserData(mem, userfun)
      flag = ARKDlsSetJacFn(mem, jac)
    else
        jac = nothing
    end

    if alg.prec !== nothing
        function getpercfun(userfun::T) where T
            @cfunction(precsolve,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             N_Vector,
                             N_Vector,Float64,Float64,Int,
                             Ref{T}))
        end
        precfun = getpercfun(userfun)

        function getpsetupfun(userfun::T) where T
            @cfunction(precsetup,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             Int,
                             Ptr{Int},Float64,Ref{T}))
        end
        psetupfun = alg.psetup === nothing ? C_NULL : getpsetupfun(userfun)

        ARKSpilsSetPreconditioner(mem, psetupfun, precfun)
    end

    callbacks_internal == nothing ? tmp = nothing : tmp = similar(u0)
    callbacks_internal == nothing ? uprev = nothing : uprev = similar(u0)
    tout = [tspan[1]]

    if save_start
      save_value!(ures,u0,uType,sizeu)
      if dense
        f!(_u0,u0,prob.p,tspan[1])
        save_value!(dures,utmp,uType,sizeu)
      end
    end

    sol = DiffEqBase.build_solution(prob, alg, ts, ures,
                   dense = dense,
                   interp = dense ? DiffEqBase.HermiteInterpolation(ts,ures,dures) :
                                    DiffEqBase.LinearInterpolation(ts,ures),
                   timeseries_errors = timeseries_errors,
                   destats = DiffEqBase.DEStats(0),
                   calculate_error = false)
    opts = DEOptions(saveat_internal,tstops_internal,save_everystep,dense,
                     timeseries_errors,dense_errors,save_on,save_end,
                     callbacks_internal,abstol,reltol,verbose,advance_to_tstop,stop_at_next_tstop,
                     progress,progress_name,progress_message)
    integrator = ARKODEIntegrator(utmp,prob.p,t0,t0,mem,_LS,_A,_MLS,_M,sol,alg,f!,userfun,jac,opts,
                       tout,tdir,sizeu,false,tmp,uprev,Cint(flag),false,0,1,callback_cache,0.)

    initialize_callbacks!(integrator)
    integrator
end # function solve

function tstop_saveat_disc_handling(tstops,saveat,tdir,tspan,tType)

  if isempty(tstops) # TODO: Specialize more
    tstops_vec = [tspan[2]]
  else
    tstops_vec = vec(collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x≤tdir*tspan[end],Iterators.flatten((tstops,tspan[end])))))
  end

  if tdir>0
    tstops_internal = BinaryMinHeap(tstops_vec)
  else
    tstops_internal = BinaryMaxHeap(tstops_vec)
  end

  if typeof(saveat) <: Number
    if (tspan[1]:saveat:tspan[end])[end] == tspan[end]
      saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:tspan[end]))
    else
      saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat)))
    end
  elseif isempty(saveat)
    saveat_vec = saveat
  else
    saveat_vec = vec(collect(tType,Iterators.filter(x->tdir*tspan[1]<tdir*x<tdir*tspan[end],saveat)))
  end

  if tdir>0
    saveat_internal = BinaryMinHeap(saveat_vec)
  else
    saveat_internal = BinaryMaxHeap(saveat_vec)
  end

  tstops_internal,saveat_internal
end

## Solve for DAEs uses IDA

function DiffEqBase.__init(
    prob::DiffEqBase.AbstractDAEProblem{uType, duType, tupType, isinplace},
    alg::SundialsDAEAlgorithm{LinearSolver},
    timeseries=[], ts=[], ks=[];

    verbose=true,
    dt=nothing, dtmax=0.0,
    save_on=true, save_start=true,
    callback=nothing, abstol=1/10^6, reltol=1/10^3,
    saveat=Float64[], tstops=Float64[], maxiters=Int(1e5),
    timeseries_errors=true,
    dense_errors = false,
    save_everystep=isempty(saveat), dense=save_everystep,
    save_timeseries=nothing, save_end = true,
    progress=false,progress_name="ODE",
    progress_message = DiffEqBase.ODE_DEFAULT_PROG_MESSAGE,
    advance_to_tstop = false, stop_at_next_tstop = false,
    userdata=nothing,
    kwargs...) where {uType, duType, tupType, isinplace, LinearSolver}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnida)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction)
            if DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    if typeof(reltol) <: AbstractArray
        error("Sundials only allows scalar reltol.")
    end

    progress && @logmsg(-1,progress_name,_id=_id = :Sundials,progress=0)

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length(callbacks_internal)
    if max_len_cb isa VectorContinuousCallback
      callback_cache = DiffEqBase.CallbackCache(max_len_cb.len,uBottomEltype,uBottomEltype)
    else
      callback_cache = nothing
    end

    tspan = prob.tspan
    t0 = tspan[1]

    tdir = sign(tspan[2]-tspan[1])

    tstops_internal, saveat_internal =
      tstop_saveat_disc_handling(tstops,saveat,tdir,tspan,tType)

    if typeof(prob.u0) <: Number
        u0 = [prob.u0]
    else
        u0 = vec(copy(prob.u0))
    end

    if typeof(prob.du0) <: Number
        du0 = [prob.du0]
    else
        du0 = vec(copy(prob.du0))
    end

    sizeu = size(prob.u0)
    sizedu = size(prob.du0)

    ### Fix the more general function to Sundials allowed style
    if !isinplace && typeof(prob.u0)<:Number
        f! = (out, du, u, p, t) -> (out .= prob.f(first(du),first(u), p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:Vector{Float64}
        f! = (out, du, u, p, t) -> (out .= prob.f(du, u, p, t); Cint(0))
    elseif !isinplace && typeof(prob.u0)<:AbstractArray
        f! = (out, du, u, p, t) -> (out .= vec(
                            prob.f(reshape(du, sizedu), reshape(u, sizeu), p, t)
                                 );Cint(0))
    elseif typeof(prob.u0)<:Vector{Float64}
        f! = prob.f
    else # Then it's an in-place function on an abstract array
        f! = (out, du, u, p, t) -> (prob.f(reshape(out, sizeu), reshape(du, sizedu),
                                    reshape(u, sizeu), p, t); Cint(0))
    end

    mem_ptr = IDACreate()
    (mem_ptr == C_NULL) && error("Failed to allocate IDA solver object")
    mem = Handle(mem_ptr)

    !verbose && IDASetErrHandlerFn(mem,@cfunction(null_error_handler, Nothing,
                                    (Cint, Char,
                                    Char, Ptr{Cvoid})),C_NULL)

    ures = Vector{uType}()
    dures = Vector{uType}()
    ts   = [t0]

    _u0 = copy(u0)
    utmp = NVector(_u0)
    _du0 = copy(du0)
    dutmp = NVector(_du0)
    rtest = zeros(length(u0))

    userfun = FunJac(f!,prob.f.jac,prob.p,nothing,prob.f.jac_prototype,alg.prec,alg.psetup,_u0,_du0,rtest)

    u0nv = NVector(u0)

    function getcfun(userfun::T) where T
        @cfunction(idasolfun,
                         Cint, (realtype, N_Vector, N_Vector,
                                N_Vector, Ref{T}))
    end
    cfun = getcfun(userfun)
    flag = IDAInit(mem, cfun,
                    t0, convert(N_Vector, u0),
                    convert(N_Vector, du0))
    dt != nothing && (flag = IDASetInitStep(mem, dt))
    flag = IDASetUserData(mem, userfun)
    flag = IDASetMaxStep(mem, dtmax)
    if typeof(abstol) <: Array
        flag = IDASVtolerances(mem, reltol, abstol)
    else
        flag = IDASStolerances(mem, reltol, abstol)
    end
    flag = IDASetMaxNumSteps(mem, maxiters)
    flag = IDASetMaxOrd(mem,alg.max_order)
    flag = IDASetMaxErrTestFails(mem,alg.max_error_test_failures)
    flag = IDASetNonlinConvCoef(mem,alg.nonlinear_convergence_coefficient)
    flag = IDASetMaxNonlinIters(mem,alg.max_nonlinear_iters)
    flag = IDASetMaxConvFails(mem,alg.max_convergence_failures)
    flag = IDASetNonlinConvCoefIC(mem,alg.nonlinear_convergence_coefficient_ic)
    flag = IDASetMaxNumStepsIC(mem,alg.max_num_steps_ic)
    flag = IDASetMaxNumJacsIC(mem,alg.max_num_jacs_ic)
    flag = IDASetMaxNumItersIC(mem,alg.max_num_iters_ic)
    #flag = IDASetMaxBacksIC(mem,alg.max_num_backs_ic) # Needs newer version?
    flag = IDASetLineSearchOffIC(mem,alg.use_linesearch_ic)

    if LinearSolver == :Dense
        A = SUNDenseMatrix(length(u0),length(u0))
        LS = SUNDenseLinearSolver(u0,A)
        flag = IDADlsSetLinearSolver(mem, LS, A)
        _A = MatrixHandle(A,DenseMatrix())
        _LS = LinSolHandle(LS,Dense())
    elseif LinearSolver == :Band
        A = SUNBandMatrix(length(u0), alg.jac_upper, alg.jac_lower,
                                   alg.stored_upper)
        LS = SUNBandLinearSolver(u0,A)
        flag = IDADlsSetLinearSolver(mem, LS, A)
        _A = MatrixHandle(A,BandMatrix())
        _LS = LinSolHandle(LS,Band())
    elseif LinearSolver == :GMRES
        LS = SUNSPGMR(u0, alg.prec_side, alg.krylov_dim)
        flag = IDASpilsSetLinearSolver(mem, LS)
        _A = nothing
        _LS = LinSolHandle(LS,SPGMR())
    elseif LinearSolver == :FGMRES
        LS = SUNSPFGMR(u0, alg.prec_side, alg.krylov_dim)
        flag = IDASpilsSetLinearSolver(mem, LS)
        _A = nothing
        _LS = LinSolHandle(LS,SPFGMR())
    elseif LinearSolver == :BCG
        LS = SUNSPBCGS(u0, alg.prec_side, alg.krylov_dim)
        flag = IDASpilsSetLinearSolver(mem, LS)
        _A = nothing
        _LS = LinSolHandle(LS,SPBCGS())
    elseif LinearSolver == :PCG
        LS = SUNPCG(u0, alg.prec_side, alg.krylov_dim)
        flag = IDASpilsSetLinearSolver(mem, LS)
        _A = nothing
        _LS = LinSolHandle(LS,PCG())
    elseif LinearSolver == :TFQMR
        LS = SUNSPTFQMR(u0, alg.prec_side, alg.krylov_dim)
        flag = IDASpilsSetLinearSolver(mem, LS)
        _A = nothing
        _LS = LinSolHandle(LS,PTFQMR())
    elseif LinearSolver == :KLU
        nnz = length(nonzeros(prob.f.jac_prototype))
        A = SUNSparseMatrix(length(u0),length(u0), nnz, CSC_MAT)
        LS = SUNKLU(u0, A)
        flag = IDADlsSetLinearSolver(mem, LS, A)
        _A = MatrixHandle(A,SparseMatrix())
        _LS = LinSolHandle(LS,KLU())
    end

    if typeof(prob.f.jac_prototype) <: DiffEqBase.AbstractDiffEqLinearOperator
        function getcfunjtimes(userfun::T) where T
            @cfunction(idajactimes,
                           Cint,
                           (realtype,
                            N_Vector,N_Vector,N_Vector,N_Vector,N_Vector,
                            realtype,
                            Ref{T},
                            N_Vector,N_Vector))
        end
        jtimes = getcfunjtimes(userfun)
        IDASpilsSetJacTimes(mem, C_NULL, jtimes)
    end

    if alg.prec !== nothing
        function getpercfun(userfun::T) where T
            @cfunction(idaprecsolve,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             N_Vector,
                             N_Vector,N_Vector,Float64,Float64,Int,
                             Ref{T}))
        end
        precfun = getpercfun(userfun)

        function getpsetupfun(userfun::T) where T
            @cfunction(idaprecsetup,
                            Cint,
                            (Float64,
                             N_Vector,
                             N_Vector,
                             N_Vector,
                             Float64,Ref{T}))
        end
        psetupfun = alg.psetup === nothing ? C_NULL : getpsetupfun(userfun)

        IDASpilsSetPreconditioner(mem, psetupfun, precfun)
    end

    if DiffEqBase.has_jac(prob.f)
      function getcfunjacc(userfun::T) where T
          @cfunction(idajac,
                     Cint,
                     (realtype,
                      realtype,
                      N_Vector,
                      N_Vector,
                      N_Vector,
                      SUNMatrix,
                      Ref{T},
                      N_Vector,
                      N_Vector,
                      N_Vector))
      end
      jac = getcfunjacc(userfun)
      flag = IDASetUserData(mem, userfun)
      flag = IDADlsSetJacFn(mem, jac)
    else
      jac = nothing
    end

    tout = [tspan[1]]

    f!(rtest, du0, u0, prob.p, t0)
    if any(abs.(rtest) .>= reltol)
        if prob.differential_vars === nothing && !alg.init_all
            error("Must supply differential_vars argument to DAEProblem constructor to use IDA initial value solver.")
        end
        prob.differential_vars != nothing && (flag = IDASetId(mem, collect(Float64, prob.differential_vars)))

        if dt != nothing
            _t = float(dt)
        else
            _t = float(tspan[2])
        end
        if alg.init_all
            init_type = IDA_Y_INIT
        else
            init_type = IDA_YA_YDP_INIT
        end
        flag = IDACalcIC(mem, init_type, _t)
    end

    if save_start
      save_value!(ures,u0,uType,sizeu)
      if dense
        save_value!(dures,du0,uType,sizedu) # Does this need to update for IDACalcIC?
      end
    end

    callbacks_internal == nothing ? tmp = nothing : tmp = similar(u0)
    callbacks_internal == nothing ? uprev = nothing : uprev = similar(u0)

    if flag >= 0
        retcode = :Default
    else
        retcode = :InitialFailure
    end

    sol = DiffEqBase.build_solution(prob, alg, ts, ures,
                   dense = dense,
                   interp = dense ? DiffEqBase.HermiteInterpolation(ts,ures,dures) :
                                    DiffEqBase.LinearInterpolation(ts,ures),
                   calculate_error = false,
                   timeseries_errors = timeseries_errors,
                   retcode = retcode,
                   destats = DiffEqBase.DEStats(0),
                   dense_errors = dense_errors)

    opts = DEOptions(saveat_internal,tstops_internal,save_everystep,dense,
                    timeseries_errors,dense_errors,save_on,save_end,
                    callbacks_internal,abstol,reltol,verbose,advance_to_tstop,stop_at_next_tstop,
                    progress,progress_name,progress_message)

    integrator = IDAIntegrator(utmp,dutmp,prob.p,t0,t0,mem,_LS,_A,sol,alg,f!,userfun,jac,opts,
                   tout,tdir,sizeu,sizedu,false,tmp,uprev,Cint(flag),false,0,1,callback_cache,0.)

    initialize_callbacks!(integrator)
    integrator
end # function solve

## Common calls

function interpret_sundials_retcode(flag)
  flag >= 0 && return :Success
  flag == -1 && return :MaxIters
  (flag == -2 || flag == -3) && return :Unstable
  flag == -4 && return :ConvergenceFailure
  return :Failure
end

function solver_step(integrator::CVODEIntegrator,tstop)
    integrator.flag = CVode(integrator.mem, tstop, integrator.u, integrator.tout, CV_ONE_STEP)
    if integrator.opts.progress
      @logmsg(-1,
      integrator.opts.progress_name,
      _id = :Sundials,
      message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
      progress=integrator.t/integrator.sol.prob.tspan[2])
    end
end
function solver_step(integrator::ARKODEIntegrator,tstop)
    integrator.flag = ARKode(integrator.mem, tstop, integrator.u, integrator.tout, ARK_ONE_STEP)
    if integrator.opts.progress
      @logmsg(-1,
      integrator.opts.progress_name,
      _id = :Sundials,
      message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
      progress=integrator.t/integrator.sol.prob.tspan[2])
    end
end
function solver_step(integrator::IDAIntegrator,tstop)
    integrator.flag = IDASolve(integrator.mem, tstop, integrator.tout,
                               integrator.u, integrator.du, IDA_ONE_STEP)
    if integrator.opts.progress
      @logmsg(-1,
      integrator.opts.progress_name,
      _id = :Sundials,
      message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
      progress=integrator.t/integrator.sol.prob.tspan[2])
    end
end

function set_stop_time(integrator::CVODEIntegrator,tstop)
    CVodeSetStopTime(integrator.mem,tstop)
end
function set_stop_time(integrator::ARKODEIntegrator,tstop)
    ARKodeSetStopTime(integrator.mem,tstop)
end
function set_stop_time(integrator::IDAIntegrator,tstop)
    IDASetStopTime(integrator.mem,tstop)
end

function DiffEqBase.solve!(integrator::AbstractSundialsIntegrator)
    uType = eltype(integrator.sol.u)
    while !isempty(integrator.opts.tstops)
        # Sundials can have floating point issues approaching a tstop if
        # there is a modifying event each
        while integrator.tdir*(integrator.t-top(integrator.opts.tstops)) < -1e6eps()
            tstop = top(integrator.opts.tstops)
            set_stop_time(integrator,tstop)
            integrator.tprev = integrator.t
            if !(typeof(integrator.opts.callback.continuous_callbacks)<:Tuple{})
                integrator.uprev .= integrator.u
            end
            integrator.userfun.p = integrator.p
            solver_step(integrator,tstop)
            integrator.t = first(integrator.tout)
            integrator.flag < 0 && break
            handle_callbacks!(integrator)
            integrator.flag < 0 && break
            if isempty(integrator.opts.tstops)
              break
            end
        end
        (integrator.flag < 0) && break
        handle_tstop!(integrator)
    end

    if integrator.opts.save_end && integrator.sol.t[end] != integrator.t
        save_value!(integrator.sol.u,integrator.u,uType,integrator.sizeu)
        push!(integrator.sol.t, integrator.t)
        if integrator.opts.dense
          integrator(integrator.u,integrator.t,Val{1})
          save_value!(integrator.sol.interp.du,integrator.u,uType,integrator.sizeu)
        end
    end

    if integrator.opts.progress
      @logmsg(-1,
      integrator.opts.progress_name,
      _id = :Sundials,
      message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
      progress="done")
    end

    fill_destats!(integrator)
    empty!(integrator.mem)
    integrator.A != nothing && empty!(integrator.A)
    integrator.LS != nothing && empty!(integrator.LS)

    if DiffEqBase.has_analytic(integrator.sol.prob.f)
        DiffEqBase.calculate_solution_errors!(integrator.sol;
        timeseries_errors=integrator.opts.timeseries_errors,
        dense_errors=integrator.opts.dense_errors)
    end

    if integrator.sol.retcode != :Default
      return integrator.sol
    end
    integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol,interpret_sundials_retcode(integrator.flag))
    nothing
end

function handle_tstop!(integrator::AbstractSundialsIntegrator)
    tstops = integrator.opts.tstops
    if !isempty(tstops)
      if integrator.tdir*(integrator.t-top(integrator.opts.tstops)) > -1e6eps()
          pop!(tstops)
          t = integrator.t
          integrator.just_hit_tstop = true
      end
    end
end

function fill_destats!(integrator::AbstractSundialsIntegrator)
end

function fill_destats!(integrator::CVODEIntegrator)
    destats = integrator.sol.destats
    mem = integrator.mem
    tmp = Ref(Clong(-1))
    CVodeGetNumRhsEvals(mem,tmp)
    destats.nf = tmp[]
    CVodeGetNumLinSolvSetups(mem,tmp)
    destats.nw = tmp[]
    CVodeGetNumErrTestFails(mem,tmp)
    destats.nreject = tmp[]
    CVodeGetNumSteps(mem,tmp)
    destats.naccept = tmp[] - destats.nreject
    CVodeGetNumNonlinSolvIters(mem,tmp)
    destats.nnonliniter = tmp[]
    CVodeGetNumNonlinSolvConvFails(mem,tmp)
    destats.nnonlinconvfail = tmp[]
    if method_choice(integrator.alg) == :Newton
        CVDlsGetNumJacEvals(mem,tmp)
        destats.njacs = tmp[]
    end
end

function fill_destats!(integrator::ARKODEIntegrator)
    destats = integrator.sol.destats
    mem = integrator.mem
    tmp = Ref(Clong(-1))
    tmp2 = Ref(Clong(-1))
    ARKodeGetNumRhsEvals(mem,tmp,tmp2)
    destats.nf = tmp[]
    destats.nf2 = tmp2[]
    ARKodeGetNumLinSolvSetups(mem,tmp)
    destats.nw = tmp[]
    ARKodeGetNumErrTestFails(mem,tmp)
    destats.nreject = tmp[]
    ARKodeGetNumSteps(mem,tmp)
    destats.naccept = tmp[] - destats.nreject
    ARKodeGetNumNonlinSolvIters(mem,tmp)
    destats.nnonliniter = tmp[]
    ARKodeGetNumNonlinSolvConvFails(mem,tmp)
    destats.nnonlinconvfail = tmp[]
    if method_choice(integrator.alg) == :Newton
        ARKDlsGetNumJacEvals(mem,tmp)
        destats.njacs = tmp[]
    end
end

function fill_destats!(integrator::IDAIntegrator)
    destats = integrator.sol.destats
    mem = integrator.mem
    tmp = Ref(Clong(-1))
    IDAGetNumResEvals(mem,tmp)
    destats.nf = tmp[]
    IDAGetNumLinSolvSetups(mem,tmp)
    destats.nw = tmp[]
    IDAGetNumErrTestFails(mem,tmp)
    destats.nreject = tmp[]
    IDAGetNumSteps(mem,tmp)
    destats.naccept = tmp[] - destats.nreject
    IDAGetNumNonlinSolvIters(mem,tmp)
    destats.nnonliniter = tmp[]
    IDAGetNumNonlinSolvConvFails(mem,tmp)
    destats.nnonlinconvfail = tmp[]
    if method_choice(integrator.alg) == :Newton
        IDADlsGetNumJacEvals(mem,tmp)
        destats.njacs = tmp[]
    end
end

function initialize_callbacks!(integrator, initialize_save = true)
  t = integrator.t
  u = integrator.u
  callbacks = integrator.opts.callback
  integrator.u_modified = true

  u_modified = initialize!(callbacks,u,t,integrator)

  # if the user modifies u, we need to fix current values
  if u_modified

    handle_callback_modifiers!(integrator)

    if initialize_save &&
      (any((c)->c.save_positions[2],callbacks.discrete_callbacks) ||
      any((c)->c.save_positions[2],callbacks.continuous_callbacks))
      savevalues!(integrator,true)
    end
  end

  # reset this as it is now handled so the integrators should proceed as normal
  integrator.u_modified = false
end


text:
# ### Load packages
using GeoPhyInv

# ### A surface acquisition
mgrid=repeat(range(-1000.0,stop=1000.0,length=201), 2) # test spatial grid
#ageom_surf=GeoPhyInv.Gallery.AGeom(mgrid,:surf,nss=5,nr=30);

#

ext: jl


text:
const function2latex = Dict{Symbol, String}(
    ## Greek alphabet
    :alpha => "\\alpha",
    :beta => "\\beta",
    :gamma => "\\gamma",
    :delta => "\\delta",
    :epsilon => "\\epsilon",
    :zeta => "\\zeta",
    :eta => "\\eta",
    :theta => "\\theta",
    :iota => "\\iota",
    :kappa => "\\kappa",
    :lambda => "\\lambda",
    :mu => "\\mu",
    :nu => "\\nu",
    :xi => "\\xi",
    :pi => "\\pi",
    :rho => "\\rho",
    :sigma => "\\sigma",
    :tau => "\\tau",
    :upsilon => "\\upsilon",
    :phi => "\\phi",
    :chi => "\\chi",
    :psi => "\\psi",
    :omega => "\\omega",
    :Gamma => "\\Gamma",
    :Delta => "\\Delta",
    :Theta => "\\Theta",
    :Lambda => "\\Lambda",
    :Xi => "\\Xi",
    :Pi => "\\Pi",
    :Sigma => "\\Sigma",
    :Upsilon => "\\Upsilon",
    :Phi => "\\Phi",
    :Psi => "\\Psi",
    :Omega => "\\Omega",
    ## Trinogometry
	:sin =>   "\\mathrm{sen}",
    :cos =>   "\\mathrm{cos}",
    :tan =>   "\\mathrm{tg}",
	:sind =>  "\\mathrm{sen}",
    :cosd =>  "\\mathrm{cos}",
    :tand =>  "\\mathrm{tg}",
    :cot =>   "\\mathrm{cotg}",
    :sec =>   "\\mathrm{sec}",
    :csc =>   "\\mathrm{cos\\,sec}",
	:sinh =>  "\\mathrm{senh}",
    :cosh =>  "\\mathrm{cosh}",
    :tanh =>  "\\mathrm{tanh}",
    :coth =>  "\\mathrm{coth}",
	:asin =>  "\\mathrm{arc\\,sen}",
    :acos =>  "\\mathrm{arc\\,cos}",
    :atan =>  "\\mathrm{arc\\,tg}",
    :atan2 => "\\mathrm{atan2}",
    :asinh => "\\mathrm{arcsinh}",
    :sinc =>  "\\mathrm{senc}",
    :acosh => "\\mathrm{arc\\,cosh}",
    :cosc =>  "\\mathrm{cosc}",
    :atanh => "\\mathrm{arc\\,tanh}",
    :acot =>  "\\mathrm{arc\\,cot}",
    :acoth => "\\mathrm{arc\\,coth}",
    :asec =>  "\\mathrm{arc\\,sec}",
    :sech =>  "\\mathrm{sech}",
    :asech => "\\mathrm{arc\\,sech}",
    :acsc =>  "\\mathrm{arc\\,csc}",
    :csch =>  "\\mathrm{csch}",
    :acsch => "\\mathrm{arc\\,csch}",
    ## Misc
	:log =>   "\\mathrm{ln}",
    :log10 => "\\log_{10}",
    :log2 =>   "\\log_{2}",
    :gamma => "\\Gamma", # The Gamma function
    )

const trigonometric_functions = [
    :sin,
    :cos,
    :tan,
	:sind,
	:cosd,
	:tand,
    :cot,
    :sec,
    :csc,
    :sinh,
    :cosh,
    :tanh,
    :coth,
    :asin,
    :acos,
    :atan,
    :atan2,
    :asinh,
    :sinc,
    :acosh,
    :cosc,
    :atanh,
    :acot,
    :acoth,
    :asec,
    :sech,
    :asech,
    :acsc,
    :csch,
    :acsch,
    ]

const trigonometric_functions_degrees = [
	:sind,
	:cosd,
	:tand,
	]


text:
#TODO keep max/min, prevent relim

using Munkres
using PyCall
pygui(:tk) #prevent conflict with gtk
using PyPlot
# remove all keybindings
D = PyCall.PyDict(matplotlib["rcParams"])
for k in keys(D); contains(k, "keymap") && (D[k] = Any[]) end
ioff()

#####

# this visualization module expects functions specific to Galerkin (projection)

"""
todo
"""
type GalerkinViz <: Visualization
	parent::Session
	figBif
	figSol
	idxLineBranch::Dict
	idxBranchLines::Dict
	color::Dict
	activeSolutionMark
end

function Base.show(V::GalerkinViz)
	figure(V.figBif[:number])
	clf()
	PyPlot.draw()
	V.idxBranchLines = Dict{Any,Any}()
	V.idxLineBranch = Dict{Any,Any}()
	map(b -> try pushBranch(V,b) end, V.parent.P)
	V.figBif["show"]()
	V.figSol["show"]()
end

function GalerkinViz(parent::Session)
	figBif = figure(figsize=(8,6), dpi=80, facecolor="w")
	figBif["canvas"]["set_window_title"]("Bifurcation Plot")
	subplot(111)
	subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1)

	figSol = figure(figsize=(4,3), dpi=80, facecolor="w")
	figSol["canvas"]["set_window_title"]("Solution Plot")
	subplot(111, projection="3d")
	subplots_adjust(left=0.01, right=0.99, top=0.99, bottom=0.01)
	figSol[:canvas][:toolbar][:pack_forget]()

	V = GalerkinViz(parent, figBif, figSol, Dict{Any,Any}(), Dict{Any,Any}(), Dict{Any,Any}(), Void)

	figBif["canvas"]["mpl_connect"]("pick_event", ev -> handlerPick(V,ev))
	figBif["canvas"]["mpl_connect"]("key_press_event", ev -> @schedule handlerKey(V,ev))

	#connect to proj events
	P = parent.P

	observe(P, :pushBranch) do b
		pushBranch(V, b) end

	observe(P, :activeSolutionChanged) do
		s = V.parent.P.activeSolution #TODO no direct access...
		isa(s, Solution) && plotSolution(V, s.data)
		isa(s, Vector{Float64}) && plotSolution(V, s)
	end

	observe(P, :pushSolution) do b,s
		addToBranch(push!, V, b, s) end

	observe(P, :unshiftSolution) do b,s
		addToBranch(unshift!, V, b, s) end

	observe(P, :popSolution) do b,s
		delFromBranch(pop!, V, b) end

	observe(P, :shiftSolution) do b,s
		delFromBranch(shift!, V, b) end

	observe(P, :delBranch) do b
		delBranch(V, b) end

	return V
end

@noinline function handlerPick(V::GalerkinViz, ev) b = get(V.idxLineBranch, hash(ev[:artist]), nothing)
	b == nothing && return

	try
		i = ev[:ind][1]+1
		setActiveSolution(V.parent.P, b[i])
	catch e
		println(e)
	end
	return Void
end

@noinline function handlerKey(V::GalerkinViz, ev)
	# println(ev[:key]) #debug

	figure(V.figBif[:number])

	if ev[:key] == "backspace"
		!isa(V.parent.P.activeSolution, Solution) && return Void
		s,b = V.parent.P.activeSolution, V.parent.P.activeSolution.parent
		if length(b) ≤ 2
			delBranch(V, b)
			deleteat!(V.parent.P, findfirst(V.parent.P, b)) #TODO fix
		elseif s==b[end]
			setActiveSolution(V.parent.P, b[end-1])
			pop!(V.parent.P.activeSolution.parent)
		elseif s==b[1]
			setActiveSolution(V.parent.P, b[2])
			shift!(V.parent.P.activeSolution.parent)
		end
	elseif ev[:key] == "ctrl+delete" #TODO fix
		!isa(V.parent.P.activeSolution, Solution) && return Void
		b = V.parent.P.activeSolution.parent
		delBranch(V, b)
		deleteat!(V.parent.P, findfirst(V.parent.P, b)) #TODO fix
	elseif ev[:key] == "r"
		figure(V.figBif[:number])
		gca()[:relim]()
		autoscale()
	elseif ev[:key] == "ctrl+r"
		# figure(V.figBif[:number])
		# clf()
		# PyPlot.draw()
		# sleep(.1)
		# V.idxBranchLines = Dict{Any,Any}()
		# V.idxLineBranch = Dict{Any,Any}()
		# map(b -> pushBranch(V,b), V.parent.P)
		show(V)
	elseif ev[:key] == "l"
		!isa(V.parent.P.activeSolution, Solution) && return Void
		b = V.parent.P.activeSolution.parent
		for l in V.idxBranchLines[b]
			if l[:get_picker]() != 5
				l[:set_picker](5)
				l[:set_marker](".")
			else
				l[:set_picker](0)
				l[:set_marker]("None")
			end
		end
	elseif ev[:key] == "ctrl+l"
		for L in values(V.idxBranchLines)
			for l in L
				l[:set_picker](0)
				l[:set_marker]("None")
			end
		end
	elseif ev[:key] == "ctrl+L"
		for L in values(V.idxBranchLines)
			for l in L
				l[:set_picker](5)
				l[:set_marker](".")
			end
		end
	elseif ev[:key] == "ctrl+h"
		!isa(V.parent.P.activeSolution, Solution) && return Void
		b = V.parent.P.activeSolution.parent
		for l in V.idxBranchLines[b]
			l[:set_picker](0)
			l[:set_visible](false)
		end
	elseif ev[:key] == "ctrl+H"
		for L in values(V.idxBranchLines)
			for l in L
				l[:set_picker](5)
				l[:set_visible](true)
			end
		end
	elseif ev[:key] in [ "$(i)" for i in 0:9 ]
		const colMap = Dict(
			"0" => "#000000", "1" => "#FF0000",
			"2" => "#00FF00", "3" => "#0000FF",
			"4" => "#00FFFF", "5" => "#FF00FF",
			"6" => "#FFFF00", "7" => "#800000",
			"8" => "#008000", "9" => "#000080"
			)
		if typeof(V.parent.P.activeSolution) == Solution
			for l in V.idxBranchLines[V.parent.P.activeSolution.parent]
				V.color[V.parent.P.activeSolution.parent] = colMap[ev[:key]]
				l[:set_color](colMap[ev[:key]])
			end
		end
	elseif ev[:key] == "b"
		!isa(V.parent.P.activeSolution, Solution) && return
		b = V.parent.P.activeSolution.parent
		plotBranch(V, b)
	elseif ev[:key] == " "
		step(V.parent.cont)
	elseif ev[:key] == "e"
		!isa(V.parent.P.activeSolution, Solution) && return
		s = V.parent.P.activeSolution
		b = s.parent
		setActiveSolution(V.parent.P, s==b[1] ? b[end] : b[1])
	elseif ev[:key] in ["up","down"]
		P = V.parent.P
		s = P.activeSolution
		!isa(s, Solution) && return
		b = s.parent
		i = findfirst(P,b)-1
		b′ = P[ mod( ev[:key] == "up" ? i-1 : i+1, length(P)) + 1 ]
		setActiveSolution(P, b′[1])
	elseif ev[:key] in ["left","right"]
		P = V.parent.P
		s = P.activeSolution
		!isa(s, Solution) && return
		b = s.parent
		i = findfirst(b,s)-1
		s′ = b[ mod( ev[:key] == "left" ? i-1 : i+1 , length(b)) + 1 ]
		setActiveSolution(P, s′)
	elseif ev[:key] == "d"
		s = V.parent.P.activeSolution
		isa(s, Solution) && plotSolution(V, s.data; split=false)
		isa(s, Vector{Float64}) && plotSolution(V, s; split=false)
	elseif ev[:key] == "D"
		s = V.parent.P.activeSolution
		isa(s, Solution) && plotSolution(V, s.data; split=true)
		isa(s, Vector{Float64}) && plotSolution(V, s; split=true)
	end

	PyPlot.draw()
	return Void
end

function pushBranch(V::GalerkinViz, b::Branch)
	isempty(b) && error("empty branch")

	x = map(last, b)

	function reducer(A,a)
		sA,sa = size(A,2), size(a,2)

		if sA ≠ sa
			warn("in-branch change in number of projected points")
			tA,ta = vec(A[end,:]), vec(a)

			# find repetitions of elements of A s.t. norm(A[is]-B) is minimal
			function align(A,B)
				m,n = length(A), length(B)
				@assert m < n
				d(a,b) = norm(a-b)
				M = fill(Inf,m,n)
				M[1,1] = d(A[1],B[1])
				for j in 2:n-m+1 M[1,j] = M[1,j-1] + d(A[1],B[j]) end
				for i in 2:m, j in i:n-m+i
					M[i,j] = min(M[i-1,j-1], M[i,j-1]) + d(A[i],B[j])
				end
				is = [m]
				for j in n-1:-1:1
					i = is[1]
					unshift!(is, i ≠ 1 ? M[i-1,j] ≤ M[i,j] ? i-1 : i : 1)
				end
				return is
			end

			if sA>sa
				a = ta[align(ta,tA)]'
			else
				A = A[:,align(tA,ta)]
			end
		end

		return vcat(A,a)
	end

	function mapper(x)
		tmp = projection(V.parent.core, x)
		length(tmp) == 0 && warn("empty projection")
		return tmp'
	end

	y = mapreduce(mapper, reducer, b)

	figure(V.figBif[:number])
	hold(true)
	lines = plot(x, y, picker=5, color=get(V.color, b, "k"), marker=".", markersize=3)
	PyPlot.draw()

	map(l -> (V.idxLineBranch[hash(l)]=b), lines)
	V.idxBranchLines[b] = lines

	return lines
end

function addToBranch(op, V::GalerkinViz, b::Branch, s::Solution)
	!haskey(V.idxBranchLines, b) && error("branch missing")
	lines = V.idxBranchLines[b]

	x = last(s) #TODO extend relay?
	Y = projection(V.parent.core, s)

	for (l,y) in zip(lines, Y)
		l[:set_xdata](op(l[:get_xdata](), x))
		l[:set_ydata](op(l[:get_ydata](), y))
	end

	figure(V.figBif[:number])
	PyPlot.draw()
	return Void
end

#TODO handle last solution removed?
function delFromBranch(op, V::GalerkinViz, b::Branch)
	for l in V.idxBranchLines[b]
		x,y = l[:get_xdata](), l[:get_ydata]()
		op(x); op(y)
		l[:set_xdata](x); l[:set_ydata](y)
	end

	figure(V.figBif[:number])
	PyPlot.draw()
	return Void
end

# plotSolution(V::GalerkinViz, S::Solution) = plotSolution(B,S.data) # handled by convert? prob not..
function plotSolution(V::GalerkinViz, v::Vector{Float64}; split=false)
	figure(V.figBif[:number])
	hold(true)
	try V.activeSolutionMark[:remove]() end
	x = v[end]
	y = projection(V.parent.core, v)
	V.activeSolutionMark = scatter(fill(x, size(y)), y; facecolors="none", edgecolors="k", marker="o")
	PyPlot.draw()

	figure(V.figSol[:number])
	cla()
	gca(projection="3d")
	axis("equal")
	hold(false)
	m = length(v)÷6
	t = linspace(0, 2pi, 2m*4)
	w = interp(unwrap(v)[1])(t)
	split && (w[:,3] += linspace(0, maximum(w[:,3])-minimum(w[:,3]), length(w[:,3])))
	plot(w[:,1], w[:,2], w[:,3], color="k")

	hold(true)
	try cbPlotSolution(v)
	catch e; println(e) end

	PyPlot.draw()

	return Void
end

function delBranch(V::GalerkinViz, b::Branch)
	!haskey(V.idxBranchLines, b) && return
	lines = V.idxBranchLines[b]
	delete!(V.idxBranchLines, b)
	map(lines) do l
		l[:remove]()
		delete!(V.idxLineBranch, hash(l))
	end

	figure(V.figSol[:number])
	PyPlot.draw()
end

# plot a whole branch in
function plotBranch(V::GalerkinViz, B::Branch)
	figure(V.figSol[:number])
	clf()
	gca(projection="3d")
	hold(true)

	for v in B
		t = linspace(0, 2pi, length(v)÷3 * 6)
		w = interp(unwrap(v)[1])(t)
		plot(w[:,1], w[:,2], w[:,3], color="k", alpha=.1)
	end

	v = B[1]
	t = linspace(0, 2pi, length(v)÷3 * 6)
	w = interp(unwrap(v)[1])(t)
	plot(w[:,1], w[:,2], w[:,3], color="y")

	v = B[end]
	t = linspace(0, 2pi, length(v)÷3 * 6)
	w = interp(unwrap(v)[1])(t)
	plot(w[:,1], w[:,2], w[:,3], color="r")

	PyPlot.draw()
	return
end

#=
@noinline function hungarian(M)
	#make square
	addrows, addcols = 0,0
	if size(M,2) < size(M,1)
		addcols = size(M,1)-size(M,2)
		M = hcat(M, fill(maximum(M), size(M,1), addcols))
	elseif size(M,1) < size(M,2)
		addrows = size(M,2)-size(M,1)
		M = vcat(M, fill(maximum(M), addrows, size(M,2)))
	end

	L = size(M,1)

	#subtract row/col minima
	mapslices(c->c-minimum(c), M, [1])
	mapslices(r->r-minimum(r), M, [2])

	while true
		println(0)

		#count zeros
		is,js = zeros(L),zeros(L)
		zs = 0
		for i in 1:L, j in 1:L
			if M[i,j] == 0
				is[i] += 1
				js[j] += 1
				zs += 1
			end
		end

		#sorted idx
		seqi, seqj = sortperm(is, rev=true), sortperm(js, rev=true)

		println(1)

		#mapping with definite candidates set
		# mapping = zeros(Int64, L)
		# for k in L:-1:1
		# 	is[seqi[k]] == 1 && (mapping[seqi[k]] = findfirst(M[seqi[k],:], 0))
		# 	js[seqj[k]] == 1 && (mapping[findfirst(M[seqj[k],:]], 0)] = k)
		# 	is[seqi[k]] > 1 && js[seqj[k]] > 1 && break
		# end
		#
		# #check for solution
		# idx = find(x->x==0, mapping)
		# opts = [ find(x->x==0, M[i,:]) for i in idx]
		# K = length(opts)
		# count = zeros(Int64, K)
		# while true
		# 	#construct candidate solution, check
		# 	candidate = map((o,c)->o[c+1], opts, count)
		# 	flagBadCandidate = false
		# 	for c in 1:K
		# 		for d in c+1:K
		# 			candidate[c] == candidate[d] && (flagBadCandidate = true) && break
		# 		end
		# 		flagBadCandidate && break
		# 	end
		#
		# 	#good solution found
		# 	if !flagBadCandidate
		# 		mapping[idx] = candidate
		# 		return mapping
		# 	end
		#
		# 	#prepare next candidate
		# 	k = length(opts)
		# 	while true
		# 		count[k] = (count[k] + 1) % length(opts[k])
		# 		count[k] ≠ 0 && break
		# 		k -= 1
		# 		k == 0 && break
		# 	end
		# 	k == 0 && break
		# end

		#construct lines through zeros
		lis, ljs = [], []
		while zs > 0
			i,j = L,L
			while is[seqi[i]] == 1
				push!(lis, seqi[i])
				push!(lis, seqi[i])
			end
			if is[i] ≥ js[j]
				zs -= is[i]
				is[i] = 0
				seqi = circshift(seqi, [-1])
				push!(lis, i)
				for j in L:-1:1
					M[i,seqj[j]] ≠ 0 && continue
					js[seqj[j]] -= 1
					while j < L && js[seqj[j]] < js[seqj[j+1]]
						seqj[j],seqj[j+1] = seqj[j+1],seqj[j]
						j += 1
					end
				end
			else
				zs -= js[j]
				js[j] = 0
				seqj = circshift(seqj, [-1])
				push!(ljs, j)
				for i in L:-1:1
					M[seqi[i],j] ≠ 0 && continue
					is[seqi[i]] -= 1
					while i < L && is[seqi[i]] < is[seqi[i+1]]
						seqi[i],seqi[i+1] = seqi[i+1],seqi[i]
						i += 1
					end
				end
			end
		end

		println(2, " length(lis)+length(ljs) > L : ", length(lis)+length(ljs) > L)

		if length(lis)+length(ljs) == L
			println(ceil(Integer, M))
			opts = [ find(x->x==0, M[i,:]) for i in 1:L ]
			cnt = ones(Int64, L)
			for k in 1:prod(map(length, opts))
				mapping = map(getindex, opts, cnt)
				good = true
				for i in 1:L, j in i+1:L; mapping[i] == mapping[j] && (good = false; break) end
				good && return mapping[1:L-addcols]
				for i in 1:L
					cnt[i] = cnt[i] % length(opts[i]) + 1
					cnt[i] ≠ 1 && break
				end
			end
		end

		println(3)

		#prepare next round
		minunc = minimum(M[setdiff(1:L, lis), setdiff(1:L, ljs)])
		M[lis,:] += minunc
		M[:,ljs] += minunc
		M -= minunc

		println(4)
	end
end
hungarian(A,B) = hungarian([ norm(a-b) for a in A, b in B ])

hungarian(rand(5,5))
hungarian(rand(4,5))
hungarian(rand(5,4))

hungarian([1,4,4], [1,4])
=#


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"The `Future` module implements future behavior of already existing functions,
which will replace the current version in a future release of Julia."
module Future

## copy!

"""
    Future.copy!(dst, src) -> dst

Copy `src` into `dst`.
For collections of the same type, copy the elements of `src` into `dst`,
discarding any pre-existing elements in `dst`.
Usually, `dst == src` holds after the call.
"""
copy!(dst::AbstractSet, src::AbstractSet) = union!(empty!(dst), src)
copy!(dst::AbstractDict, src::AbstractDict) = merge!(empty!(dst), src)
copy!(dst::AbstractVector, src::AbstractVector) = append!(empty!(dst), src)

function copy!(dst::AbstractArray, src::AbstractArray)
    size(dst) == size(src) || throw(ArgumentError(
        "arrays must have the same size for copy! (consider using `copyto!`)"))
    copyto!(dst, src)
end

end # module Future

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
struct CartesianSpace{N}
    dims::NTuple{N, Int}
    offsets::NTuple{N, Int}
end

Base.length(s::CartesianSpace) = prod(s.dims)
Base.eltype(::CartesianSpace{N}) where N = NTuple{N, Int}

function Base.iterate(s::CartesianSpace, st=1)
    st > length(s) && return
    return s[st], st + 1
end

function Base.getindex(s::CartesianSpace{1}, idx::Int)
    return idx + s.offsets[1]
end

function Base.getindex(s::CartesianSpace{2}, idx::Int)
    i2, i1 = divrem(idx - 1, s.dims[1])
    return i1 + s.offsets[1], i2 + s.offsets[2]
end

function Base.getindex(s::CartesianSpace{3}, idx::Int)
    i2, i1 = divrem(idx - 1, s.dims[1])
    i3, i2 = divrem(i2, s.dims[2])
    return i1 + s.offsets[1], i2 + s.offsets[2], i3 + s.offsets[3]
end

function Base.getindex(s::CartesianSpace{4}, idx::Int)
    i2, i1 = divrem(idx - 1, s.dims[1])
    i3, i2 = divrem(i2, s.dims[2])
    i4, i3 = divrem(i3, s.dims[3])
    return i1 + s.offsets[1], i2 + s.offsets[2], i3 + s.offsets[3], i4 + s.offsets[4]
end

function CartesianSpace(xs::AbstractUnitRange...)
    CartesianSpace(length.(xs), first.(xs))
end

struct CartesianPartition{N}
    space::CartesianSpace{N}
    from::Int
    to::Int
end

function schedule_task(s::CartesianSpace, tid::Int, nthreads::Int)
    len, rem = divrem(length(s), nthreads)
    f = (tid - 1) * len + 1
    l = f + len - 1

    if rem > 0
        if tid <= rem
            f = f + (tid - 1)
            l = l + tid
        else
            f = f + rem
            l = l + rem
        end
    end
    return CartesianPartition(s, f, l)
end

# NOTE:
# The following is a ~5% slower alternative to schedule_task.
# Not strictly equivalent, e.g., for length(s)=100, nthreads=7,
# schedule_task gives (15, 15, 14, 14, 14, 14, 14),
# but schedule_task_2 gives (14, 14, 14, 15, 14, 14, 15).
#
# function schedule_task_2(s::CartesianSpace, tid::Int, nthreads::Int)
#     return CartesianPartition(s,
#         length(s) * (tid-1) ÷ nthreads + 1,
#         length(s) * tid ÷ nthreads)
# end

Base.length(it::CartesianPartition) = it.to - it.from + 1
Base.eltype(it::CartesianPartition) = eltype(it.space)
Base.getindex(it::CartesianPartition, idx::Int) = it.space[idx]
function Base.iterate(it::CartesianPartition, st::Int=it.from)
    st > it.to && return
    return it.space[st], st + 1
end

CartesianSpace(::AbstractVector, subspace) = CartesianSpace(1:length(subspace))
CartesianSpace(S::AbstractMatrix, subspace) = CartesianSpace(1:size(S, 1), 1:length(subspace))
function CartesianSpace(S::AbstractMatrix{Complex{T}}, subspace) where {T <: Base.HWReal}
    return CartesianSpace(0:(size(S,1)-1) >>> 3, 1:length(subspace))
end


text:
function get_trace_bounds(ts::Int64, te::Int64, t0::Int64, t1::Int64, Δ::Int64, nx::Int64)
  i0 = 1
  i1 = nx
  while t0 < ts
    (i0 >= i1) && break
    t0 += Δ
    i0 += 1
  end

  while t1 > te
    (i1 <= i0) && break
    t1 -= Δ
    i1 -= 1
  end

  return i0, i1, t0
end

function get_trace_bound(t0::Int64, ts::Int64, Δ::Int64, nx::Int64)
  i0 = 1
  while ts < t0
    (i0 >= nx) && break
    ts  += Δ
    i0 += 1
  end

  return i0
end


text:
using ExtremeLearning

# Testing a simple binary adder

const nhidden = 100

x = [1.0 1.0; 0.0 1.0; 0.0 0.0; 1.0 0.0]
y = [0.0 1.0; 1.0 0.0; 0.0 0.0; 1.0 0.0]

elm = ELM(nhidden,x,y)
fit!(elm, x, y)

new_x = [1.0 1.0; 0.0 1.0; 0.0 0.0; 1.0 0.0]
y_pred = predict(elm, new_x)

@test y_pred[1, 1] < 0.1
@test y_pred[1, 2] > 0.9
@test y_pred[2, 1] > 0.9
@test y_pred[2, 2] < 0.1
@test y_pred[3, 1] < 0.1
@test y_pred[3, 2] < 0.1
@test y_pred[4, 1] > 0.9
@test y_pred[4, 2] < 0.1


text:
function make_vargs(vargs)
    arr = Tuple[]
    for (n,v) in vargs
        push!(arr, (Symbol(n),v))
    end
    arr
end

function isvalidcmd(cmd)
    isempty(cmd) && return false
    Base.is_id_start_char(cmd[1]) || return false
    for c in cmd
        Base.is_id_char(c) || return false
    end
    true
end

function parsepostdata(req, query)
    post_data = ""
    req_body = getfield(req, :body)

    if !isempty(req_body)
        idx = findfirst(isequal(UInt8('\0'), req_body))
        idx = (idx == nothing) ? lastindex(req_body) : (idx - 1)
        post_data = (idx == 0) ? String(req_body) : String(req_body[1:(idx-1)])
    end

    ('=' in post_data) || (return query)
    merge(query, HTTP.queryparams(post_data))
end

"""
Handles multipart form data.
Transfers all content as String.
Binary files can be uplaoded by encoding them with base64 first.
"""
const LFLF = UInt8['\n', '\n']
const CRLFCRLF = UInt8['\r', '\n', '\r', '\n']
function searcharr(haystack, needle, startpos)
    @static if VERSION < v"0.7.0-"
        search(haystack, needle, startpos)
    else
        Lh = length(haystack)
        Ln = length(needle)
        if Lh > Ln
            for idx in startpos:(Lh-Ln+1)
                r = idx:(idx+Ln-1)
                (view(haystack, r) == needle) && (return r)
            end
        end
        0:-1
    end
end
function parsepostdata(req, data_dict, multipart_boundary)
    data = getfield(req, :body)
    boundary_end = "--" * multipart_boundary * "--"
    boundary = "--" * multipart_boundary * "\r\n"
    Lbound = length(boundary)
    Ldata = length(data)

    boundbytes = convert(Vector{UInt8}, codeunits(boundary))
    boundbytes_end = convert(Vector{UInt8}, codeunits(boundary_end))

    boundloc = searcharr(data, boundbytes, 1)
    endpos = startpos = last(boundloc) + 1
    parts = Vector{Vector{UInt8}}()
    isend = false

    while !isend && ((startpos + Lbound) < Ldata)
        boundloc = searcharr(data, boundbytes, startpos)
        if isempty(boundloc)
            boundloc = searcharr(data, boundbytes_end, startpos)
            isend = true
        end
        if !isempty(boundloc)
            endpos = first(boundloc) - 3 # skip \r\n too

            part = data[startpos:endpos]
            push!(parts, part)

            startpos = last(boundloc) + 1
        end
    end

    for part in parts
        hdrloc1 = searcharr(part, LFLF, 1)
        hdrloc2 = searcharr(part, CRLFCRLF, 1)

        if length(hdrloc1) == 0
            hdrloc = hdrloc2
        elseif length(hdrloc2) == 0
            hdrloc = hdrloc1
        else
            hdrloc = (first(hdrloc1) < first(hdrloc2)) ? hdrloc1 : hdrloc2
        end

        parthdr = String(part[1:(first(hdrloc)-1)])
        partdata = part[(last(hdrloc)+1):end]

        collect_part_data(data_dict, parthdr, partdata)
    end
    data_dict
end

function collect_part_data(data_dict, parthdr, partdata)
    hdrdict = parse_part_headers(parthdr)

    # convert all to base64 for now
    #content_type = header(hdrdict, "Content-Type", "application/octet-stream")
    #istext = startswith(content_type, "text")

    content_disposition = header(hdrdict, "Content-Disposition")
    (isempty(content_disposition) || !startswith(content_disposition, "form-data;")) && return
    for attr in split(content_disposition, ';')
        ('=' in attr) || continue
        attr = strip(attr)
        (n,v) = split(attr, '='; limit=2)
        if n == "name"
            startswith(v, '"') && endswith(v, '"') && (v = v[2:(end-1)])
            data_dict[v] = base64encode(partdata)
        end
    end
end

function parse_part_headers(parthdr)
    hdrdict = Dict{String,String}()
    for hdr in split(parthdr, "\n")
        (':' in hdr) || continue
        (n,v) = split(hdr, ':'; limit=2)
        hdrdict[String(strip(n))] = String(strip(v))
    end
    hdrdict
end

header(req::HTTP.Request, name, default="") = HTTP.header(req, name, default)
function header(hdrdict::Dict, name, default="")
    for (n,v) in hdrdict
        (lowercase(n) == lowercase(name)) && (return v)
    end
    default
end

function get_multipart_form_boundary(req::HTTP.Request)
    content_type = header(req, "Content-Type")
    isempty(content_type) && (return nothing)
    parts = split(content_type, ";")
    (length(parts) < 2) && (return nothing)
    (lowercase(strip(parts[1])) == "multipart/form-data") || (return nothing)
    parts = split(strip(parts[2]), "=")
    (length(parts) < 2) && (return nothing)
    (lowercase(strip(parts[1])) == "boundary") || (return nothing)
    parts[2]
end

function http_handler(apis::Channel{APIInvoker{T,F}}, preproc::Function, req::HTTP.Request) where {T,F}
    @info("processing", target=getfield(req, :target))
    res = HTTP.Response(500)

    try
        comps = split(getfield(req, :target), '?', limit=2, keepempty=false)
        if isempty(comps)
            res = HTTP.Response(404)
        else
            res = preproc(req)
            if res === nothing
                comps = split(getfield(req, :target), '?', limit=2, keepempty=false)
                path = popfirst!(comps)
                data_dict = isempty(comps) ? Dict{String,String}() : HTTP.queryparams(comps[1])
                multipart_boundary = get_multipart_form_boundary(req)
                if multipart_boundary === nothing
                    data_dict = parsepostdata(req, data_dict)
                else
                    data_dict = parsepostdata(req, data_dict, multipart_boundary)
                end
                args = map(String, split(path, '/', keepempty=false))

                if isempty(args) || !isvalidcmd(args[1])
                    res = HTTP.Response(404)
                else
                    cmd = popfirst!(args)
                    @info("waiting for a handler")
                    api = take!(apis)
                    try
                        if isempty(data_dict)
                            @debug("calling", cmd, args)
                            res = httpresponse(api.format, apicall(api, cmd, args...))
                        else
                            vargs = make_vargs(data_dict)
                            @debug("calling", cmd, args, vargs)
                            res = httpresponse(api.format, apicall(api, cmd, args...; vargs...))
                        end
                    finally
                        put!(apis, api)
                    end
                end
            end
        end
    catch e
        res = HTTP.Response(500)
        Base.showerror(stderr, e, catch_backtrace())
        @error("Exception in handler: ", e)
    end
    @debug("response", res)
    return res
end

default_preproc(req::HTTP.Request) = nothing

# add a multipart form handler, provide default
struct HttpRpcServer{T,F}
    api::Channel{APIInvoker{T,F}}
    handler::HTTP.RequestHandlerFunction
end

HttpRpcServer(api::APIInvoker{T,F}, preproc::Function=default_preproc) where {T,F} = HttpRpcServer([api], preproc)
function HttpRpcServer(apis::Vector{APIInvoker{T,F}}, preproc::Function=default_preproc) where {T,F}
    api = Channel{APIInvoker{T,F}}(length(apis))
    for member in apis
        put!(api, member)
    end

    handler_fn = (req)->JuliaWebAPI.http_handler(api, preproc, req)
    handler = HTTP.RequestHandlerFunction(handler_fn)
    HttpRpcServer{T,F}(api, handler)
end

run_http(api::Union{Vector{APIInvoker{T,F}},APIInvoker{T,F}}, port::Int, preproc::Function=default_preproc; kwargs...) where {T,F} = run_http(HttpRpcServer(api, preproc), port; kwargs...)
function run_http(httprpc::HttpRpcServer{T,F}, port::Int; kwargs...) where {T,F}
    @info("running HTTP RPC server...")
    HTTP.listen(ip"0.0.0.0", port; kwargs...) do req
        HTTP.handle(httprpc.handler, req)
    end
end


text:
exp_name = "Competitors"

data_dirs = [
             # "ALOI", # too many observations
             "Annthyroid",
             # "Arrhythmia", # too many attributes
             "Cardiotocography",
             "Glass",
             "HeartDisease",
             "Hepatitis",
             "Ionosphere",
             # "InternetAds", # too many attributes
             # "KDDCup99", # too many observations
             "Lymphography",
             "PageBlocks",
             "Parkinson",
             "PenDigits",
             "Pima",
             "Shuttle",
             "SpamBase",
             "Stamps",
             "Waveform",
             "WBC",
             "WDBC",
             "WPBC"]

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
oracle_param = Dict{Symbol, Any}( :type => PoolOracle, :param => Dict{Symbol, Any}())

quality_metrics = Dict(:kappa => cohens_kappa, :f1 => f1_score)
selection_metric = :kappa

n_iter=50

learning_methods = vcat(
    # DFN author recommendation
    [DFNHeuristic(KernelLearning.DFNGammaStrategy(),
                  FixedCStrategy(0.05),
                  n_iter)],
    # DFN with random grid
    [DFNHeuristic(KernelLearning.DFNGammaStrategy(),
                  KernelLearning.RandomSampleGrid(solver,
                                                 range(0.001, stop=0.1, length=20),
                                                 nothing,
                                                 quality_metrics[selection_metric]),
                 n_iter)],
    # Wang parametrization from paper
    [WangHeuristic(init_strategy= KernelLearning.WangCombinedStrategy(solver,
                                                                     [10^x for x in -4.0:1.0:4.0],
                                                                     [0.01, 0.05],
                                                                     SVDD.f1_scoring)),
     # Wang extended parameter range
     WangHeuristic(init_strategy= KernelLearning.WangCombinedStrategy(solver,
                                                                      10.0.^range(-.1, stop=2, length=20),
                                                                      range(0.001, stop=0.1, length=20),
                                                                      SVDD.f1_scoring))],
    # Fully supervised baseline
    [SupervisedGridSearch(gamma_range = 10.0.^range(-.1, stop=2, length=20),
                          C_range = range(0.001, stop=0.1, length=20))],
    # QMS
    [QMSHeuristic(p=0.05)]
)

seeds = [42 * i for i in 1:5]

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
using Revise, ProjectedGPE

ħ = 1.0545718e-34
kB = 1.38064852e-23
amu = 1.660339040e-27
Bohr = 5.29e-11
#Rb87 mass and scattering length
m = 86.909180527*amu
as = 100*Bohr
#trap frequencies
ωx = 2π
ωy = 4ωx
ωz = 0.
#choice of time, length, energy units
t0 = 1.0/ωy
x0 = sqrt(ħ*t0/m)
E0 = ħ/t0
#interactions
#g  = (4*pi*ħ^2*as/m)*x0^3/E0 #dimensionless 3D
g = 0.1 #test 2D
#damping parameters (dimensionless)
γ = 0.05
ℳ  = 0.0
#chemical potential (dimensionless)
μ  = 12.0
#time evolution parameters
ti = 0.0
tf = 1.0/γ  #evolve for 2 damping times
Nt = 50
t = collect(linspace(ti,tf,Nt))
dt = 0.01π/μ #integrate step size [ - sh

basis = "Hermite"
ecut = 30*ħ*ωy/E0
Ω = [ωx; ωy]*t0
cinfo = makecinfo(ecut,Ω,basis)
tinfo = maketinfo(cinfo,4)
@unpack en,P,M = cinfo ;Mx,My = M
 c0   = randn(Mx,My)+im*randn(Mx,My); c0=P.*c0
ψ = tinfo.Tx*c0*tinfo.Ty'

function nlin!(dc,c,tinfo)
    #ψ = Tx*c*Ty'
    #dc.= Tx'*(W.*abs2.(ψ).*ψ)*Ty
    c2x!(ψ,c,tinfo)
    x2c!(dc,abs2.(ψ).*ψ,tinfo)
end

function Lgp!(dc,c,p,t)
    nlin!(dc,c,tinfo)
    dc .= P.*(-im*(1-im*γ)*((en - μ).*c .+ g*dc))
end
tinfo

dc=c0
nlin!(dc,c0,tinfo)
c2x!(ψ,c0,tinfo)
x2c!(dc,abs2.(ψ).*ψ,tinfo)


text:
requiredTypes = [ ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

# _______________________________________________________________type

export AbstractLevel
abstract type AbstractLevel end

# _______________________________________________________________functions

function setName(lvl::AbstractLevel, name::String)
    lvl.Name = name
end

function setMinStay(lvl::AbstractLevel, min::Int)
    lvl.MinStay = min
end

function setMaxStay(lvl::AbstractLevel, max::Int)
    lvl.MaxStay = max
end

function getMinStay(lvl::AbstractLevel)
    return lvl.MinStay
end

function getMaxStay(lvl::AbstractLevel)
    return lvl.MaxStay
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
# ## Literate part
#
# Head to the [Literate.jl](https://github.com/fredrikekre/Literate.jl) repo for more information
# about the syntax but, overall, it's about writing a Julia script with lots of comments.
#
# You can also convert a Jupyter notebook into a `.jl` script using Literate and include that with
# Franklin to get a webpage that reproduces the notebook.
#
# Julia code blocks will be executed and their output shown (unless you add a `;` at the end of
# the last line of the code block.
#
# Let's start with something very simple

x = 5

# And we can keep going with simple things

foreach(println, ('*'^i for i in 1:x))

# You can use packages of course, just make sure your website folder has a `Project.toml`
# with all your dependencies listed; also if you use the GitHub deployment path, you will need to make sure the deploy script installs the right package (head to the [docs](https://franklinjl.org) for more information).

using PyPlot

figure(figsize=(8, 6))
x = range(0, 10, length=250)
y = @. sin(3x) * exp(-x/5) + log(x) * cos(x) / 10
plot(x, y)

savefig(joinpath(@OUTPUT, "example1.svg")) # hide

# in the comments the markdown is treated as Franklin markdown so you can do things like
# include images inside div blocks etc:
#
# @@small-img \figalt{Simple example}{example1.svg} @@
#
# And include maths of course
#
# $$ \exp(i\pi) + 1 = 0 $$

max_stars_repo_name: Pramodh-G/ft-academic-resume

max_issues_repo_name: Pramodh-G/ft-academic-resume

max_forks_repo_name: Pramodh-G/ft-academic-resume