text:
using LoweredCodeUtils
using InteractiveUtils
using JuliaInterpreter
using JuliaInterpreter: finish_and_return!
using Core: CodeInfo
using Base.Meta: isexpr
using Test

module Lowering
using Parameters
struct Caller end
struct Gen{T} end
end

# Stuff for https://github.com/timholy/Revise.jl/issues/422
module Lowering422
const LVec{N, T} = NTuple{N, Base.VecElement{T}}
const LT{T} = Union{LVec{<:Any, T}, T}
const FloatingTypes = Union{Float32, Float64}
end

# Stuff for https://github.com/timholy/Revise.jl/issues/550
if Base.VERSION >= v"1.1"
    try
        using CBinding
    catch
        @info "Adding CBinding to the environment for test purposes"
        using Pkg
        Pkg.add("CBinding")    # not available for Julia 1.0
    end
    eval(:(module Lowering550
    using CBinding
    end))
end

bodymethtest0(x) = 0
function bodymethtest0(x)
    y = 2x
    y + x
end
bodymethtest1(x, y=1; z="hello") = 1
bodymethtest2(x, y=Dict(1=>2); z="hello") = 2
bodymethtest3(x::T, y=Dict(1=>2); z="hello") where T<:AbstractFloat = 3
# No kw but has optional args
bodymethtest4(x, y=1) = 4
bodymethtest5(x, y=Dict(1=>2)) = 5

@testset "Signatures" begin
    signatures = Set{Any}()
    newcode = CodeInfo[]
    for ex in (:(f(x::Int8; y=0) = y),
               :(f(x::Int16; y::Int=0) = 2),
               :(f(x::Int32; y="hello", z::Int=0) = 3),
               :(f(x::Int64;) = 4),
               :(f(x::Array{Float64,K}; y::Int=0) where K = K),
               # Keyword-arg functions that have an anonymous function inside
               :(fanon(list; sorted::Bool=true,) = sorted ? sort!(list, by=x->abs(x)) : list),
               # Keyword & default positional args
               :(g(x, y="hello"; z::Int=0) = 1),
               # Return type annotations
               :(annot(x, y; z::Bool=false,)::Nothing = nothing),
               # Generated methods
               quote
                   @generated function h(x)
                       if x <: Integer
                           return :(x ^ 2)
                       else
                           return :(x)
                       end
                   end
               end,
               quote
                   function h(x::Int, y)
                       if @generated
                           return y <: Integer ? :(x*y) : :(x+y)
                       else
                           return 2x+3y
                       end
                   end
               end,
               :(@generated genkw(; b=2) = nothing),  # https://github.com/timholy/Revise.jl/issues/290
               # Generated constructors
               quote
                   function Gen{T}(x) where T
                       if @generated
                           return T <: Integer ? :(x^2) : :(2x)
                       else
                           return 7x
                       end
                   end
               end,
               # Conditional methods
               quote
                   if 0.8 > 0.2
                       fctrue(x) = 1
                   else
                       fcfalse(x) = 1
                   end
               end,
               # Call methods
               :((::Caller)(x::String) = length(x)),
               )
        Core.eval(Lowering, ex)
        frame = Frame(Lowering, ex)
        rename_framemethods!(frame)
        pc = methoddefs!(signatures, frame; define=false)
        push!(newcode, frame.framecode.src)
    end

    # Manually add the signature for the Caller constructor, since that was defined
    # outside of manual lowering
    push!(signatures, Tuple{Type{Lowering.Caller}})

    nms = names(Lowering; all=true)
    modeval, modinclude = getfield(Lowering, :eval), getfield(Lowering, :include)
    failed = []
    for fsym in nms
        f = getfield(Lowering, fsym)
        isa(f, Base.Callable) || continue
        (f === modeval || f === modinclude) && continue
        for m in methods(f)
            if m.sig ∉ signatures
                push!(failed, m.sig)
            end
        end
    end
    @test isempty(failed)
    # Ensure that all names are properly resolved
    for code in newcode
        Core.eval(Lowering, code)
    end
    nms2 = names(Lowering; all=true)
    @test nms2 == nms
    @test Lowering.f(Int8(0)) == 0
    @test Lowering.f(Int8(0); y="LCU") == "LCU"
    @test Lowering.f(Int16(0)) == Lowering.f(Int16(0), y=7) == 2
    @test Lowering.f(Int32(0)) == Lowering.f(Int32(0); y=22) == Lowering.f(Int32(0); y=:cat, z=5) == 3
    @test Lowering.f(Int64(0)) == 4
    @test Lowering.f(rand(3,3)) == Lowering.f(rand(3,3); y=5) == 2
    @test Lowering.fanon([1,3,-2]) == [1,-2,3]
    @test Lowering.g(0) == Lowering.g(0,"LCU") == Lowering.g(0; z=5) == Lowering.g(0,"LCU"; z=5) == 1
    @test Lowering.annot(0,0) === nothing
    @test Lowering.h(2) == 4
    @test Lowering.h(2.0) == 2.0
    @test Lowering.h(2, 3) == 6
    @test Lowering.h(2, 3.0) == 5.0
    @test Lowering.fctrue(0) == 1
    @test_throws UndefVarError Lowering.fcfalse(0)
    @test (Lowering.Caller())("Hello, world") == 12
    g = Lowering.Gen{Float64}
    @test g(3) == 6

    # Don't be deceived by inner methods
    signatures = []
    ex = quote
        function fouter(x)
            finner(::Float16) = 2x
            return finner(Float16(1))
        end
    end
    Core.eval(Lowering, ex)
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    methoddefs!(signatures, frame; define=false)
    @test length(signatures) == 1
    @test LoweredCodeUtils.whichtt(signatures[1]) == first(methods(Lowering.fouter))

    # Check output of methoddef!
    frame = Frame(Lowering, :(function nomethod end))
    ret = methoddef!(empty!(signatures), frame; define=true)
    @test isempty(signatures)
    @test ret === nothing
    frame = Frame(Lowering, :(function amethod() nothing end))
    ret = methoddef!(empty!(signatures), frame; define=true)
    @test !isempty(signatures)
    @test isa(ret, NTuple{2,Int})

    # Anonymous functions in method signatures
    ex = :(max_values(T::Union{map(X -> Type{X}, Base.BitIntegerSmall_types)...}) = 1 << (8*sizeof(T)))  # base/abstractset.jl
    frame = Frame(Base, ex)
    rename_framemethods!(frame)
    signatures = Set{Any}()
    methoddef!(signatures, frame; define=false)
    @test length(signatures) == 1
    @test first(signatures) == which(Base.max_values, Tuple{Type{Int16}}).sig

    # define
    ex = :(fdefine(x) = 1)
    frame = Frame(Lowering, ex)
    empty!(signatures)
    methoddefs!(signatures, frame; define=false)
    @test_throws MethodError Lowering.fdefine(0)
    frame = Frame(Lowering, ex)
    empty!(signatures)
    methoddefs!(signatures, frame; define=true)
    @test Lowering.fdefine(0) == 1

    # define with correct_name!
    ex = quote
        @generated function generated1(A::AbstractArray{T,N}, val) where {T,N}
            ex = Expr(:tuple)
            for i = 1:N
                push!(ex.args, :val)
            end
            return ex
        end
    end
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=true)
    @test length(signatures) == 2
    @test Lowering.generated1(rand(2,2), 3.2) == (3.2, 3.2)
    ex = quote
        another_kwdef(x, y=1; z="hello") = 333
    end
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=true)
    @test length(signatures) == 5
    @test Lowering.another_kwdef(0) == 333
    ex = :(@generated genkw2(; b=2) = nothing)  # https://github.com/timholy/Revise.jl/issues/290
    frame = Frame(Lowering, ex)
    # rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=true)
    @test length(signatures) == 4
    @test Lowering.genkw2() === nothing

    # Test for correct exit (example from base/namedtuples.jl)
    ex = quote
        function merge(a::NamedTuple{an}, b::NamedTuple{bn}) where {an, bn}
            if @generated
                names = merge_names(an, bn)
                types = merge_types(names, a, b)
                vals = Any[ :(getfield($(sym_in(n, bn) ? :b : :a), $(QuoteNode(n)))) for n in names ]
                :( NamedTuple{$names,$types}(($(vals...),)) )
            else
                names = merge_names(an, bn)
                types = merge_types(names, typeof(a), typeof(b))
                NamedTuple{names,types}(map(n->getfield(sym_in(n, bn) ? b : a, n), names))
            end
        end
    end
    frame = Frame(Base, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    stmt = JuliaInterpreter.pc_expr(frame)
    if !LoweredCodeUtils.ismethod(stmt)
        pc = JuliaInterpreter.next_until!(LoweredCodeUtils.ismethod, frame, true)
    end
    pc, pc3 = methoddef!(signatures, frame; define=false)  # this tests that the return isn't `nothing`
    pc, pc3 = methoddef!(signatures, frame; define=false)
    @test length(signatures) == 2  # both the GeneratedFunctionStub and the main method

    # With anonymous functions in signatures
    ex = :(const BitIntegerType = Union{map(T->Type{T}, Base.BitInteger_types)...})
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=false)
    @test !isempty(signatures)

    for m in methods(bodymethtest0)
        @test bodymethod(m) === m
    end
    @test startswith(String(bodymethod(first(methods(bodymethtest1))).name), "#")
    @test startswith(String(bodymethod(first(methods(bodymethtest2))).name), "#")
    @test startswith(String(bodymethod(first(methods(bodymethtest3))).name), "#")
    @test bodymethod(first(methods(bodymethtest4))).nargs == 3  # one extra for #self#
    @test bodymethod(first(methods(bodymethtest5))).nargs == 3
    m = @which sum([1]; dims=1)
    # Issue in https://github.com/timholy/CodeTracking.jl/pull/48
    mbody = bodymethod(m)
    @test mbody != m && mbody.file != :none
    # varargs keyword methods
    m = which(Base.print_shell_escaped, (IO, AbstractString))
    mbody = bodymethod(m)
    @test isa(mbody, Method) && mbody != m

    ex = quote
        function print_shell_escaped(io::IO, cmd::AbstractString, args::AbstractString...;
                                     special::AbstractString="")
            print_shell_word(io, cmd, special)
            for arg in args
                print(io, ' ')
                print_shell_word(io, arg, special)
            end
        end
    end
    frame = Frame(Base, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=false)
    @test length(signatures) >= 3

    ex = :(typedsig(x) = 1)
    frame = Frame(Lowering, ex)
    methoddefs!(signatures, frame; define=true)
    ex = :(typedsig(x::Int) = 2)
    frame = Frame(Lowering, ex)
    JuliaInterpreter.next_until!(LoweredCodeUtils.ismethod3, frame, true)
    empty!(signatures)
    methoddefs!(signatures, frame; define=true)
    @test first(signatures).parameters[end] == Int

    # Multiple keyword arg methods per frame
    # (Revise issue #363)
    ex = quote
        keywrd1(x; kwarg=false) = 1
        keywrd2(x; kwarg="hello") = 2
        keywrd3(x; kwarg=:stuff) = 3
    end
    Core.eval(Lowering, ex)
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    pc, pc3 = methoddef!(signatures, frame; define=false)
    @test pc < length(frame.framecode.src.code)
    kw2sig = Tuple{typeof(Lowering.keywrd2), Any}
    @test kw2sig ∉ signatures
    pc = methoddefs!(signatures, frame; define=false)
    @test pc === nothing
    @test kw2sig ∈ signatures

    # Module-scoping
    ex = :(Base.@irrational π        3.14159265358979323846  pi)
    frame = Frame(Base.MathConstants, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=false)
    @test !isempty(signatures)

    # Inner methods in structs. Comes up in, e.g., Core.Compiler.Params.
    # The body of CustomMS is an SSAValue.
    ex = quote
        struct MyStructWithMeth
            x::Int
            global function CustomMS(;x=1)
                return new(x)
            end
            MyStructWithMeth(x) = new(x)
        end
    end
    Core.eval(Lowering, ex)
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
    empty!(signatures)
    methoddefs!(signatures, frame; define=false)
    @test Tuple{typeof(Lowering.CustomMS)} ∈ signatures

    # https://github.com/timholy/Revise.jl/issues/398
    ex = quote
        @with_kw struct Items
            n::Int
            items::Vector{Int} = [i for i=1:n]
        end
    end
    Core.eval(Lowering, ex)
    frame = Frame(Lowering, ex)
    dct = rename_framemethods!(frame)
    ks = collect(filter(k->startswith(String(k), "#Items#"), keys(dct)))
    @test length(ks) == 2
    @test dct[ks[1]] == dct[ks[2]]
    @test isdefined(Lowering, ks[1]) || isdefined(Lowering, ks[2])
    nms = filter(sym->occursin(r"#Items#\d+#\d+", String(sym)), names(Lowering; all=true))
    @test length(nms) == 1

    # https://github.com/timholy/Revise.jl/issues/422
    ex = :(@generated function fneg(x::T) where T<:LT{<:FloatingTypes}
         s = """
         %2 = fneg $(llvm_type(T)) %0
         ret $(llvm_type(T)) %2
         """
         return :(
             $(Expr(:meta, :inline));
             Base.llvmcall($s, T, Tuple{T}, x)
         )
     end)
    empty!(signatures)
    Core.eval(Lowering422, ex)
    frame = Frame(Lowering422, ex)
    rename_framemethods!(frame)
    pc = methoddefs!(signatures, frame; define=false)
    @test typeof(Lowering422.fneg) ∈ Set(Base.unwrap_unionall(sig).parameters[1] for sig in signatures)

    # https://github.com/timholy/Revise.jl/issues/550
    if Base.VERSION >= v"1.1"
        ex = :(@cstruct S {
            val::Int8
        })
        empty!(signatures)
        Core.eval(Lowering550, ex)
        frame = Frame(Lowering550, ex)
        rename_framemethods!(frame)
        pc = methoddefs!(signatures, frame; define=false)
        @test !isempty(signatures)   # really we just need to know that `methoddefs!` completed without getting stuck
    end

    # Scoped names (https://github.com/timholy/Revise.jl/issues/568)
    ex = :(f568() = -1)
    Core.eval(Lowering, ex)
    @test Lowering.f568() == -1
    empty!(signatures)
    ex = :(f568() = -2)
    frame = Frame(Lowering, ex)
    pcstop = findfirst(LoweredCodeUtils.ismethod3, frame.framecode.src.code)
    pc = 1
    while pc < pcstop
        pc = JuliaInterpreter.step_expr!(finish_and_return!, frame, true)
    end
    pc = methoddef!(finish_and_return!, signatures, frame, pc; define=true)
    @test Tuple{typeof(Lowering.f568)} ∈ signatures
    @test Lowering.f568() == -2

    # Undefined names
    # This comes from FileWatching; WindowsRawSocket is only defined on Windows
    ex = quote
        if Sys.iswindows()
            using Base: WindowsRawSocket
            function wait(socket::WindowsRawSocket; readable=false, writable=false)
                fdw = _FDWatcher(socket, readable, writable)
                try
                    return wait(fdw, readable=readable, writable=writable)
                finally
                    close(fdw, readable, writable)
                end
            end
        end
    end
    frame = Frame(Lowering, ex)
    rename_framemethods!(frame)
end


text:
struct PWLFunction{D}
    x::Vector{NTuple{D,Float64}}
    z::Vector{Float64}
    T::Vector{Vector{Int}}
    meta::Dict
end

function PWLFunction{D}(x::Vector{NTuple{D}}, z::Vector, T::Vector{Vector}, meta::Dict) where {D}
    @assert length(x) == length(z)
    for t in T
        @assert minimum(t) > 0 && maximum(t) <= length(x)
    end
    return PWLFunction{D}(x, z, T, meta)
end

PWLFunction(x, z, T) = PWLFunction(x, z, T, Dict())

const UnivariatePWLFunction = PWLFunction{1}

function UnivariatePWLFunction(x, z)
    @assert issorted(x)
    return PWLFunction(Tuple{Float64}[(xx,) for xx in x], convert(Vector{Float64}, z), [[i,i+1] for i in 1:length(x)-1])
end

function UnivariatePWLFunction(x, fz::Function)
    @assert issorted(x)
    return PWLFunction(Tuple{Float64}[(xx,) for xx in x], map(t->convert(Float64,fz(t)), x), [[i,i+1] for i in 1:length(x)-1])
end

const BivariatePWLFunction = PWLFunction{2}

function BivariatePWLFunction(x, y, fz::Function; pattern=:BestFit, seed=hash((length(x),length(y))))
    @assert issorted(x)
    @assert issorted(y)
    X = vec(collect(Base.product(x,y)))
    # X = vec(Tuple{Float64,Float64}[(_x,_y) for _x in x, _y in y])
    Z = map(t -> convert(Float64,fz(t...)), X)
    T = Vector{Vector{Int}}()
    m = length(x)
    n = length(y)

    mt = MersenneTwister(seed)
    # run for each square on [x[i],x[i+1]] × [y[i],y[i+1]]
    for i in 1:length(x)-1, j in 1:length(y)-1
        SWt, NWt, NEt, SEt = LinearIndices((m,n))[i,j], LinearIndices((m,n))[i,j+1], LinearIndices((m,n))[i+1,j+1], LinearIndices((m,n))[i+1,j]
        xL, xU, yL, yU = x[i], x[i+1], y[j], y[j+1]
        @assert xL == X[SWt][1] == X[NWt][1]
        @assert xU == X[SEt][1] == X[NEt][1]
        @assert yL == X[SWt][2] == X[SEt][2]
        @assert yU == X[NWt][2] == X[NEt][2]
        SW, NW, NE, SE = Z[SWt], Z[NWt], Z[NEt], Z[SEt]
        mid1 = 0.5*(SW+NE)
        mid2 = 0.5*(NW+SE)

        if pattern == :Upper
            if mid1 > mid2
                t1 = [SWt,NWt,NEt]
                t2 = [SWt,NEt,SEt]
            else
                t1 = [SWt,NWt,SEt]
                t2 = [SEt,NWt,NEt]
            end
        elseif pattern == :Lower
            if mid1 > mid2
                t1 = [SWt,NWt,SEt]
                t2 = [SEt,NWt,NEt]
            else
                t1 = [SWt,NWt,NEt]
                t2 = [SWt,NEt,SEt]
            end
        elseif pattern == :BestFit
            mid3 = fz(0.5*(xL+xU), 0.5*(yL+yU))
            if abs(mid1-mid3) < abs(mid2-mid3)
                t1 = [SWt,NWt,NEt]
                t2 = [SWt,NEt,SEt]
            else
                t1 = [SWt,NWt,SEt]
                t2 = [SEt,NWt,NEt]
            end
        elseif pattern == :UnionJack
            t1 = [SWt,SEt]
            t2 = [NWt,NEt]
            if iseven(i+j)
                push!(t1, NWt)
                push!(t2, SEt)
            else
                push!(t1, NEt)
                push!(t2, SWt)
            end
        elseif pattern == :K1
            t1 = [SEt,SWt,NWt]
            t2 = [NWt,NEt,SEt]
        elseif pattern == :Random
            if rand(mt, Bool)
                t1 = [NWt,NEt,SEt]
                t2 = [SEt,SWt,NWt]
            else
                t1 = [SWt,NWt,NEt]
                t2 = [NEt,SEt,SWt]
            end
        else
            error("pattern $pattern not currently supported")
        end

        push!(T, t1)
        push!(T, t2)
    end

    return PWLFunction{2}(X, Z, T, Dict(:structure=>pattern))
end

ext: jl


text:
import Random: seed!

# use our own random seed for msg_id so that we
# don't alter the user-visible random state (issue #336)
const IJulia_RNG = seed!(Random.MersenneTwister(0))
import UUIDs
uuid4() = repr(UUIDs.uuid4(IJulia_RNG))

const orig_stdin  = Ref{IO}()
const orig_stdout = Ref{IO}()
const orig_stderr = Ref{IO}()
const SOFTSCOPE = Ref{Bool}()
function __init__()
    seed!(IJulia_RNG)
    orig_stdin[]  = stdin
    orig_stdout[] = stdout
    orig_stderr[] = stderr
    SOFTSCOPE[] = lowercase(get(ENV, "IJULIA_SOFTSCOPE", "yes")) in ("yes", "true")
end

# the following constants need to be initialized in init().
const publish = Ref{Socket}()
const raw_input = Ref{Socket}()
const requests = Ref{Socket}()
const control = Ref{Socket}()
const heartbeat = Ref{Socket}()
const profile = Dict{String,Any}()
const read_stdout = Ref{Base.PipeEndpoint}()
const read_stderr = Ref{Base.PipeEndpoint}()
const socket_locks = Dict{Socket,ReentrantLock}()

# similar to Pkg.REPLMode.MiniREPL, a minimal REPL-like emulator
# for use with Pkg.do_cmd.  We have to roll our own to
# make sure it uses the redirected stdout, and because
# we don't have terminal support.
import REPL
struct MiniREPL <: REPL.AbstractREPL
    display::TextDisplay
end
REPL.REPLDisplay(repl::MiniREPL) = repl.display
const minirepl = Ref{MiniREPL}()

function init(args)
    inited && error("IJulia is already running")
    if length(args) > 0
        merge!(profile, open(JSON.parse,args[1]))
        verbose && println("PROFILE = $profile")
        global connection_file = args[1]
    else
        # generate profile and save
        let port0 = 5678
            merge!(profile, Dict{String,Any}(
                "ip" => "127.0.0.1",
                "transport" => "tcp",
                "stdin_port" => port0,
                "control_port" => port0+1,
                "hb_port" => port0+2,
                "shell_port" => port0+3,
                "iopub_port" => port0+4,
                "key" => uuid4()
            ))
            fname = "profile-$(getpid()).json"
            global connection_file = "$(pwd())/$fname"
            println("connect ipython with --existing $connection_file")
            open(fname, "w") do f
                JSON.print(f, profile)
            end
        end
    end

    if !isempty(profile["key"])
        signature_scheme = get(profile, "signature_scheme", "hmac-sha256")
        isempty(signature_scheme) && (signature_scheme = "hmac-sha256")
        sigschm = split(signature_scheme, "-")
        if sigschm[1] != "hmac" || length(sigschm) != 2
            error("unrecognized signature_scheme $signature_scheme")
        end
        hmacstate[] = MbedTLS.MD(getfield(MbedTLS, Symbol("MD_", uppercase(sigschm[2]))),
                                 profile["key"])
    end

    publish[] = Socket(PUB)
    raw_input[] = Socket(ROUTER)
    requests[] = Socket(ROUTER)
    control[] = Socket(ROUTER)
    heartbeat[] = Socket(ROUTER)
    bind(publish[], "$(profile["transport"])://$(profile["ip"]):$(profile["iopub_port"])")
    bind(requests[], "$(profile["transport"])://$(profile["ip"]):$(profile["shell_port"])")
    bind(control[], "$(profile["transport"])://$(profile["ip"]):$(profile["control_port"])")
    bind(raw_input[], "$(profile["transport"])://$(profile["ip"]):$(profile["stdin_port"])")
    bind(heartbeat[], "$(profile["transport"])://$(profile["ip"]):$(profile["hb_port"])")

    # associate a lock with each socket so that multi-part messages
    # on a given socket don't get inter-mingled between tasks.
    for s in (publish[], raw_input[], requests[], control[], heartbeat[])
        socket_locks[s] = ReentrantLock()
    end

    start_heartbeat(heartbeat[])
    if capture_stdout
        read_stdout[], = redirect_stdout()
        redirect_stdout(IJuliaStdio(stdout,"stdout"))
    end
    if capture_stderr
        read_stderr[], = redirect_stderr()
        redirect_stderr(IJuliaStdio(stderr,"stderr"))
    end
    redirect_stdin(IJuliaStdio(stdin,"stdin"))
    minirepl[] = MiniREPL(TextDisplay(stdout))

    logger = Base.CoreLogging.SimpleLogger(Base.stderr)
    Base.CoreLogging.global_logger(logger)

    send_status("starting")
    global inited = true
end


text:
using DeepNets
using Test

@testset "DeepNets.jl" begin
    # Write your tests here.
end


text:
using Documenter
using ScanImageTiffReader

push!(LOAD_PATH, "../src")

makedocs(
    modules = [ScanImageTiffReader],
    doctest = true,
    format = :html,
    sitename = "ScanImage Tiff Reader",
    pages = [
        "index.md"
    ]
)


text:
using Test
using TDD

@testset "interval.jl" begin
    iv = TDD.Interval(1,32)
    @test minimum(iv) == 1
    @test maximum(iv) == 32
    @test (4 in iv) == true
    @test (40 in iv) == false
    @test (-1.5 in iv) == false
    @test isempty(iv) == false
    @test issubset([1,2,3], iv) == true
    @test issubset([1,2,3,40], iv) == false
    @test issubset([], iv) == true
    @test issubset(1:32, iv) == true
    @test intersect([1,2,3], iv) == [1,2,3]
    @test intersect([40], iv) == []
    @test intersect([], iv) == []

    io = IOBuffer()
    show(io, iv)
    str = String(take!(io))
    @test str == "\u301a1:32\u301b"

    iv = TDD.Interval()
    @test isempty(iv) == true
    @test minimum(iv) === nothing
    @test maximum(iv) === nothing
    @test (4 in iv) == false
    @test issubset([], iv) == true
    @test issubset(1:32, iv) == false
    @test intersect([1,2,3], iv) == []
    @test intersect([], iv) == []

    io = IOBuffer()
    show(io, iv)
    str = String(take!(io))
    @test str == "∅"

end


text:
module CIFAR10Reader

export

    readdata!,
    readdata

const NROW = 32
const NCOL = 32
const NCHAN = 3
const NBYTE = NROW * NCOL * NCHAN + 1 # "+ 1" for label
const CHUNK_SIZE = 10_000

function readnext!(buffer::Array{UInt8}, io::IO)
    y = Int(read(io, UInt8))
    read!(io, buffer)
    buffer, y
end

function readdata!(buffer::Array{UInt8}, io::IO, index::Integer)
    seek(io, (index - 1) * NBYTE)
    readnext!(buffer, io)
end

function readdata(io::IO, index::Integer)
    buffer = Array{UInt8}(undef, NROW, NCOL, NCHAN)
    readdata!(buffer, io, index)
end

function readdata(io::IO)
    X = Array{UInt8}(undef, NROW, NCOL, NCHAN, CHUNK_SIZE)
    Y = Array{Int}(undef, CHUNK_SIZE)
    buffer = Array{UInt8}(undef, NROW, NCOL, NCHAN)
    @inbounds for index in 1:CHUNK_SIZE
        _, ty = readnext!(buffer, io)
        copyto!(view(X,:,:,:,index), buffer)
        Y[index] = ty
    end
    X, Y
end

function readdata(file::AbstractString, index::Integer)
    open(file, "r") do io
        readdata(io, index)
    end::Tuple{Array{UInt8,3},Int}
end

function readdata(file::AbstractString)
    open(file, "r") do io
        readdata(io)
    end::Tuple{Array{UInt8,4},Vector{Int}}
end

end


text:
clima_formatter_options = (
    indent = 4,
    margin = 80,
    always_for_in = true,
    whitespace_typedefs = true,
    whitespace_ops_in_indices = true,
    remove_extra_newlines = false,
)


text:
using GenGlobal
using Test
using LinearAlgebra

using GenGlobal.testModule: set_globalx, get_globalx

@show dump(@GenGlobal somevar)

@testset "GenGlobal Tests" begin

    x = 1.0
    a = Matrix(1.0I,3,3)

    set_globalx(x)
    @test x === get_globalx()

    set_globalx(a)
    a === get_globalx()

    @test get_globalx() != x
    set_globalx(x)
    @test get_globalx() == x
    @test get_globalx() === x
    @test get_globalx() != a

end


text:
module TestChess

using Test
import Moves
import Util
import Magic

@testset "moves" verbose=true begin
    @test Moves.rookMoves(0,1) == parse(UInt64,
    "00000001"*
    "00000001"*
    "00000001"*
    "00000001"*
    "00000001"*
    "00000001"*
    "00000001"*
    "11111110",base=2)

    @test Moves.rookMoves(0,64) == parse(UInt64,
    "01111111"*
    "10000000"*
    "10000000"*
    "10000000"*
    "10000000"*
    "10000000"*
    "10000000"*
    "10000000",base=2)

    @test Moves.bishopMoves(0,1) == parse(UInt64,
    "10000000"*
    "01000000"*
    "00100000"*
    "00010000"*
    "00001000"*
    "00000100"*
    "00000010"*
    "00000000",base=2)

    @test Moves.bishopMoves(0,64) == parse(UInt64,
    "00000000"*
    "01000000"*
    "00100000"*
    "00010000"*
    "00001000"*
    "00000100"*
    "00000010"*
    "00000001",base=2)

    @test Moves.queenMoves(0,32) == parse(UInt64,
    "10001000"*
    "10010000"*
    "10100000"*
    "11000000"*
    "01111111"*
    "11000000"*
    "10100000"*
    "10010000",base=2)
end

@testset "check " verbose = true begin
    board = Util.startPositions()
    @test Moves.isCheck(board,'B') == false

    board = Util.emptyPositions()
    board["BP"] = Util.setBit(board["BP"], 11,1)
    board["WK"] = Util.setBit(board["WK"],4,1)
    @test Moves.isCheck(board,'W') == true

    board = Util.emptyPositions()
    board["WN"] = Util.setBit(board["WN"], 52,1)
    board["BK"] = Util.setBit(board["BK"],37,1)
    @test Moves.isCheck(board,'B') == true

    board = Util.emptyPositions()
    board["BQ"] = Util.setBit(board["BQ"], 14,1)
    board["WK"] = Util.setBit(board["WK"],8,1)
    @test Moves.isCheck(board,'B') == false

    board = Util.emptyPositions()
    board["BQ"] = Util.setBit(board["BQ"], 1,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    board["BR"] = Util.setBit(board["BR"], 9,1)
    @test Moves.isCheckMate(board,'W',[]) == 1

    #white in checkmate surrounded by white pieces
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],13,1)
    board["WP"] = Util.setBit(board["WP"],12,1)
    board["WB"] = Util.setBit(board["WB"],6,1)
    board["WQ"] = Util.setBit(board["WQ"],4,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    board["BQ"] = Util.setBit(board["BQ"],32,1)
    @test Moves.isCheckMate(board,'W',[]) == 1

    #white in check, king could take a piece, that piece protected
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],14,1)
    board["BR"] = Util.setBit(board["BR"],13,1)
    board["BB"] = Util.setBit(board["BB"],29,1)
    board["BQ"] = Util.setBit(board["BQ"],4,1)
    board["WK"] = Util.setBit(board["WK"],6,1)
    @test Moves.isCheckMate(board,'W',[]) == 1

    #white in check but piece can intervene
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],14,1)
    board["WR"] = Util.setBit(board["WR"],13,1)
    board["BB"] = Util.setBit(board["BB"],29,1)
    board["BQ"] = Util.setBit(board["BQ"],4,1)
    board["WK"] = Util.setBit(board["WK"],6,1)
    @test Moves.isCheckMate(board,'W',[]) == 0
    @test Moves.isCheck(board,'W') == true

    #stalemate
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],50,1)
    board["WQ"] = Util.setBit(board["WR"],42,1)
    board["BK"] = Util.setBit(board["BB"],58,1)
    @test Moves.isCheckMate(board,'B',[]) == 2
    @test Moves.isCheck(board,'B') == false

    #capture out of check
    board = Util.emptyPositions()
    board["BQ"] = Util.setBit(board["BQ"],60,1)
    board["WQ"] = Util.setBit(board["WQ"],59,1)
    board["BK"] = Util.setBit(board["BK"],61,1)
    @test Moves.isCheck(board,'B') == false

    #king capture out of checkmate
    board = Util.emptyPositions()
    board["BQ"] = Util.setBit(0,63,1)
    board["BR"] = Util.setBit(0,49,1)
    board["WK"] = Util.setBit(0,64,1)
    @test Moves.isCheck(board,'W') == true
    @test Moves.isCheckMate(board,'W',[]) == 0

    #cannot move into check

    board = Util.emptyPositions()
    board["BK"] = Util.setBit(0,60,1)
    board["WQ"] = Util.setBit(0,46,1)
    board["BP"] = Util.setBit(0,53,1)
    Util.prettyPrintBoard(board)
    @test length(Moves.getMoves(board,'B',[])) == 5

    board = Dict("BR" => 0x8100000000000000, "WK" => 0x0000000000000008, "WB" => 0x0000000000000004, "WQ" => 0x0000000004000000, "BQ" => 0x0000000000000020, "BB" => 0x2000100000000000, "BN" => 0x4200000000000000, "BK" => 0x1000000000000000, "WR" => 0x0000000000000001, "WN" => 0x0000000000000002, "BP" => 0x00e7001000000000, "WP" => 0x00000042110c8000)
    @test Moves.isCheckMate(board,'W',[]) == 0

end

@testset "pawn promotion" verbose = true begin

    #pawn promotion forward to queen
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],56,1)
    board = Moves.makeMoveUCI("h7h8Q",board,'W')
    @test board["WQ"] > 0
    @test board["WP"] == 0

    #pawn promotion forward to queen
    board = Util.emptyPositions()
    board["BP"] = Util.setBit(board["BP"],9,1)
    @test "a2a1Q" in Moves.getMoves(board,'B',[])
end

@testset "en passant" verbose = true begin

    #simple en passant of white pawn
    history = []
    board = Util.emptyPositions()
    board["WP"] = Util.setBit(board["WP"],37,1)
    board["BP"] = Util.setBit(board["BP"],38,1)
    push!(history,"e3e5")
    moves = Moves.getMoves(board,'B',history)
    @test "f5e4" in moves

    # #not en passant of black pawn
    # history = []
    # board = emptyPositions()
    # board["WP"] = Util.setBit(board["WP"],37,1)
    # board["BP"] = Util.setBit(board["BP"],38,1)
    # push!(history,"e3e5")
    # moves = getMovesUCI(board,'B',history)
    # @test "" not in moves
end

@testset "castling" verbose = true begin
    #simple en passant of white pawn
    history = []

    #black castling king side
    board = Util.emptyPositions()
    board["BR"] = Util.setBit(board["BR"],64,1)
    board["BK"] = Util.setBit(board["BK"],61,1)
    occ = reduce((x,y) -> x | y, values(board))
    @test "e8g8" in Moves.castleMoves(board,'B',history,occ)

    #white castle queenside
    board = Util.emptyPositions()
    board["WR"] = Util.setBit(board["WR"],1,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    occ = reduce((x,y) -> x | y, values(board))
    @test "e1c1" in Moves.castleMoves(board,'W',history,occ)

    #rook already moved
    board = Util.emptyPositions()
    board["WR"] = Util.setBit(board["WR"],1,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    occ = reduce((x,y) -> x | y, values(board))
    history = ["a3a1"]
    @test length(Moves.castleMoves(board,'W',history,occ)) == 0

    #king already moved
    board = Util.emptyPositions()
    board["WR"] = Util.setBit(board["WR"],1,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    occ = reduce((x,y) -> x | y, values(board))
    history = ["e2e1"]
    @test length(Moves.castleMoves(board,'W',history,occ)) == 0

    #piece in the way
    board = Util.emptyPositions()
    board["WR"] = Util.setBit(board["WR"],1,1)
    board["WK"] = Util.setBit(board["WK"],5,1)
    board["WB"] = Util.setBit(board["WB"],3,1)
    occ = reduce((x,y) -> x | y, values(board))
    history = []
    @test length(Moves.castleMoves(board,'W',history,occ)) == 0

    #castle through a check
    board = Util.emptyPositions()
    board["BR"] = Util.setBit(board["BR"],60,1)
    board["BK"] = Util.setBit(board["BK"],64,1)
    board["WQ"] = Util.setBit(board["WQ"],9,1)
    occ = reduce((x,y) -> x | y, values(board))
    history = []
    @test length(Moves.castleMoves(board,'B',history,occ)) == 0
end

@testset "search" verbose = true begin
end

@testset "eval" verbose = true begin

    board = Util.emptyPositions()
    board["WR"] = Util.setBit(0,16,1)
    board["BR"] = Util.setBit(0,64,1)
    board["BB"] = Util.setBit(0,9,1)
    @test count_ones(Moves.getAttackBoard(board,'W',false) & Util.getBlackPieces(board)) == 2

end

@testset "util" verbose = true begin
    #least sig bit
#    @test 18  == leastSignificantBit(1130822006472704)
#     @test leastSignificantBit(1130822006734848) == 10

end

end

max_stars_repo_name: grahamclyne/chessEngine

max_issues_repo_name: grahamclyne/chessEngine

max_forks_repo_name: grahamclyne/chessEngine


text:
using .Models

# Get the tensor network of <peps|peps'>
function inner_network(peps::PEPS, peps_prime::PEPS)
  return vcat(vcat(peps.data...), vcat(peps_prime.data...))
end

# Get the tensor network of <peps|mpo|peps'>
# The local MPO specifies the 2-site term of the Hamiltonian
function inner_network(
  peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, mpo::MPO, coordinates::Array
)
  @assert(length(mpo) == length(coordinates))
  network = vcat(peps.data...)
  dimy, dimx = size(peps.data)
  for ii in 1:dimx
    for jj in 1:dimy
      if (jj, ii) in coordinates
        index = findall(x -> x == (jj, ii), coordinates)
        @assert(length(index) == 1)
        network = vcat(network, [mpo.data[index[1]]])
        network = vcat(network, [peps_prime_ham.data[jj, ii]])
      else
        network = vcat(network, [peps_prime.data[jj, ii]])
      end
    end
  end
  return network
end

"""Generate an array of networks representing inner products, <p|H_1|p>, ..., <p|H_n|p>
Parameters
----------
peps: a peps network with datatype PEPS
peps_prime: prime of peps used for inner products
peps_prime_ham: prime of peps used for calculating expectation values
Hs: An array of MPO operators with datatype LocalMPO
Returns
-------
An array of networks.
"""
function inner_networks(peps::PEPS, peps_prime::PEPS, peps_prime_ham::PEPS, Hs::Array)
  network_list = Vector{Vector{ITensor}}()
  for H_term in Hs
    if H_term isa Models.LocalMPO
      coords = [H_term.coord1, H_term.coord2]
    elseif H_term isa Models.LineMPO
      if H_term.coord[1] isa Colon
        coords = [(i, H_term.coord[2]) for i in 1:length(H_term.mpo)]
      else
        coords = [(H_term.coord[1], i) for i in 1:length(H_term.mpo)]
      end
    end
    inner = inner_network(peps, peps_prime, peps_prime_ham, H_term.mpo, coords)
    network_list = vcat(network_list, [inner])
  end
  return network_list
end

function inner_network(peps::PEPS, peps_prime::PEPS, ::typeof(tree))
  Ny, Nx = size(peps.data)
  network = SubNetwork(peps.data[1, :]..., peps_prime.data[1, :]...)
  for i in 2:Ny
    network = SubNetwork(network, peps.data[i, :]..., peps_prime.data[i, :]...)
  end
  return network
end

function inner_network(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  mpo::MPO,
  coordinate::Tuple{<:Integer,Colon},
  ::typeof(tree),
)
  Ny, Nx = size(peps.data)
  network1, network2 = SubNetwork(), SubNetwork()
  if coordinate[1] > 1
    network1 = SubNetwork(peps.data[1, :]..., peps_prime.data[1, :]...)
    for i in 2:(coordinate[1] - 1)
      network1 = SubNetwork(network1, peps.data[i, :]..., peps_prime.data[i, :]...)
    end
  end
  if coordinate[1] < Ny
    network2 = SubNetwork(peps.data[Ny, :]..., peps_prime.data[Ny, :]...)
    for i in reverse((coordinate[1] + 1):(Ny - 1))
      network2 = SubNetwork(network2, peps.data[i, :]..., peps_prime.data[i, :]...)
    end
  end
  if coordinate[1] == 1
    return SubNetwork(
      network2, peps.data[1, :]..., peps_prime_ham.data[1, :]..., mpo.data...
    )
  elseif coordinate[1] == Ny
    return SubNetwork(
      network1, peps.data[Ny, :]..., peps_prime_ham.data[Ny, :]..., mpo.data...
    )
  else
    index = coordinate[1]
    return SubNetwork(
      network1,
      network2,
      peps.data[index, :]...,
      peps_prime_ham.data[index, :]...,
      mpo.data...,
    )
  end
end

function inner_network(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  mpo::MPO,
  coordinate::Tuple{Colon,<:Integer},
  ::typeof(tree),
)
  Ny, Nx = size(peps.data)
  network1, network2 = SubNetwork(), SubNetwork()
  if coordinate[2] > 1
    network1 = SubNetwork(peps.data[:, 1]..., peps_prime.data[:, 1]...)
    for i in 2:(coordinate[2] - 1)
      network1 = SubNetwork(network1, peps.data[:, i]..., peps_prime.data[:, i]...)
    end
  end
  if coordinate[2] < Nx
    network2 = SubNetwork(peps.data[:, Nx]..., peps_prime.data[:, Nx]...)
    for i in reverse((coordinate[2] + 1):(Nx - 1))
      network2 = SubNetwork(network2, peps.data[:, i]..., peps_prime.data[:, i]...)
    end
  end
  if coordinate[2] == 1
    return SubNetwork(
      network2, peps.data[:, 1]..., peps_prime_ham.data[:, 1]..., mpo.data...
    )
  elseif coordinate[2] == Nx
    return SubNetwork(
      network1, peps.data[:, Nx]..., peps_prime_ham.data[:, Nx]..., mpo.data...
    )
  else
    index = coordinate[2]
    return SubNetwork(
      network1,
      network2,
      peps.data[:, index]...,
      peps_prime_ham.data[:, index]...,
      mpo.data...,
    )
  end
end

function inner_networks(
  peps::PEPS,
  peps_prime::PEPS,
  peps_prime_ham::PEPS,
  Hs::Vector{Models.LineMPO},
  ::typeof(tree),
)
  function generate_each_network(H)
    return inner_network(peps, peps_prime, peps_prime_ham, H.mpo, H.coord, tree)
  end
  return [generate_each_network(H) for H in Hs]
end


text:
"""
    AbstractRipsFiltration{I<:Signed, T} <: AbstractFiltration{I, T}

An abstract Vietoris-Rips filtration. Its subtypes can only overload
[`adjacency_matrix`](@ref) and get default implementations for the rest of the filtration
interface.

# Example

```jldoctest
julia> struct MyRips <: Ripserer.AbstractRipsFiltration{Int, Float16} end

julia> Ripserer.adjacency_matrix(::MyRips) = [0 1 1; 1 0 1; 1 1 0]

julia> ripserer(MyRips())
2-element Vector{PersistenceDiagramsBase.PersistenceDiagram}:
 3-element 0-dimensional PersistenceDiagram
 0-element 1-dimensional PersistenceDiagram

```
"""
abstract type AbstractRipsFiltration{I<:Signed,T} <: AbstractFiltration{I,T} end

nv(rips::AbstractRipsFiltration) = size(adjacency_matrix(rips), 1)

function births(rips::AbstractRipsFiltration)
    adj = adjacency_matrix(rips)
    return view(adj, diagind(adj))
end

simplex_type(::Type{<:AbstractRipsFiltration{I,T}}, D) where {I,T} = Simplex{D,T,I}

function edges(rips::AbstractRipsFiltration)
    if issparse(adjacency_matrix(rips))
        _sparse_edges(rips)
    else
        _full_edges(rips)
    end
end

function _full_edges(rips::AbstractRipsFiltration)
    result = edge_type(rips)[]
    @inbounds for u in 1:nv(rips), v in (u + 1):nv(rips)
        sx = unsafe_simplex(rips, Val(1), (v, u))
        !isnothing(sx) && push!(result, sx)
    end
    return result
end

function _sparse_edges(rips::AbstractRipsFiltration)
    adj = adjacency_matrix(rips)
    result = edge_type(rips)[]
    rows = rowvals(adj)
    for u in 1:nv(rips)
        for i in nzrange(adj, u)
            v = rows[i]
            if v > u
                sx = unsafe_simplex(rips, Val(1), (v, u))
                !isnothing(sx) && push!(result, sx)
            end
        end
    end
    return result
end

@inline @propagate_inbounds function unsafe_simplex(
    ::Type{S}, rips::AbstractRipsFiltration{I,T}, vertices
) where {I,T,S<:Simplex}
    if dim(S) == 0
        return S(vertices[1], births(rips)[vertices[1]])
    else
        adj = adjacency_matrix(rips)
        n = length(vertices)
        diameter = typemin(T)
        for i in 1:n, j in (i + 1):n
            e = adj[vertices[j], vertices[i]]
            (iszero(e) || e > threshold(rips)) && return nothing
            diameter = ifelse(e > diameter, e, diameter)
        end
        return S(index(vertices), diameter)
    end
end

@inline @propagate_inbounds function unsafe_cofacet(
    ::Type{S}, rips::AbstractRipsFiltration{I,T}, simplex, cofacet_vertices, new_vertex
) where {I,T,S<:Simplex}
    diameter = birth(simplex)
    adj = adjacency_matrix(rips)
    for v in cofacet_vertices
        v == new_vertex && continue
        # Even though this looks like a tight loop, v changes way more often than us, so
        # this is the faster order of indexing by new_vertex and v.
        e = adj[new_vertex, v]
        (iszero(e) || e > threshold(rips)) && return nothing
        diameter = ifelse(e > diameter, e, diameter)
    end
    return S(index(cofacet_vertices), diameter)
end

# This definition is used in SparseCoboundary
@inline @propagate_inbounds function unsafe_cofacet(
    ::Type{S}, rips::AbstractRipsFiltration{I}, simplex, cofacet_vertices, _, new_edges
) where {I,S<:Simplex}
    diameter = birth(simplex)
    for e in new_edges
        e > threshold(rips) && return nothing
        diameter = ifelse(e > diameter, e, diameter)
    end
    return S(index(cofacet_vertices), diameter)
end

function coboundary(
    rips::AbstractRipsFiltration, sx::AbstractSimplex, ::Val{A}=Val(true)
) where {A}
    if adjacency_matrix(rips) isa SparseMatrixCSC
        return SparseCoboundary{A}(rips, sx)
    else
        return Coboundary{A}(rips, sx)
    end
end

distance_matrix(rips::AbstractRipsFiltration) = adjacency_matrix(rips)

function _check_distance_matrix(dist::SparseMatrixCSC)
    n = LinearAlgebra.checksquare(dist)
    vals = nonzeros(dist)
    rows = rowvals(dist)
    for i in 1:n
        vertex_birth = dist[i, i]
        for j in nzrange(dist, i)
            i == rows[j] && continue
            edge_birth = vals[j]
            iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))
            edge_birth < vertex_birth && throw(
                ArgumentError(
                    "edges with birth value lower than vertex births in input matrix"
                ),
            )
            edge_birth ≠ dist[i, rows[j]] &&
                throw(ArgumentError("input matrix not symmetric"))
        end
    end
end

function _check_distance_matrix(dist::AbstractMatrix)
    n = LinearAlgebra.checksquare(dist)
    for i in 1:n
        for j in (i + 1):n
            vertex_birth = max(dist[j, j], dist[i, i])
            edge_birth = dist[j, i]
            iszero(edge_birth) && throw(ArgumentError("zero edges in input matrix"))
            edge_birth < vertex_birth && throw(
                ArgumentError(
                    "edges with birth value lower than vertex births in input matrix"
                ),
            )
            edge_birth ≠ dist[i, j] && throw(ArgumentError("input matrix not symmetric"))
        end
    end
end

"""
    Rips{I, T} <: AbstractRipsFiltration{I, T}

This type represents a filtration of Vietoris-Rips complexes.

Diagonal items in the input matrix are treated as vertex birth times.

Zero values are not allowed due to how sparse matrices work in Julia. If you need zero birth
times, try offseting all values by a constant.

Threshold defaults to the radius of the input space. When using low `threshold`s, consider
using the `sparse=true` keyword argument. It will give the same result, but may be much
faster.

# Constructors

* `Rips(distance_matrix; threshold=nothing)`
* `Rips(points; metric=Euclidean(1e-12), threshold=nothing)`
* `Rips{I}(args...)`: `I` sets the size of integer used to represent simplices.

# Examples

```jldoctest
julia> data = [(sin(t), cos(t)) for t in range(0, 2π, length=101)][1:end-1];

julia> ripserer(Rips(data))
2-element Vector{PersistenceDiagramsBase.PersistenceDiagram}:
 100-element 0-dimensional PersistenceDiagram
 1-element 1-dimensional PersistenceDiagram

julia> ripserer(Rips(data, threshold=1.7))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0628, ∞)

julia> using Distances

julia> ripserer(Rips(data, metric=Cityblock()))[2]
1-element 1-dimensional PersistenceDiagram:
 [0.0888, 2.0)

```
"""
struct Rips{I,T,A<:AbstractMatrix{T}} <: AbstractRipsFiltration{I,T}
    adj::A
    threshold::T
end

function Rips{I}(
    dists::AbstractMatrix{T}; threshold=nothing, sparse=false, verbose=false
) where {I,T}
    _check_distance_matrix(dists)
    thresh = isnothing(threshold) ? radius(dists) : T(threshold)
    if sparse
        adj = SparseArrays.sparse(dists)
    else
        adj = copy(dists)
    end
    if issparse(adj)
        adj isa SparseMatrixCSC ||
            throw(ArgumentError("only SparseMatrixCSC sparse matrices supported"))
        SparseArrays.fkeep!(adj, (_, _, v) -> v ≤ thresh)
    end
    return Rips{I,T,typeof(adj)}(adj, thresh)
end
function Rips{I}(points::AbstractVector; metric=Euclidean(1e-12), kwargs...) where {I}
    if !allunique(points)
        @warn "points not unique"
        points = unique(points)
    end
    return Rips{I}(distances(points, metric); kwargs...)
end
function Rips(dist_or_points; kwargs...)
    return Rips{Int}(dist_or_points; kwargs...)
end

function Base.show(io::IO, rips::Rips{I,T}) where {I,T}
    return print(io, "Rips{$I, $T}(nv=$(nv(rips)), sparse=$(issparse(rips.adj)))")
end

threshold(rips::Rips) = rips.threshold
adjacency_matrix(rips::Rips) = rips.adj


text:
using Test, Documenter
using MeshArrays

MeshArrays.GRID_LL360_download()
MeshArrays.GRID_LLC90_download()
MeshArrays.GRID_CS32_download()

p=dirname(pathof(MeshArrays))
include(joinpath(p,"../examples/Demos.jl"))

@testset "MeshArrays tests:" begin
    for nTopo=1:4
        if nTopo==1; grTopo="CubeSphere"; nFaces=6; N=200;
        elseif nTopo==2; grTopo="LatLonCap"; nFaces=5; N=200;
        elseif nTopo==3; grTopo="PeriodicChannel"; nFaces=1; N=400;
        elseif nTopo==4; grTopo="PeriodicDomain"; nFaces=1; N=400;
        end;
        Npt=nFaces*N*N
        γ,Γ=MeshArrays.GridOfOnes(grTopo,nFaces,N;option="full")
        @test γ.class == grTopo
        Rini= 0.; Rend= 0.;
        (Rini,Rend,DXCsm,DYCsm)=demo2(Γ);
        @test isa(Rend,MeshArray)
        @test sum(isfinite.(Rend)) == Npt
        Sini=sqrt(sum(Rini*Rini)/(Npt-1.0))
        Send=sqrt(sum(Rend*Rend)/(Npt-1.0))
        #println([Sini Send])
        @test isapprox(Sini,1.000; atol=1e-2)
        @test isapprox(Send,0.093; atol=1e-2)
        (dRdx,dRdy)=gradient(Rend,Γ)
        exchange(Rend,2)
        exchange(dRdx,dRdy,1)
        (dRdx_e,dRdy_e)=exchange(dRdx,dRdy,1)
        @test isa(dRdx_e,MeshArray)
    end
end

@testset "Vertical Dimension:" begin
    γ=GridSpec("PeriodicChannel",MeshArrays.GRID_LL360)
    Γ=GridLoad(γ;option="full")
    θ=Float64.(Γ.hFacC)
    nk=length(Γ.RC)
    #[θ[:,k]=cosd.((nk-k)/nk*Γ.YC) for k in 1:nk]
    #why fail? [θ[:,k]=(nk-k) .+ cosd.(0.5*Γ.YC[:]) for k in 1:nk]
    [θ[1,k]=0.01*(nk-k) .+ cosd.(Γ.YC[1]) for k in 1:nk]
    θ[findall(Γ.hFacC.==0.0)].=NaN
    d=isosurface(θ,1.1,Γ)

    @test isapprox(d[1][180,90],-2204.8384919029777)
end

@testset "Transport computations:" begin
    #Load grid and transport / vector field
    γ=GridSpec("LatLonCap",MeshArrays.GRID_LLC90)
    Tx=γ.read(MeshArrays.GRID_LLC90*"TrspX.bin",MeshArray(γ,Float32))
    Ty=γ.read(MeshArrays.GRID_LLC90*"TrspY.bin",MeshArray(γ,Float32))
    Γ=GridLoad(γ;option="full")

    μ =Γ.hFacC[:,1]
    μ[findall(μ.>0.0)].=1.0
    μ[findall(μ.==0.0)].=NaN

    #Meridional transport integral
    uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"])
    L=-85.0:5.0:85.0; LC=LatitudeCircles(L,Γ)
    T=Array{Float64,1}(undef,length(LC))
    [T[i]=1e-6*ThroughFlow(uv,LC[i],Γ) for i=1:length(LC)]

    #See: OceanTransports/helper_functions.jl
    #u,v,uC,vC=rotate_uv(uv,Γ);

    #Transport convergence
    TrspCon=μ.*convergence(Tx,Ty)
    #Scalar potential
    TrspPot=ScalarPotential(TrspCon)
    #Divergent transport component
    (TxD,TyD)=gradient(TrspPot,Γ)
    TxD=TxD.*Γ.DXC
    TyD=TyD.*Γ.DYC
    #Rotational transport component
    TxR = Tx-TxD
    TyR = Ty-TyD
    #vector Potential
    TrspPsi=VectorPotential(TxR,TyR,Γ);

    tst=extrema(write(TrspPsi))

    @test prod(isapprox.(tst,(-8.08f7,2.19f8)))
end

@testset "gcmfaces type:" begin
    γ=GridSpec("CubeSphere",MeshArrays.GRID_CS32)

    MeshArrays.gcmfaces(γ)
    MeshArrays.gcmfaces(γ,Float32)
    tmp=MeshArrays.gcmfaces(γ,Float32,3)
    tmp[3,1,2]=1.0
    view(tmp,:,1,1:2)
    MeshArrays.gcmfaces(γ,tmp.f)
    MeshArrays.fsize(tmp)
    MeshArrays.fsize(tmp,2)
    MeshArrays.fsize(tmp.f)
    MeshArrays.fsize(tmp.f,2)
    size(tmp)
    size(tmp,3)
    show(tmp)

    x=tmp[1:10,1,1:2]; y=x[2]; x[3]=1.0
    view(x,1:3,:,1)
    MeshArrays.gcmsubset(γ,x.f,x.fSize,x.aSize,x.i,x.iSize)
    MeshArrays.fsize(x.f)
    MeshArrays.fsize(x.f,2)
    size(x)
    size(x,3)

    MeshArrays.fijind(tmp,10)
    MeshArrays.fsize(tmp)
    @test isa(tmp,MeshArrays.gcmfaces)

    tmp=MeshArray(γ)
    tmp1=findall(tmp.>0)
    tmp[tmp1]
    tmp[tmp1].=1.0
    size(tmp1)
    tmp1[2]
    view(tmp1,:)
    show(tmp1)
    similar(tmp1)
    #tmp[tmp1]

    show(tmp)
    MeshArrays.getindexetc(tmp,2)
    MeshArray(γ,tmp.f,meta=tmp.meta)
    MeshArray(γ,meta=tmp.meta)
    MeshArray(γ,Float32,3,4)
    MeshArray(γ,Float32,tmp.fSize,tmp.fIndex,2,3)

    tmp1=MeshArray(γ,Float32,3)
    MeshArray(γ,tmp1.f,meta=tmp1.meta)
    MeshArrays.getindexetc(tmp1,2,1)
end

@testset "doctests" begin
    doctest(MeshArrays; manual = false)
end


text:
import Base: isnan
import DataFrames: DataFrame, DataArray, isna, dropna, array, ncol, convert, NA, NAtype

export GLRM, observations, expand_categoricals!, NaNs_to_NAs!, NAs_to_0s!

include("fit_dataframe_w_type_imputation.jl")

probabilistic_losses = Dict{Symbol, Any}(
    :real        => QuadLoss,
    :bool        => LogisticLoss,
    :ord         => MultinomialOrdinalLoss,
    :cat         => MultinomialLoss
)

function GLRM(df::DataFrame, k::Int, datatypes::Array{Symbol,1};
              loss_map = probabilistic_losses,
              rx = QuadReg(.01), ry = QuadReg(.01),
              offset = true, scale = false, prob_scale = true,
              transform_data_to_numbers = true, NaNs_to_NAs = true)

    # check input
    if ncol(df)!=length(datatypes)
        error("third argument (datatypes) must have one entry for each column of data frame.")
    end
    # validate input
    # XXX check for undefined entry?
    for dt in datatypes
        if !(dt in keys(loss_map))
            error("data types must be either :real, :bool, :ord, or :cat, not $dt")
        end
    end

    # clean up dataframe if needed
    A = copy(df)
    if NaNs_to_NAs
        NaNs_to_NAs!(A)
    end

    # define loss functions for each column
    losses = Array(Loss, ncol(A))
    for j=1:ncol(df)
        losstype = loss_map[datatypes[j]]
        if transform_data_to_numbers
            map_to_numbers!(A, j, datatypes[j])
        end
        losses[j] = pick_loss(losstype, A[j])
    end

    # identify which entries in data frame have been observed (ie are not N/A)
    obs = observations(df)

    # form model
    rys = Array(Regularizer, length(losses))
    for i=1:length(losses)
        if isa(losses[i], MultinomialOrdinalLoss) || isa(losses[i], OrdisticLoss)
            rys[i] = OrdinalReg(copy(ry))
        else
            rys[i] = copy(ry)
        end
    end
    glrm = GLRM(A, losses, rx, rys, k, obs=obs, offset=offset, scale=scale)

    # scale model so it really computes the MAP estimator of the parameters
    if prob_scale
        prob_scale!(glrm)
    end

    return glrm
end

## transform data to numbers

function map_to_numbers!(df, j::Int, datatype::Symbol)
    # easy case
    if datatype == :real
        if all(xi -> isa(xi, Number), df[j][!isna(df[j])])
            return df[j]
        else
            error("column contains non-numerical values")
        end
    end

    # harder cases
    col = copy(df[j])
    levels = Set(col[!isna(col)])
    if datatype == :bool
        if length(levels)>2
            error("Boolean variable should have at most two levels; instead, got:\n$levels")
        end
        colmap = Dict{Any,Int}(zip(sort(collect(levels)), [-1,1][1:length(levels)]))
    elseif datatype == :cat || datatype == :ord
        colmap = Dict{Any,Int}(zip(sort(collect(levels)), 1:length(levels)))
    else
        error("datatype $datatype not recognized")
    end
    df[j] = DataArray(Int, length(df[j]))
    for i in 1:length(col)
        if !isna(col[i])
            df[j][i] = colmap[col[i]]
        end
    end
    return df[j]
end

function map_to_numbers!(df, j::Int, loss::Type{QuadLoss})
    if all(xi -> isa(xi, Number), df[j][!isna(df[j])])
        return df[j]
    else
        error("column contains non-numerical values")
    end
end

function map_to_numbers!(df, j::Int, loss::Type{LogisticLoss})
    col = copy(df[j])
    levels = Set(col[!isna(col)])
    if length(levels)>2
        error("Boolean variable should have at most two levels")
    end
    colmap = Dict{Any,Int}(zip(sort(collect(levels)), [-1,1][1:length(levels)]))
    df[j] = DataArray(Int, length(df[j]))
    for i in 1:length(col)
        if !isna(col[i])
            df[j][i] = colmap[col[i]]
        end
    end
    return df[j]
end

function map_to_numbers!(df, j::Int, loss::Type{MultinomialLoss})
    col = copy(df[j])
    levels = Set(col[!isna(col)])
    colmap = Dict{Any,Int}(zip(sort(collect(levels)), 1:length(levels)))
    df[j] = DataArray(Int, length(df[j]))
    for i in 1:length(col)
        if !isna(col[i])
            df[j][i] = colmap[col[i]]
        end
    end
    return df[j]
end

function map_to_numbers!(df, j::Int, loss::Type{MultinomialOrdinalLoss})
    col = copy(df[j])
    levels = Set(col[!isna(col)])
    colmap = Dict{Any,Int}(zip(sort(collect(levels)), 1:length(levels)))
    df[j] = DataArray(Int, length(df[j]))
    for i in 1:length(col)
        if !isna(col[i])
            df[j][i] = colmap[col[i]]
        end
    end
    return df[j]
end

## sanity check the choice of loss

# this default definition could be tighter: only needs to be defined for arguments of types that subtype Loss
function pick_loss(l, col)
    return l()
end

function pick_loss(l::Type{LogisticLoss}, col)
    if all(xi -> isna(xi) || xi in [-1,1], col)
        return l()
    else
        error("LogisticLoss can only be used on data taking values in {-1, 1}")
    end
end

function pick_loss(l::Type{MultinomialLoss}, col)
    if all(xi -> isna(xi) || (isa(xi, Int) && xi >= 1), col)
        return l(maximum(col[!isna(col)]))
    else
        error("MultinomialLoss can only be used on data taking positive integer values")
    end
end

function pick_loss(l::Type{MultinomialOrdinalLoss}, col)
    if all(xi -> isna(xi) || (isa(xi, Int) && xi >= 1), col)
        return l(maximum(col[!isna(col)]))
    else
        error("MultinomialOrdinalLoss can only be used on data taking positive integer values")
    end
end

observations(da::DataArray) = df_observations(da)
observations(df::DataFrame) = df_observations(df)
function df_observations(da)
    obs = @compat Tuple{Int, Int}[]
    m,n = size(da)
    for j=1:n # follow column-major order. First element of index in innermost loop
        for i=1:m
            if !isna(da[i,j])
                push!(obs,(i,j))
            end
        end
    end
    return obs
end

# NAs in the data frame will be replaced by the number `z`
function df2array(df::DataFrame, z::Number)
    A = zeros(size(df))
    for i=1:size(A,2)
        if issubtype(typeof(df[:,i]), Array)
            A[:,i] = df[:,i]
        elseif typeof(df[i]) == Bool
            A[:,i] = convert(Array, (2*df[i]-1), z)
        else
            A[:,i] = convert(Array, df[i], z)
        end
    end
    return A
end
df2array(df::DataFrame) = df2array(df, 0)

function get_reals(df::DataFrame)
    m,n = size(df)
    reals = [typeof(df[i])<:AbstractArray{Float64,1} for i in 1:n]
    n1 = sum(reals)
    losses = Array(Loss,n1)
    for i=1:n1
        losses[i] = QuadLoss()
    end
    return reals, losses
end

function get_bools(df::DataFrame)
    m,n = size(df)
    bools = [isa(df[i], AbstractArray{Bool,1}) for i in 1:n]
    n1 = sum(bools)
    losses = Array(Loss,n1)
    for i=1:n1
        losses[i] = HingeLoss()
    end
    return bools, losses
end

function get_ordinals(df::DataFrame)
    m,n = size(df)
    # there must be a better way to check types...
    ordinals = [(isa(df[i], AbstractArray{Int,1}) ||
                 isa(df[i], AbstractArray{Int32,1}) ||
                 isa(df[i], AbstractArray{Int64,1})) for i in 1:n]
    nord = sum(ordinals)
    ord_idx = (1:size(df,2))[ordinals]
    maxs = zeros(nord,1)
    mins = zeros(nord,1)
    for i in 1:nord
        col = df[ord_idx[i]]
        try
            maxs[i] = maximum(dropna(col))
            mins[i] = minimum(dropna(col))
        end
    end

    # set losses and regularizers
    losses = Array(Loss,nord)
    for i=1:nord
        losses[i] = OrdinalHingeLoss(mins[i],maxs[i])
    end
    return ordinals, losses
end

# expand categorical columns, given as column indices, into one boolean column for each level
function expand_categoricals!(df::DataFrame,categoricals::Array{Int,1})
    # map from names to indices; not used: categoricalidxs = map(y->df.colindex[y], categoricals)
    # create one boolean column for each level of categorical column
    colnames = names(df)
    for col in categoricals
        levels = sort(unique(df[:,col]))
        for level in levels
            if !isna(level)
                colname = symbol(string(colnames[col])*"="*string(level))
                df[colname] = (df[:,col] .== level)
            end
        end
    end
    # remove the original categorical columns
    for cat in sort(categoricals, rev=true)
        delete!(df, cat)
    end
    return df
end
function expand_categoricals!(df::DataFrame,categoricals::UnitRange{Int})
    expand_categoricals!(df, Int[i for i in categoricals])
end
# expand categoricals given as names of columns rather than column indices
function expand_categoricals!(df::DataFrame,categoricals::Array)
    # map from names to indices
    categoricalidxs = map(y->df.colindex[y], categoricals)
    return expand_categoricals!(df, categoricalidxs)
end

# convert NaNs to NAs
isnan(x::NAtype) = false
isnan(x::AbstractString) = false
function NaNs_to_NAs!(df::DataFrame)
    m,n = size(df)
    for j=1:n # follow column-major order. First element of index in innermost loop
        for i=1:m
            if isnan(df[i,j])
                df[i,j] = NA
            end
        end
    end
    return df
end

function NAs_to_0s!(df::DataFrame)
    m,n = size(df)
    for j=1:n # follow column-major order. First element of index in innermost loop
        for i=1:m
            if isna(df[i,j])
                df[i,j] = 0
            end
        end
    end
    return df
end


text:
"""
Given parEst, an estimate the ergm parameters, sample from that ergm nSample times and re-estimate the same ergm on the obtained sample of observations.

Return a distribution of estimated parameter, aka parametric bootstrap following "Exponential Random Graph Models with Big Networks: Maximum Pseudolikelihood Estimation and the Parametric Bootstrap" of Christian S. Schmid, Bruce A. Desmarais.
"""
function get_par_boot_ergm_distrib(model, parEst, N; nSample = 100)

    mat_sample = StaticNets.sample_ergm(model, N, parEst, nSample)

    estimates = zeros(length(parEst), nSample)

   for i in 1:nSample
        estimates[:,i] = StaticNets.estimate(model,  mat_sample[:,:,i])
    end
    return estimates
end
export get_par_boot_ergm_distrib

"""
Get the confidence intervals from parametric bootstrap, as described in "Exponential Random Graph Models with Big Networks: Maximum Pseudolikelihood Estimation and the Parametric Bootstrap" of Christian S. Schmid, Bruce A. Desmarais.
"""
function get_conf_int_ergm_par_boot(model, N, parEst, quantilesVals; nSample = 100)

    par_b_dist = get_par_boot_ergm_distrib(model, parEst, N; nSample = nSample)

    confInt = reduce(vcat, [quantile(par_b_dist[p,:], quantilesVals)' for p in 1:model.nErgmPar])
end
export get_conf_int_ergm_par_boot

is_between(x, interval) = interval[1] <= x <= interval[2]

function get_coverage_conf_int_par_boot_ergm(model, N, parDgp, quantilesVals; nSampleDgp=100, nSampleBands = 100)

    # sample from dgp
    matsSample = StaticNets.sample_ergm(model, N, parDgp, nSampleDgp)

    estimates = reduce(hcat,[StaticNets.estimate(model, A) for A in eachslice(matsSample, dims=3)])

    coverFlags =falses(size(estimates))

    for n in 1:nSampleDgp
        confInt = get_conf_int_ergm_par_boot(model, N, estimates[:, n], quantilesVals; nSample=nSampleBands)

        for p in 1:model.nErgmPar
            coverFlags[p,n] = is_between(parDgp[p], confInt[p,:])
        end
    end

    return coverFlags
end
export get_coverage_conf_int_par_boot_ergm


text:
#=

Predictor

## Introduction

In an ideal world the predictor only uses an order (2,1) Padé approximant to generate
an initial guess for the Newton corrector.
The choice of Padé approximants as predictors is motived by the results in [1].
In particular, the ability to provide a trust region of the prediction method to
significantly decrease the likelihood of path jumping.
additionally, sice Padé approximants are rational functions, these perform very well
for diverging paths and close to singularities.

Now, for an Padé approximant of order (2,1) the first 3 derivatives are needed.
The most robust way to achieve this it to use automatic differentiation.
However, due to significant compile times, we also have to consider the case
that the Padé approximants are computed by numerical differentiation using the scheme
derived in [2]. The numerical differentiation scheme only uses O(h^2) approximations
since we are mostly constrained by the available precision.

However, the smaller the trust region is (i.e. the closer we are to a singularity) the
smaller we need to take the discretization in the numerical differentiation.
Thus, the error in the numerical differentiation is inversely corrolated to the trust
region.
Additionally, the third derivative needs to be sufficiently accurate. Otherwise
the predictor doesn't generate sufficiently accurate intial guesses resulting
in very small step sizes.
Therefore it can necessary to drop down to an Padé approximant of order (1,1).
Usually the arithmetic is sufficient to obtain for the second derivative a sufficiently
accurate approximation.

If also the second derivative is not available we use cubic hermite interpolation.
For diverging paths this works much much better than (explicit) Runge-Kutta methods.
The reason that we do not use expl. RK methods is that, then when numerical differentiation
has to give up, we are in a regime where the ODE of the path is *stiff*. Thus, expl.
RK just don't yield satisfying results.

## Endgame and coordinates

Additionally, the predictor computes in s-coordinates, that is for the path x(t)
the predictor internally computes with y(s) = x(s^m) where m is the *winding number*.
If m = 1 then this is the same as computing with x(t). However, during the endgame
`m` is possibly set.
All APIs expect input in t-coordinates and return results in t-coordinates.

[1]: A Robust Numerical Path Tracking Algorithm for Polynomial Homotopy Continuation,
 Telen, Van Barel, Verschelde https://arxiv.org/abs/1909.04984

[2]:Mackens, Wolfgang. "Numerical differentiation of implicitly defined space curves."
   Computing 41.3 (1989): 237-260.
=#

module PredictionMethod
@enum methods begin
    Pade21
    Hermite
end
end
using .PredictionMethod: PredictionMethod

"""
    Predict(H::AbstractHomotopy)

Construct a predictor. The predicor uses a Padé approximant of order (2,1) to produce
an initial guess for the Newton corrector.
For this, the predictor tries to compute the local derivatives up to order 3.
This is done using automatic differentiation for the derivatives up to order `3`.
Otherwise numerical differentiation is used.
"""
Base.@kwdef mutable struct Predictor
    method::PredictionMethod.methods = PredictionMethod.Pade21
    order::Int = 4
    use_hermite::Bool = true
    trust_region::Float64 = Inf
    local_error::Float64 = Inf
    cond_H_ẋ::Float64 = Inf

    tx⁰::TaylorVector{1,ComplexF64}
    tx¹::TaylorVector{2,ComplexF64}
    tx²::TaylorVector{3,ComplexF64}
    tx³::TaylorVector{4,ComplexF64}
    t::ComplexF64 = complex(NaN)
    tx_norm::Vector{Float64} = zeros(4)

    # finite diff
    xtemp::Vector{ComplexF64}
    u::Vector{ComplexF64}
    u₁::Vector{ComplexF64}
    u₂::Vector{ComplexF64}

    # for interpolation
    prev_tx¹::TaylorVector{2,ComplexF64}
    prev_t::ComplexF64 = complex(NaN)

    # endgame predictor
    winding_number::Int = 1
    s::ComplexF64 = complex(NaN)
    prev_s::ComplexF64 = complex(NaN)
    ty¹::TaylorVector{2,ComplexF64}
    prev_ty¹::TaylorVector{2,ComplexF64}
end

function Predictor(H::AbstractHomotopy)
    m, n = size(H)
    tx³ = TaylorVector{4}(ComplexF64, n)
    Predictor(
        tx⁰ = TaylorVector{1}(tx³),
        tx¹ = TaylorVector{2}(tx³),
        tx² = TaylorVector{3}(tx³),
        tx³ = tx³,
        xtemp = zeros(ComplexF64, n),
        u = zeros(ComplexF64, m),
        u₁ = zeros(ComplexF64, m),
        u₂ = zeros(ComplexF64, m),
        prev_tx¹ = TaylorVector{2}(ComplexF64, n),
        ty¹ = TaylorVector{2}(ComplexF64, n),
        prev_ty¹ = TaylorVector{2}(ComplexF64, n),
    )
end

function init!(predictor::Predictor)
    predictor.cond_H_ẋ = 1.0
    predictor.winding_number = 1

    predictor.t = predictor.prev_t = NaN
    predictor.s = predictor.prev_s = NaN
    predictor.trust_region = predictor.local_error = NaN
    predictor.use_hermite = true
    predictor
end

order(predictor::Predictor) = predictor.order
local_error(predictor::Predictor) = predictor.local_error
trust_region(predictor::Predictor) = predictor.trust_region

function winding_number!(P::Predictor, m)
    P.winding_number = m
    P
end

"""
    update!(
        P::Predictor,
        H::AbstractHomotopy,
        x,
        t,
        jacobian::Jacobian,
        norm::AbstractNorm,
        x̂ = nothing,
    )

Upadte the predictor with the new solution `(x,t)` of `H(x,t) = 0`.
This computes new local derivatives and chooses an appriopriate prediction method.
"""
function update!(
    predictor::Predictor,
    H::AbstractHomotopy,
    x,
    t,
    J::Jacobian,
    norm::AbstractNorm,
    x̂ = nothing, # the predicted value
)
    # The general strategy is as follows:
    #
    # If m = winding_number > 1 then we only use a hermite predictor.
    #
    # Otherwise we use a Padé predictor of order (2,1). For this we need the first three
    # derivativatives of the path x(t).
    # We obtain accurate derivatives by using automatic differentiation (AD) with the taylor!
    # api.

    @unpack u, tx⁰, tx¹, tx², tx³, ty¹, xtemp, tx_norm = predictor

    x⁰, x¹, x², x³ = vectors(tx³)
    y⁰, y¹ = vectors(ty¹)

    m = predictor.winding_number
    @inbounds for i = 1:length(xtemp)
        predictor.prev_tx¹[i, 1] = predictor.tx¹[i, 1]
        predictor.prev_tx¹[i, 2] = predictor.tx¹[i, 2]
    end

    predictor.prev_t, predictor.t = predictor.t, t

    if m > 1
        predictor.prev_s = predictor.s
        predictor.s = t_to_s_plane(t, m)
    end

    # we use the accuracy of the previous prediction for the local error estimate
    if isnothing(x̂)
        predictor.local_error = NaN
    else
        Δs = fast_abs(t - predictor.prev_t)
        predictor.local_error = norm(x̂, x) / Δs^predictor.order
    end

    x⁰ .= x
    tx_norm[1] = norm(x)
    if m > 1
        y⁰ .= x
    end

    # Compute ẋ always using AD
    taylor!(u, Val(1), H, x, t, true)
    # @show u[1:2]
    u .= .-u
    LA.ldiv!(xtemp, J, u)
    # Check error made in the linear algebra
    δ = fixed_precision_iterative_refinement!(xtemp, workspace(J), u, norm)
    predictor.cond_H_ẋ = cond_H_ẋ = δ / eps()
    # @show cond_H_ẋ
    tol_δ₁ = 1e-10

    δ > tol_δ₁ && iterative_refinement!(xtemp, J, u, InfNorm(); tol = tol_δ₁, max_iters = 5)
    tx_norm[2] = norm(xtemp)

    x¹ .= xtemp
    if m > 1
        if m == 2
            μ = 2 * predictor.s
        else
            μ = m * predictor.s^(m - 1)
        end
        y¹ .= μ .* xtemp

        predictor.method = PredictionMethod.Hermite
        predictor.order = 4
        predictor.trust_region = tx_norm[1] / tx_norm[2]
        if isnan(predictor.local_error)
            predictor.local_error = (tx_norm[2] / tx_norm[1])^3
        end
        return predictor
    end

    taylor!(u, Val(2), H, tx¹, t, true)
    u .= .-u
    LA.ldiv!(xtemp, J, u)
    tol_δ₂ = 1e-10
    δ > tol_δ₂ && iterative_refinement!(xtemp, J, u, norm; tol = tol_δ₂, max_iters = 4)
    tx_norm[3] = norm(xtemp)
    x² .= xtemp

    taylor!(u, Val(3), H, tx², t, true)
    u .= .-u
    LA.ldiv!(xtemp, J, u)
    tol_δ₃ = 1e-4
    δ > tol_δ₃ && iterative_refinement!(xtemp, J, u, norm; tol = tol_δ₃, max_iters = 3)
    tx_norm[4] = norm(xtemp)
    x³ .= xtemp

    τ = Inf
    for (i, (x, x¹, x², x³)) in enumerate(tx³)
        c, c¹, c², c³ = fast_abs.((x, x¹, x², x³))
        λ = max(1e-6, c¹)
        c¹ /= λ
        c² /= λ^2
        c³ /= λ^3
        tol = 1e-14 * max(c¹, c², c³)
        if !((c¹ ≤ tol && c² ≤ tol && c³ ≤ tol) || c² ≤ tol)
            τᵢ = (c² / c³) / λ
            if τᵢ < τ
                τ = τᵢ
            end
        end
    end
    if !isfinite(τ)
        τ = tx_norm[3] / tx_norm[4]
    end
    if !isfinite(τ)
        τ = tx_norm[1] / maximum(tx_norm)
    end

    predictor.method = PredictionMethod.Pade21
    predictor.order = 4
    predictor.trust_region = τ
    if isnan(predictor.local_error)
        predictor.local_error = ((inv(τ))^2)^2
    end

    predictor
end

predict!(x̂, pred::Predictor, H, x, t, Δt) = predict!(x̂, pred, pred.method, H, x, t, Δt)
function predict!(x̂, pred::Predictor, method::PredictionMethod.methods, H, x, t, Δt)
    m = pred.winding_number

    if method == PredictionMethod.Pade21
        λ = pred.trust_region
        λ² = λ^2
        λ³ = λ^3
        tol = 1e-12
        for (i, (x, x¹, x², x³)) in enumerate(pred.tx³)
            c, c¹, c², c³ = fast_abs.((x, x¹, x², x³))
            # check if only taylor series is used
            τ = tol * √(c^2 + (c¹ * λ)^2 + (c² * λ²)^2 + (c³ * λ³)^2)
            if c³ * λ³ ≤ τ || c² * λ² ≤ τ
                x̂[i] = x + Δt * (x¹ + Δt * x²)
            else
                δᵢ = 1 - Δt * x³ / x²
                x̂[i] = x + Δt * (x¹ + Δt * x² / δᵢ)
            end
        end
    elseif method == PredictionMethod.Hermite
        prev_s = t_to_s_plane(pred.prev_t, m)
        s = t_to_s_plane(t, m)
        s′ = t_to_s_plane(t + Δt, m)
        if m == 2
            prev_sm = 2 * prev_s
            sm = 2 * s
        else
            prev_sm = m * prev_s^(m - 1)
            sm = m * s^(m - 1)
        end
        for i = 1:length(pred.prev_tx¹)
            pred.prev_ty¹[i, 1] = pred.prev_tx¹[i, 1]
            pred.prev_ty¹[i, 2] = prev_sm * pred.prev_tx¹[i, 2]
        end
        for i = 1:length(pred.prev_tx¹)
            pred.ty¹[i, 1] = pred.tx¹[i, 1]
            pred.ty¹[i, 2] = sm * pred.tx¹[i, 2]
        end
        cubic_hermite!(x̂, pred.prev_ty¹, prev_s, pred.ty¹, s, s′)
    end

    x̂
end

"""
    t_to_s_plane(t::Complex, m::Int; branch::Int = 0)

Convert ``t = r * \\exp(√(-1)θ)`` to ``t = r^{1/m} * \\exp(√(-1)θ/m)`` where for
``r^{1/m}`` the positive real root is choosen.
The `branch` keyword allows to choose a different branch cut.
"""
function t_to_s_plane(t::Complex, m::Int; branch::Int = 0)
    r = fast_abs(t)
    if isreal(t) && real(t) > 0
        if branch > 0
            complex(nthroot(r, m), 2branch * π / m)
        else
            complex(nthroot(r, m))
        end
    else
        # angle returns a value in [-π,π], but we want [0,π]
        θ = mod2pi(angle(t)) + branch * 2π
        return nthroot(r, m) * cis(θ / m)
    end
end

function cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    if isreal(t₀) && isreal(t₁) && isreal(t)
        _cubic_hermite!(x̂, tx¹₀, real(t₀), tx¹₁, real(t₁), real(t))
    else
        _cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    end
end
@inline function _cubic_hermite!(x̂, tx¹₀, t₀, tx¹₁, t₁, t)
    s = (t - t₀) / (t₁ - t₀)
    h₀₀ = (1 + 2s) * (1 - s)^2
    h₁₀ = (t - t₀) * (1 - s)^2
    h₀₁ = s^2 * (3 - 2s)
    h₁₁ = (t - t₀) * s * (s - 1)
    @inbounds for i in eachindex(x̂)
        x̂[i] = h₀₀ * tx¹₀[i, 1] + h₁₀ * tx¹₀[i, 2] + h₀₁ * tx¹₁[i, 1] + h₁₁ * tx¹₁[i, 2]
    end
    x̂
end


text:
set_theory = [
    Axiom("Definition of the universal set", P ⊂ Ω, TRUE),
    Axiom("Definition of the universal set", x ∈ Ω, TRUE),
    Axiom("Definition of the empty set", ∅ ⊂ A, TRUE),
    Axiom("Definition of the empty set", x ∈ ∅, FALSE),
    Axiom("Equal sets are subsets", A ⊂ A, TRUE),

    Axiom("Definition of Element Of", x ∈ (A ∩ B), (x ∈ A) ∧ (x ∈ B)),
    Axiom("Definition of Element Of", (x ∈ A) ∧ (x ∈ B), x ∈ (A ∩ B)),
    Axiom("Definition of Element Of", x ∈ (A ∪ B), (x ∈ A) ∨ (x ∈ B)),
    Axiom("Definition of Element Of", (x ∈ A) ∨ (x ∈ B), x ∈ (A ∪ B)),

    Axiom("Flatten Condition Set", ConditionSet(y, y ∈ A), A),  # todo: change to f(x)
    Axiom("Expand Condition Set", A, ConditionSet(_y, _y ∈ A)),  # todo: change to f(x)
    Axiom("Definition of Subset", A ⊂ B, (_y ∈ A) ⟹ (_y ∈ B)),
    Axiom("Definition of Subset", (_y ∈ A) ⟹ (_y ∈ B), A ⊂ B),

    Axiom("Definition of the universal set", x ∈ (P ∪ Ω), x ∈ Ω),

    #Theorem("Intersections are subsets of Unions", x ∈ (A ∩ B), x ∈ (A ∪ B), Proof([])),
    #Theorem("Intersections are subsets of Unions", (A ∩ B) ⊂ (A ∪ B), TRUE, Proof([])),
]


text:
using AdversarialDriving
using POMDPs, POMDPPolicies, POMDPSimulators
using Distributions
using GridInterpolations, LocalFunctionApproximation, LocalApproximationValueIteration

## setup training params
Np, Nv = 3,3 #25, 15
Nsteps = 2#50

## Construct a disturbance model for the adversaries
da_std = rand(Exponential(0.5))
goal_toggle_p = 10. ^rand(Uniform(-5., -1.))
blinker_toggle_p = 10. ^rand(Uniform(-5., -1.))
v_des = rand(Uniform(15., 30.))
per_timestep_penalty = rand([0, 1e-4, 1e-3, 5e-3])
disturbances = Sampleable[Normal(0,da_std), Bernoulli(goal_toggle_p), Bernoulli(blinker_toggle_p), Bernoulli(0.), Bernoulli(0.)]

## construct the MDP
sut_agent = BlinkerVehicleAgent(rand_up_left(id=1, s_dist=Uniform(15.,25.), v_dist=Uniform(10., 29)), TIDM(Tint_TIDM_template))
right_adv = BlinkerVehicleAgent(rand_right(id=3, s_dist=Uniform(15.,35.), v_dist=Uniform(15., 29.)), TIDM(Tint_TIDM_template), disturbance_model = disturbances)
mdp = DrivingMDP(sut_agent, [right_adv], Tint_roadway, 0.2, γ = 0.95, per_timestep_penalty = per_timestep_penalty, v_des = v_des)

## Solve using local approximation value iteration
goals = Float64.(Tint_goals[laneid(adversaries(mdp)[1].get_initial_entity())])
grid = RectangleGrid(range(0., stop=100.,length=Np), range(0, stop=30., length=Nv), goals, [0.0, 1.0],
                      range(15., stop=100., length=Np), range(0., stop=30., length=Nv), [5.0], [1.0])

interp = LocalGIFunctionApproximator(grid)
solver = LocalApproximationValueIterationSolver(interp, is_mdp_generative = true, n_generative_samples = 5, verbose = true, max_iterations = Nsteps)
policy = solve(solver, mdp)

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
using MolHandler
using Test

include("test_coordinate.jl")
include("test_methods.jl")


text:
# This file was generated, do not modify it. # hide
G1 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks1 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=ParetoRecall()))  # ParetoRecall is the default
@elapsed index!(G1; callbacks=callbacks1)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using Test
using MPI
using CLIMA.Mesh.BrickMesh

function main()
  MPI.Init()
  comm = MPI.COMM_WORLD

  (elemtovert, elemtocorner, facecode) =
    brickmesh((2.0:5.0,4.0:6.0), (false,true);
              part=MPI.Comm_rank(comm)+1,
              numparts=MPI.Comm_size(comm))

  code = centroidtocode(comm, elemtocorner)
  (d, nelem) = size(code)

  root = 0
  counts = MPI.Allgather(Cint(length(code)), comm)
  code_all = MPI.Gatherv(code, counts, root, comm)

  if MPI.Comm_rank(comm) == root
    code_all = reshape(code_all, d, div(sum(counts),d))

    code_expect = UInt64[0x0000000000000000 0x1555555555555555 0xffffffffffffffff 0x5555555555555555 0x6aaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa; 0x0000000000000000 0x5555555555555555 0xffffffffffffffff 0x5555555555555555 0xaaaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa]
    @test code_all == code_expect
  end
end

main()

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# ---- Decision Tree ----------------------------------------------------------

abstract Action
type Decision{T<:Action}
    guards::Vector{Guard}
    action::T
end

type Branch <: Action
    pass::Decision
    fail::Decision
end

type Method <: Action
    args::Vector{Value}
    f::Function
end

code_pass(results::Dict{Node,Any}, a::Branch) = code_dispatch(results, a.pass)
code_fail(results::Dict{Node,Any}, a::Branch) = code_dispatch(results, a.fail)

function code_pass(results::Dict{Node,Any}, a::Method)
    expr(:call, a.f, {evaluate!(results, arg) for arg in a.args}...)
end
code_fail(::Dict{Node,Any}, ::Method) = :(error("no matching pattern found"))

function code_dispatch(results::Dict{Node,Any}, d::Decision)
    pred = quot(true)
    results_fail = copy(results)
    for g in d.guards
        ex = evaluate!(results, g)
        pred = pred == quot(true) ? ex : (:($pred && $ex))
    end
    if pred == quot(true)
        code_pass(results, d.action)
    else
        pass = code_pass(results,      d.action)
        fail = code_fail(results_fail, d.action)
        :( $ex ? $pass : $fail )
    end
end

# ---- MethodTable ------------------------------------------------------------

using Graph

type Sig
    p::Pattern
    gt::Set{Sig}
    lt::Set{Sig}

    Sig(sig::Pattern) = new(sig, Set{Sig}(), Set{Sig}())
end

type MethodTable
    top::Sig
    bottom::Sig
    sigs::Set{Sig}

    function MethodTable()
        top, bottom = Sig(toppat), Sig(nullpat)
        add(top.gt, bottom)
        add(bottom.lt, top)
        new(top, bottom, Set{Sig}(top, bottom))
    end
end

for (below, above, higher_eq) in ((:gt, :lt, >=), (:lt, :gt, <=)=)
    below, above, visit! = quot(below), quot(above), symbol("visit_$(below)!")
    @eval function $visit!(seen::Set{Sig}, s::Sig, at::Sig)
        if has(seen, at); return; end
        add(seen, at)
        if $higher_eq(s.p, at.p)  # s >= at
            if $higher_eq(at.p, s.p); return; end  # s == at
            # s > at
            add(s.($below)), at)
            add(at.($above), s)
        else
            for below in at.($below); $visit!(seen, s, below); end
        end
    end
end

function add(mt::MethodTable, p::Pattern)
    s = Sig(p)
    add(mt.signatures, s) # todo: what if it's already in the DAG?
    visit_gt!(Set{Sig}(), s, mt.top)
    visit_lt!(Set{Sig}(), s, mt.bottom)
    for above in s.lt;  del_each(above.gt, s.gt)  end
    for below in s.gt;  del_each(below.lt, s.lt)  end
end

# ---- Create Decision Tree ---------------------------------------------------

firstitem(iter) = next(iter, start(iter))[1]

function dtree(mt::MethodTable)
    dtree(mt.top, mt.signatures)
end

function visit_gt!(seen::Set{Sig}, s::Sig)
    if has(seen, s); return; end
    add(seen, s)
    for below in s.gt; visit_gt!(seen, below); end
end

function dtree(top::Sig, sigs::Set{Sig})
    pivot = firstitem(top.gt) # should always have at least the bottom below
    if pivot.p === nullpat
        # todo: fill in guards
        # todo: fill method
        Decision(Guard[], Method(Value[], ()->()))
    else
        below = Set{Sig}()
        visit_gt!(below, pivot)
        pass = dtree(pivot, sigs & below)
        fail = dtree(top,   sigs - below)
        # todo: fill in guards
        Decision(Guard[], Branch(pass, fail))
    end
end


text:
module AddDateToTasks

import SearchLight.Migrations: add_column, remove_column

function up()
    add_column(:tasks, :deadline, :date)
end

function down()
    remove_column(:tasks, :deadline)
end

end

max_stars_repo_name: ky-s/MyGenieApp

max_issues_repo_name: ky-s/MyGenieApp

max_forks_repo_name: ky-s/MyGenieApp


text:
struct ODEInterpolatingAdjointSensitivityFunction{C<:AdjointDiffCache,Alg<:InterpolatingAdjoint,
                                                  uType,SType,CPS,pType,fType<:DiffEqBase.AbstractDiffEqFunction} <: SensitivityFunction
  diffcache::C
  sensealg::Alg
  discrete::Bool
  y::uType
  sol::SType
  checkpoint_sol::CPS
  prob::pType
  f::fType
  noiseterm::Bool
end

mutable struct CheckpointSolution{S,I,T,T2}
  cpsol::S # solution in a checkpoint interval
  intervals::I # checkpoint intervals
  cursor::Int # sol.prob.tspan = intervals[cursor]
  tols::T
  tstops::T2 # for callbacks
end

function ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,f,checkpoints,tols,tstops=nothing;noiseterm=false)
  tspan = reverse(sol.prob.tspan)
  checkpointing = ischeckpointing(sensealg, sol)
  (checkpointing && checkpoints === nothing) && error("checkpoints must be passed when checkpointing is enabled.")

  checkpoint_sol = if checkpointing
    intervals = map(tuple, @view(checkpoints[1:end-1]), @view(checkpoints[2:end]))
    interval_end = intervals[end][end]
    tspan[1] > interval_end && push!(intervals, (interval_end, tspan[1]))
    cursor = lastindex(intervals)
    interval = intervals[cursor]

    if typeof(sol.prob) <: Union{SDEProblem,RODEProblem}
      # replicated noise
      _sol = deepcopy(sol)
      idx1 = searchsortedfirst(_sol.W.t, interval[1]-1000eps(interval[1]))
      if typeof(sol.W) <: DiffEqNoiseProcess.NoiseProcess
        sol.W.save_everystep = false
        _sol.W.save_everystep = false
        forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=idx1)
      elseif typeof(sol.W) <: DiffEqNoiseProcess.NoiseGrid
        #idx2 = searchsortedfirst(_sol.W.t, interval[2]+1000eps(interval[1]))
        forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.W.t[idx1:end], _sol.W.W[idx1:end])
      else
        error("NoiseProcess type not implemented.")
      end
      dt = choose_dt((_sol.W.t[idx1]-_sol.W.t[idx1+1]), _sol.W.t, interval)

      cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1]), noise=forwardnoise),
         sol.alg, save_noise=false; dt=dt, tstops=_sol.t[idx1:end] ,tols...)
    else

      if tstops === nothing
        cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1])),sol.alg; tols...)
      else
        if maximum(interval[1] .< tstops .< interval[2])
          # callback might have changed p
          _p = reset_p(prob.kwargs[:callback], interval)
          cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1])),tstops=tstops,
                        p=_p, sol.alg; tols...)
        else
          cpsol = solve(remake(sol.prob, tspan=interval, u0=sol(interval[1])),tstops=tstops, sol.alg; tols...)
        end
      end
    end
    CheckpointSolution(cpsol, intervals, cursor, tols, tstops)

  else
    nothing
  end

  diffcache, y = adjointdiffcache(g,sensealg,discrete,sol,dg,f;quad=false,noiseterm=noiseterm)

  return ODEInterpolatingAdjointSensitivityFunction(diffcache,sensealg,
                                                    discrete,y,sol,
                                                    checkpoint_sol,sol.prob,f,noiseterm)
end

function findcursor(intervals, t)
  # equivalent with `findfirst(x->x[1] <= t <= x[2], intervals)`
  lt(x, t) = <(x[2], t)
  return searchsortedfirst(intervals, t, lt=lt)
end

function choose_dt(dt, ts, interval)
  if dt < 1000eps(interval[2])
    if length(ts) > 2
      dt = ts[end-1]-ts[end-2]
      if dt < 1000eps(interval[2])
        dt = interval[2] - interval[1]
      end
    else
      dt = interval[2] - interval[1]
    end
  end
  return dt
end

# u = λ'
# add tstop on all the checkpoints
function (S::ODEInterpolatingAdjointSensitivityFunction)(du,u,p,t)
  @unpack sol,checkpoint_sol, discrete, prob, f = S

  λ,grad,y,dλ,dgrad,dy = split_states(du,u,t,S)

  if S.noiseterm
    if length(u) == length(du)
      vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad)
    elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
      vecjacobian!(dλ, y, λ, p, t, S)
      jacNoise!(λ, y, p, t, S, dgrad=dgrad)
    else
      jacNoise!(λ, y, p, t, S, dgrad=dgrad, dλ=dλ)
    end
  else
    vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad)
  end

  dλ .*= -one(eltype(λ))

  discrete || accumulate_cost!(dλ, y, p, t, S, dgrad)
  return nothing
end

function (S::ODEInterpolatingAdjointSensitivityFunction)(du,u,p,t,W)
  @unpack sol,checkpoint_sol, discrete, prob, f = S

  λ,grad,y,dλ,dgrad,dy = split_states(du,u,t,S)

  vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad, W=W)

  dλ .*= -one(eltype(λ))

  discrete || accumulate_cost!(dλ, y, p, t, S, dgrad)
  return nothing
end

function split_states(du,u,t,S::TS;update=true) where TS<:ODEInterpolatingAdjointSensitivityFunction
  @unpack sol, y, checkpoint_sol, discrete, prob, f = S
  idx = length(y)

  if update
    if checkpoint_sol === nothing
      if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
        y = sol(t, continuity=:right)
      else
        sol(y,t, continuity=:right)
      end
    else
      intervals = checkpoint_sol.intervals
      interval = intervals[checkpoint_sol.cursor]
      if !(interval[1] <= t <= interval[2])
        cursor′ = findcursor(intervals, t)
        interval = intervals[cursor′]
        cpsol_t = checkpoint_sol.cpsol.t
        if typeof(t) <: ForwardDiff.Dual && eltype(S.y) <: AbstractFloat
          y = sol(interval[1])
        else
          sol(y, interval[1])
        end
        if typeof(sol.prob) <: Union{SDEProblem,RODEProblem}
          #idx1 = searchsortedfirst(sol.t, interval[1])
          _sol = deepcopy(sol)
          idx1 = searchsortedfirst(_sol.t, interval[1]-100eps(interval[1]))
          idx2 = searchsortedfirst(_sol.t, interval[2]+100eps(interval[2]))
          idx_noise = searchsortedfirst(_sol.W.t, interval[1]-100eps(interval[1]))
          if typeof(sol.W) <: DiffEqNoiseProcess.NoiseProcess
            _sol.W.save_everystep = false
            forwardnoise = DiffEqNoiseProcess.NoiseWrapper(_sol.W, indx=idx_noise)
          elseif typeof(sol.W) <: DiffEqNoiseProcess.NoiseGrid
            forwardnoise = DiffEqNoiseProcess.NoiseGrid(_sol.W.t[idx_noise:end], _sol.W.W[idx_noise:end])
          else
            error("NoiseProcess type not implemented.")
          end
          prob′ = remake(prob, tspan=intervals[cursor′], u0=y, noise=forwardnoise)
          dt = choose_dt(abs(cpsol_t[1]-cpsol_t[2]), cpsol_t, interval)
          cpsol′ = solve(prob′, sol.alg, noise=forwardnoise, save_noise=false; dt=dt, tstops=_sol.t[idx1:idx2], checkpoint_sol.tols...)
        else
          if checkpoint_sol.tstops===nothing
            prob′ = remake(prob, tspan=intervals[cursor′], u0=y)
            cpsol′ = solve(prob′, sol.alg; dt=abs(cpsol_t[end] - cpsol_t[end-1]), checkpoint_sol.tols...)
          else
            if maximum(interval[1] .< checkpoint_sol.tstops .< interval[2])
              # callback might have changed p
              _p = reset_p(prob.kwargs[:callback], interval)
              prob′ = remake(prob, tspan=intervals[cursor′], u0=y, p=_p)
              cpsol′ = solve(prob′, sol.alg; dt=abs(cpsol_t[end] - cpsol_t[end-1]), tstops=checkpoint_sol.tstops, checkpoint_sol.tols...)
            else
              prob′ = remake(prob, tspan=intervals[cursor′], u0=y)
              cpsol′ = solve(prob′, sol.alg; dt=abs(cpsol_t[end] - cpsol_t[end-1]), tstops=checkpoint_sol.tstops, checkpoint_sol.tols...)
            end
          end
        end
        checkpoint_sol.cpsol = cpsol′
        checkpoint_sol.cursor = cursor′
      end
      checkpoint_sol.cpsol(y, t, continuity=:right)
    end
  end

  λ     = @view u[1:idx]
  grad  = @view u[idx+1:end]

  if length(u) == length(du)
    dλ    = @view du[1:idx]
    dgrad = @view du[idx+1:end]

  elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
    idx1 = [length(u)*(i-1)+i for i in 1:idx] # for diagonal indices of [1:idx,1:idx]

    dλ    = @view du[idx1]
    dgrad = @view du[idx+1:end,1:idx]

  elseif typeof(du) <: AbstractMatrix
    # non-diagonal noise and noise mixing case
    dλ    = @view du[1:idx,1:idx]
    dgrad = @view du[idx+1:end,1:idx]
  end

  λ,grad,y,dλ,dgrad,nothing
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g::G,t=nothing,dg::DG=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     kwargs...) where {G,DG}
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  # remove duplicates from checkpoints
  if ischeckpointing(sensealg) &&  (length(unique(checkpoints)) != length(checkpoints))
    _checkpoints, duplicate_iterator_times = separate_nonunique(checkpoints)
    tstops = duplicate_iterator_times[1]
    checkpoints = filter(x->x ∉ tstops, _checkpoints)
  else
    tstops = nothing
  end

  numstates = length(u0)
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)

  len = numstates+numparams

  λ = p === nothing || p === DiffEqBase.NullParameters() ? similar(u0) : one(eltype(u0)) .* similar(p, len)
  λ .= false

  sense = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,f,
                                                     checkpoints,
                                                     (reltol=reltol,abstol=abstol),
                                                     tstops)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    adjmm = copy(sol.prob.f.mass_matrix')
    zzz = similar(adjmm, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    # using concrate I is slightly more efficient
    II = Diagonal(I, numparams)
    mm = [adjmm       zzz
          copy(zzz')   II]
  end

  jac_prototype = sol.prob.f.jac_prototype
  if !sense.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    _adjoint_jac_prototype = copy(jac_prototype')
    zzz = similar(_adjoint_jac_prototype, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    II = Diagonal(I, numparams)
    adjoint_jac_prototype = [_adjoint_jac_prototype zzz
                             copy(zzz')             II]
  end

  odefun = ODEFunction(sense, mass_matrix=mm, jac_prototype=adjoint_jac_prototype)
  return ODEProblem(odefun,z0,tspan,p,callback=cb)
end

@noinline function SDEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     diffusion_jac=nothing, diffusion_paramjac=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  # remove duplicates from checkpoints
  if ischeckpointing(sensealg) && (length(unique(checkpoints)) != length(checkpoints))
    _checkpoints, duplicate_iterator_times = separate_nonunique(checkpoints)
    tstops = duplicate_iterator_times[1]
    checkpoints = filter(x->x ∉ tstops, _checkpoints)
  else
    tstops = nothing
  end

  numstates = length(u0)
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)

  len = numstates+numparams

  λ = one(eltype(u0)) .* similar(p, len)
  λ .= false

  sense_drift = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,sol.prob.f,
                                                     checkpoints,(reltol=reltol,abstol=abstol))

  diffusion_function = ODEFunction(sol.prob.g, jac=diffusion_jac, paramjac=diffusion_paramjac)
  sense_diffusion = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,diffusion_function,
                                                     checkpoints,(reltol=reltol,abstol=abstol);noiseterm=true)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense_drift, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    adjmm = copy(sol.prob.f.mass_matrix')
    zzz = similar(adjmm, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    # using concrate I is slightly more efficient
    II = Diagonal(I, numparams)
    mm = [adjmm       zzz
          copy(zzz')   II]
  end

  jac_prototype = sol.prob.f.jac_prototype
  if !sense_drift.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    _adjoint_jac_prototype = copy(jac_prototype')
    zzz = similar(_adjoint_jac_prototype, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    II = Diagonal(I, numparams)
    adjoint_jac_prototype = [_adjoint_jac_prototype zzz
                             copy(zzz')             II]
  end

  sdefun = SDEFunction(sense_drift,sense_diffusion,mass_matrix=mm,jac_prototype=adjoint_jac_prototype)

  # replicated noise
  _sol = deepcopy(sol)
  backwardnoise = reverse(_sol.W)

  if StochasticDiffEq.is_diagonal_noise(sol.prob) && typeof(sol.W[end])<:Number
    # scalar noise case
    noise_matrix = nothing
  else
    noise_matrix = similar(z0,length(z0),numstates)
    noise_matrix .= false
  end

  return SDEProblem(sdefun,sense_diffusion,z0,tspan,p,
    callback=cb,
    noise=backwardnoise,
    noise_rate_prototype = noise_matrix
    )
end

@noinline function RODEAdjointProblem(sol,sensealg::InterpolatingAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     reltol=nothing, abstol=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  tspan = reverse(tspan)
  discrete = t != nothing

  # remove duplicates from checkpoints
  if ischeckpointing(sensealg) && (length(unique(checkpoints)) != length(checkpoints))
    _checkpoints, duplicate_iterator_times = separate_nonunique(checkpoints)
    tstops =  duplicate_iterator_times[1]
    checkpoints = filter(x->x ∉ tstops, _checkpoints)
  else
    tstops = nothing
  end

  numstates = length(u0)
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)

  len = numstates+numparams

  λ = p === nothing || p === DiffEqBase.NullParameters() ? similar(u0) : one(eltype(u0)) .* similar(p, len)
  λ .= false

  sense = ODEInterpolatingAdjointSensitivityFunction(g,sensealg,discrete,sol,dg,f,
                                                     checkpoints,
                                                     (reltol=reltol,abstol=abstol),
                                                     tstops)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense, g, λ, t, callback, init_cb)
  z0 = vec(zero(λ))
  original_mm = sol.prob.f.mass_matrix
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    adjmm = copy(sol.prob.f.mass_matrix')
    zzz = similar(adjmm, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    # using concrate I is slightly more efficient
    II = Diagonal(I, numparams)
    mm = [adjmm       zzz
          copy(zzz')   II]
  end

  jac_prototype = sol.prob.f.jac_prototype
  if !sense.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    _adjoint_jac_prototype = copy(jac_prototype')
    zzz = similar(_adjoint_jac_prototype, numstates, numparams)
    fill!(zzz, zero(eltype(zzz)))
    II = Diagonal(I, numparams)
    adjoint_jac_prototype = [_adjoint_jac_prototype zzz
                             copy(zzz')             II]
  end

  rodefun = RODEFunction(sense, mass_matrix=mm, jac_prototype=adjoint_jac_prototype)

  # replicated noise
  _sol = deepcopy(sol)
  backwardnoise = reverse(_sol.W)
  # make sure noise grid starts at correct time values, e.g., if sol.W.t is longer than sol.t
  tspan[1]!=backwardnoise.t[1] && reinit!(backwardnoise,backwardnoise.t[2]-backwardnoise.t[1],t0=tspan[1])

  if StochasticDiffEq.is_diagonal_noise(sol.prob) && typeof(sol.W[end])<:Number
    # scalar noise case
    noise_matrix = nothing
  else
    noise_matrix = similar(z0,length(z0),numstates)
    noise_matrix .= false
  end

  return RODEProblem(rodefun,z0,tspan,p,callback=cb,
                    noise=backwardnoise,
                    noise_rate_prototype = noise_matrix)
end

function reset_p(CBS, interval)
  # TODO: treat continuous callbacks
  # check which events are close to tspan[1]
  ts = map(CBS.discrete_callbacks) do cb
    indx = searchsortedfirst(cb.affect!.event_times, interval[1])
    (indx, cb.affect!.event_times[indx])
  end
  perm = minimum(sortperm([t for t in getindex.(ts,2)]))
  p = deepcopy(CBS.discrete_callbacks[perm].affect!.pleft[getindex.(ts,1)[perm]])
  return p
end