text:
# HACK: work around Pkg.jl#2500
if VERSION < v"1.8-"
test_project = Base.active_project()
preferences_file = joinpath(dirname(@__DIR__), "LocalPreferences.toml")
test_preferences_file = joinpath(dirname(test_project), "LocalPreferences.toml")
if isfile(preferences_file) && !isfile(test_preferences_file)
    cp(preferences_file, test_preferences_file)
end
end

using Enzyme
using Test
using FiniteDifferences
using ForwardDiff
using Statistics
using LinearAlgebra

using Enzyme_jll
@info "Testing against" Enzyme_jll.libEnzyme

# Test against FiniteDifferences
function test_scalar(f, x; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)
    ∂x, = autodiff(Reverse, f, Active, Active(x))
    if typeof(x) <: Complex
    else
      @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end

    rm = ∂x
    if typeof(x) <: Integer
        x = Float64(x)
    end
    ∂x, = autodiff(Forward, f, Duplicated(x, one(typeof(x))))
    if typeof(x) <: Complex
      @test ∂x ≈ rm
    else
      @test isapprox(∂x, fdm(f, x); rtol=rtol, atol=atol, kwargs...)
    end
end

include("abi.jl")
include("typetree.jl")

f0(x) = 1.0 + x
@testset "Internal tests" begin
    thunk_a = Enzyme.Compiler.thunk(f0, nothing, Active, Tuple{Active{Float64}}, Val(API.DEM_ReverseModeCombined), Val(1))
    thunk_b = Enzyme.Compiler.thunk(f0, nothing, Const, Tuple{Const{Float64}}, Val(API.DEM_ReverseModeCombined), Val(1))
    thunk_c = Enzyme.Compiler.thunk(f0, nothing, Active{Float64}, Tuple{Active{Float64}}, Val(API.DEM_ReverseModeCombined), Val(1))
    thunk_d = Enzyme.Compiler.thunk(f0, nothing, Active{Float64}, Tuple{Active{Float64}}, Val(API.DEM_ReverseModeCombined), Val(1))
    @test thunk_a.adjoint !== thunk_b.adjoint
    @test_broken thunk_c.adjoint === thunk_a.adjoint
    @test thunk_c.adjoint === thunk_d.adjoint

    @test thunk_a(Active(2.0), 1.0) == (1.0,)
    @test thunk_a(Active(2.0), 2.0) == (2.0,)
    @test thunk_b(Const(2.0)) === ()

    forward, pullback = Enzyme.Compiler.thunk(f0, nothing, Active, Tuple{Active{Float64}}, Val(Enzyme.API.DEM_ReverseModeGradient), Val(1))
    # @test thunk_split.primal !== C_NULL
    # @test thunk_split.primal !== thunk_split.adjoint
    # @test thunk_a.adjoint !== thunk_split.adjoint
end

@testset "Reflection" begin
    Enzyme.Compiler.enzyme_code_typed(Active, Tuple{Active{Float64}}) do x
        x ^ 2
    end
    sprint() do io
        Enzyme.Compiler.enzyme_code_native(io, f0, Active, Tuple{Active{Float64}})
    end

    sprint() do io
        Enzyme.Compiler.enzyme_code_llvm(io, f0, Active, Tuple{Active{Float64}})
    end
end

using  Documenter
DocMeta.setdocmeta!(Enzyme, :DocTestSetup, :(using Enzyme); recursive=true)
@testset "DocTests" begin
    doctest(Enzyme; manual = false)
end

# @testset "Split Tape" begin
#     f(x) = x[1] * x[1]

#     thunk_split = Enzyme.Compiler.thunk(f, Tuple{Duplicated{Array{Float64,1}}}, Val(Enzyme.API.DEM_ReverseModeGradient))
#     @test thunk_split.primal !== C_NULL
#     @test thunk_split.primal !== thunk_split.adjoint
# end

@testset "Simple tests" begin
    f1(x) = 1.0 + x
    f2(x) = x*x
    @test autodiff(Reverse, f1, Active, Active(1.0))[1] ≈ 1.0
    @test autodiff(Forward, f1, DuplicatedNoNeed, Duplicated(1.0, 1.0))[1] ≈ 1.0
    @test autodiff(Forward, f1, Duplicated, Duplicated(1.0, 1.0))[2] ≈ 1.0
    @test autodiff(Reverse, f2, Active, Active(1.0))[1] ≈ 2.0
    @test autodiff(Forward, f2, Duplicated(1.0, 1.0))[1] ≈ 2.0
    @test autodiff(Reverse, tanh, Active, Active(1.0))[1] ≈ 0.41997434161402606939
    @test autodiff(Forward, tanh, Duplicated(1.0, 1.0))[1] ≈ 0.41997434161402606939
    @test autodiff(Reverse, tanh, Active, Active(1.0f0))[1] ≈ Float32(0.41997434161402606939)
    @test autodiff(Forward, tanh, Duplicated(1.0f0, 1.0f0))[1] ≈ Float32(0.41997434161402606939)
    test_scalar(f1, 1.0)
    test_scalar(f2, 1.0)
    test_scalar(log2, 1.0)
    test_scalar(log1p, 1.0)

    test_scalar(log10, 1.0)

    @test autodiff(Reverse, (x)->log(x), Active(2.0)) == (0.5,)
end

@testset "Simple Exception" begin
    f_simple_exc(x, i) = ccall(:jl_, Cvoid, (Any,), x[i])
    y = [1.0, 2.0]
    f_x = zero.(y)
    @test_throws BoundsError autodiff(Reverse, f_simple_exc, Duplicated(y, f_x), 0)
end

@testset "Duplicated" begin
    x = Ref(1.0)
    y = Ref(2.0)

    ∇x = Ref(0.0)
    ∇y = Ref(0.0)

    autodiff(Reverse, (a,b)->a[]*b[], Active, Duplicated(x, ∇x), Duplicated(y, ∇y))

    @test ∇y[] == 1.0
    @test ∇x[] == 2.0
end

@testset "Simple tests" begin
    g(x) = real((x + im)*(1 - im*x))
    @test first(autodiff(Reverse, g, Active, Active(2.0))) ≈ 2.0
    @test first(autodiff(Forward, g, Duplicated(2.0, 1.0))) ≈ 2.0
    @test first(autodiff(Reverse, g, Active, Active(3.0))) ≈ 2.0
    @test first(autodiff(Forward, g, Duplicated(3.0, 1.0))) ≈ 2.0
    test_scalar(g, 2.0)
    test_scalar(g, 3.0)
end

@testset "Taylor series tests" begin

# Taylor series for `-log(1-x)`
# eval at -log(1-1/2) = -log(1/2)
function euroad(f::T) where T
    g = zero(T)
    for i in 1:10^7
        g += f^i / i
    end
    return g
end

euroad′(x) = first(autodiff(Reverse, euroad, Active, Active(x)))

@test euroad(0.5) ≈ -log(0.5) # -log(1-x)
@test euroad′(0.5) ≈ 2.0 # d/dx -log(1-x) = 1/(1-x)
test_scalar(euroad, 0.5)
end

@testset "Array tests" begin

    function arsum(f::Array{T}) where T
        g = zero(T)
        for elem in f
            g += elem
        end
        return g
    end

    inp = Float64[1.0, 2.0]
    dinp = Float64[0.0, 0.0]
    autodiff(Reverse, arsum, Active, Duplicated(inp, dinp))
    @test inp ≈ Float64[1.0, 2.0]
    @test dinp ≈ Float64[1.0, 1.0]

    @test autodiff(Forward, arsum, Duplicated(inp, dinp))[1] ≈ 2.0
end

@testset "Advanced array tests" begin
    function arsum2(f::Array{T}) where T
        return sum(f)
    end
    inp = Float64[1.0, 2.0]
    dinp = Float64[0.0, 0.0]
    autodiff(Reverse, arsum2, Active, Duplicated(inp, dinp))
    @test inp ≈ Float64[1.0, 2.0]
    @test dinp ≈ Float64[1.0, 1.0]

    @test autodiff(Forward, arsum2, Duplicated(inp, dinp))[1] ≈ 2.0
end

@testset "Dict" begin
    params = Dict{Symbol, Float64}()
    dparams = Dict{Symbol, Float64}()

    params[:var] = 10.0
    dparams[:var] = 0.0

    f_dict(params, x) = params[:var] * x

    @test autodiff(Reverse, f_dict, Const(params), Active(5.0)) == (10.0,)
    @test autodiff(Reverse, f_dict, Duplicated(params, dparams), Active(5.0)) == (10.0,)
    @test dparams[:var] == 5.0

    mutable struct MD
        v::Float64
        d::Dict{Symbol, MD}
    end

    # TODO without Float64 on return
    # there is a potential phi bug
    function sum_rec(d::Dict{Symbol,MD})::Float64
        s = 0.0
        for k in keys(d)
            s += d[k].v
            s += sum_rec(d[k].d)
        end
        return s
    end

    par = Dict{Symbol, MD}()
    par[:var] = MD(10.0, Dict{Symbol, MD}())
    par[:sub] = MD(2.0, Dict{Symbol, MD}(:a=>MD(3.0, Dict{Symbol, MD}())))

    dpar = Dict{Symbol, MD}()
    dpar[:var] = MD(0.0, Dict{Symbol, MD}())
    dpar[:sub] = MD(0.0, Dict{Symbol, MD}(:a=>MD(0.0, Dict{Symbol, MD}())))

    # TODO
    # autodiff(Reverse, sum_rec, Duplicated(par, dpar))
    # @show par, dpar, sum_rec(par)
    # @test dpar[:var].v ≈ 1.0
    # @test dpar[:sub].v ≈ 1.0
    # @test dpar[:sub].d[:a].v ≈ 1.0
end

function grad_closure(f, x)
    function noretval(x,res)
        y = f(x)
        copyto!(res,y)
        return nothing
    end
    n = length(x)
    dx = zeros(n)
    y  = zeros(n)
    dy = zeros(n)
    dy[1] = 1.0

    autodiff(Reverse, noretval, Duplicated(x,dx), Duplicated(y, dy))
    return dx
end

@testset "Closure" begin
    x = [2.0,6.0]
    dx = grad_closure(x->[x[1], x[2]], x)
    @test dx == [1.0, 0.0]
end

@testset "Bithacks" begin
    function fneg(x::Float64)
        xptr = reinterpret(Int64, x)
        y = Int64(-9223372036854775808)
        out = y ⊻ xptr;
        return reinterpret(Float64, out)
    end
    @test autodiff(Reverse, fneg, Active, Active(2.0))[1] ≈ -1.0
    @test autodiff(Forward, fneg, Duplicated(2.0, 1.0))[1] ≈ -1.0
    function expor(x::Float64)
        xptr = reinterpret(Int64, x)
        y = UInt64(4607182418800017408)
        out = y | xptr;
        return reinterpret(Float64, out)
    end
    @test autodiff(Reverse, expor, Active, Active(0.42))[1] ≈ 4.0
    @test autodiff(Forward, expor, Duplicated(0.42, 1.0))[1] ≈ 4.0
end

@testset "GC" begin
    function gc_alloc(x)  # Basically g(x) = x^2
        a = Array{Float64, 1}(undef, 10)
        for n in 1:length(a)
            a[n] = x^2
        end
        return mean(a)
    end
    @test autodiff(Reverse, gc_alloc, Active, Active(5.0))[1] ≈ 10
    @test autodiff(Forward, gc_alloc, Duplicated(5.0, 1.0))[1] ≈ 10

    A = Float64[2.0, 3.0]
    B = Float64[4.0, 5.0]
    dB = Float64[0.0, 0.0]
    f = (X, Y) -> sum(X .* Y)
    Enzyme.autodiff(Reverse, f, Active, A, Duplicated(B, dB))

    function gc_copy(x)  # Basically g(x) = x^2
        a = x * ones(10)
        for n in 1:length(a)
            a[n] = x^2
        end
        return mean(a)
    end

    @test Enzyme.autodiff(Reverse, gc_copy, Active, Active(5.0))[1] ≈ 10
    @test Enzyme.autodiff(Forward, gc_copy, Duplicated(5.0, 1.0))[1] ≈ 10
end

@testset "Compare against" begin
    x = 3.0
    fd = central_fdm(5, 1)(sin, x)

    @test fd ≈ ForwardDiff.derivative(sin, x)
    @test fd ≈ first(autodiff(Reverse, sin, Active, Active(x)))
    @test fd ≈ first(autodiff(Forward, sin, Duplicated(x, 1.0)))

    x = 0.2 + sin(3.0)
    fd = central_fdm(5, 1)(asin, x)

    @test fd ≈ ForwardDiff.derivative(asin, x)
    @test fd ≈ first(autodiff(Reverse, asin, Active, Active(x)))
    @test fd ≈ first(autodiff(Forward, asin, Duplicated(x, 1.0)))
    test_scalar(asin, x)

    function foo(x)
        a = sin(x)
        b = 0.2 + a
        c = asin(b)
        return c
    end

    x = 3.0
    fd = central_fdm(5, 1)(foo, x)

    @test fd ≈ ForwardDiff.derivative(foo, x)
    @test fd ≈ first(autodiff(Reverse, foo, Active, Active(x)))
    @test fd ≈ first(autodiff(Forward, foo, Duplicated(x, 1.0)))
    test_scalar(foo, x)

    # Input type shouldn't matter
    x = 3
    @test fd ≈ ForwardDiff.derivative(foo, x)
    @test fd ≈ first(autodiff(Reverse, foo, Active, Active(x)))
    # They do matter for duplicated, which can't be auto promoted
    # @test fd ≈ first(autodiff(Forward, foo, Duplicated(x, 1)))

    f74(a, c) = a * √c
    @test √3 ≈ first(autodiff(Reverse, f74, Active, Active(2), 3))
    @test √3 ≈ first(autodiff(Forward, f74, Duplicated(2.0, 1.0), 3))
end

@testset "SinCos" begin
	function sumsincos(theta)
		a, b = sincos(theta)
		return a + b
	end
    test_scalar(sumsincos, 1.0, rtol=1e-5, atol=1e-5)
end

"""
    J(ν, z) := ∑ (−1)^k / Γ(k+1) / Γ(k+ν+1) * (z/2)^(ν+2k)
"""
function mybesselj(ν, z, atol=1e-8)
    k = 0
    s = (z/2)^ν / factorial(ν)
    out = s
    while abs(s) > atol
        k += 1
        s *= (-1) / k / (k+ν) * (z/2)^2
        out += s
    end
    out
end
mybesselj0(z) = mybesselj(0, z)
mybesselj1(z) = mybesselj(1, z)

@testset "Bessel" begin
    autodiff(Reverse, mybesselj, Active, Const(0), Active(1.0))
    autodiff(Reverse, mybesselj, Active, 0, Active(1.0))
    autodiff(Forward, mybesselj, Const(0), Duplicated(1.0, 1.0))
    autodiff(Forward, mybesselj, 0, Duplicated(1.0, 1.0))
    @testset "besselj0/besselj1" for x in (1.0, -1.0, 0.0, 0.5, 10, -17.1,) # 1.5 + 0.7im)
        test_scalar(mybesselj0, x, rtol=1e-5, atol=1e-5)
        test_scalar(mybesselj1, x, rtol=1e-5, atol=1e-5)
    end
end

## https://github.com/JuliaDiff/ChainRules.jl/tree/master/test/rulesets
if !Sys.iswindows()
    include("packages/specialfunctions.jl")
end

@testset "Threads" begin
    cmd = `$(Base.julia_cmd()) --threads=1 --startup-file=no threads.jl`
   	@test success(pipeline(cmd, stderr=stderr, stdout=stdout))
    cmd = `$(Base.julia_cmd()) --threads=2 --startup-file=no threads.jl`
   	@test success(pipeline(cmd, stderr=stderr, stdout=stdout))
end

@testset "DiffTest" begin
    include("DiffTests.jl")

    n = rand()
    x, y = rand(5, 5), rand(26)
    A, B = rand(5, 5), rand(5, 5)

    # f returns Number
    @testset "Number to Number" for f in DiffTests.NUMBER_TO_NUMBER_FUNCS
        test_scalar(f, n)
    end

    # TODO(vchuravy/wsmoses): Enable these tests
    # for f in DiffTests.VECTOR_TO_NUMBER_FUNCS
    #     @test isa(f(y), Number)
    # end

    # for f in DiffTests.MATRIX_TO_NUMBER_FUNCS
    #     @test isa(f(x), Number)
    # end

    # for f in DiffTests.TERNARY_MATRIX_TO_NUMBER_FUNCS
    #     @test isa(f(A, B, x), Number)
    # end

    # # f returns Array

    # for f in DiffTests.NUMBER_TO_ARRAY_FUNCS
    #     @test isa(f(n), Array)
    # end

    # for f in DiffTests.ARRAY_TO_ARRAY_FUNCS
    #     @test isa(f(A), Array)
    #     @test isa(f(y), Array)
    # end

    # for f in DiffTests.MATRIX_TO_MATRIX_FUNCS
    #     @test isa(f(A), Array)
    # end

    # for f in DiffTests.BINARY_MATRIX_TO_MATRIX_FUNCS
    #     @test isa(f(A, B), Array)
    # end

    # # f! returns Nothing

    # for f! in DiffTests.INPLACE_ARRAY_TO_ARRAY_FUNCS
    #     @test isa(f!(y, x), Nothing)
    # end

    # for f! in DiffTests.INPLACE_NUMBER_TO_ARRAY_FUNCS
    #     @test isa(f!(y, n), Nothing)
    # end

end

@testset "IO" begin

    function printsq(x)
        println(x)
        x*x
    end

    @test 4.6 ≈ first(autodiff(Reverse, printsq, Active, Active(2.3)))
    @test 4.6 ≈ first(autodiff(Forward, printsq, Duplicated(2.3, 1.0)))

    function tostring(x)
        string(x)
        x*x
    end

    @test 4.6 ≈ first(autodiff(Reverse, tostring, Active, Active(2.3)))
    @test 4.6 ≈ first(autodiff(Forward, tostring, Duplicated(2.3, 1.0)))
end

@testset "hmlstm" begin
    sigm(x)  = @fastmath 1 / (1 + exp(-x))
    @fastmath function hmlstm_update_c_scalar(z, zb, c, f, i, g)
        if z == 1.0f0 # FLUSH
            return sigm(i) * tanh(g)
        elseif zb == 0.0f0 # COPY
            return c
        else # UPDATE
            return sigm(f) * c + sigm(i) * tanh(g)
        end
    end

    N = 64
    Z = round.(rand(Float32, N))
    Zb = round.(rand(Float32, N))
    C = rand(Float32, N, N)
    F = rand(Float32, N, N)
    I = rand(Float32, N, N)
    G = rand(Float32, N, N)

    function broadcast_hmlstm(out, Z, Zb, C, F, I, G)
        out .= hmlstm_update_c_scalar.(Z, Zb, C, F, I, G)
        return nothing
    end

    ∇C = zeros(Float32, N, N)
    ∇F = zeros(Float32, N, N)
    ∇I = zeros(Float32, N, N)
    ∇G = zeros(Float32, N, N)

    # TODO(wsmoses): Check after updating Enzyme_jll
    # autodiff(broadcast_hmlstm, Const,
    #          Const(zeros(Float32, N, N)), Const(Z), Const(Zb),
    #          Duplicated(C, ∇C), Duplicated(F, ∇F), Duplicated(I, ∇I), Duplicated(G, ∇G))
    # fwddiff(broadcast_hmlstm, Const,
    #          Const(zeros(Float32, N, N)), Const(Z), Const(Zb),
    #          Duplicated(C, ∇C), Duplicated(F, ∇F), Duplicated(I, ∇I), Duplicated(G, ∇G))
end

genlatestsin(x)::Float64 = Base.invokelatest(sin, x)
function genlatestsinx(xp)
    x = @inbounds xp[1]
    @inbounds xp[1] = 0.0
    Base.invokelatest(sin, x)::Float64 + 1
end

function loadsin(xp)
    x = @inbounds xp[1]
    @inbounds xp[1] = 0.0
    sin(x)
end
function invsin(xp)
    xp = Base.invokelatest(convert, Vector{Float64}, xp)
    loadsin(xp)
end

@testset "generic" begin
    @test -0.4161468365471424 ≈ Enzyme.autodiff(Reverse, genlatestsin, Active, Active(2.0))[1]
    @test -0.4161468365471424 ≈ Enzyme.autodiff(Forward, genlatestsin, Duplicated(2.0, 1.0))[1]

    x = [2.0]
    dx = [0.0]
    Enzyme.autodiff(Reverse, genlatestsinx, Active, Duplicated(x, dx))
    @test 0 ≈ x[1]
    @test -0.4161468365471424 ≈ dx[1]

    x = [2.0]
    dx = [0.0]
    Enzyme.autodiff(Reverse, invsin, Active, Duplicated(x, dx))
    @test 0 ≈ x[1]
    @test -0.4161468365471424 ≈ dx[1]
end

@testset "invoke" begin
    @noinline apply(@nospecialize(func)) = func()

    function invtest(arr)
        function f()
           arr[1] *= 5.0
           nothing
        end
        apply(f)
    end

    x  = [2.0]
    dx = [1.0]

    Enzyme.autodiff(Reverse, invtest, Duplicated(x, dx))

    @test 10.0 ≈ x[1]
    @test 5.0 ≈ dx[1]
end

@testset "broadcast" begin
    A = rand(10); B = rand(10); R = similar(A)
    dA = zero(A); dB = zero(B); dR = fill!(similar(R), 1)

    function foo_bc!(R, A, B)
        R .= A .+ B
        return nothing
    end

    autodiff(Reverse, foo_bc!, Const, Duplicated(R, dR), Duplicated(A, dA), Duplicated(B, dB))

    # works since aliasing is "simple"
    autodiff(Reverse, foo_bc!, Const, Duplicated(R, dR), Duplicated(R, dR), Duplicated(B, dB))

    A = rand(10,10); B = rand(10, 10)
    dA = zero(A); dB = zero(B); dR = fill!(similar(A), 1)

    autodiff(Reverse, foo_bc!, Const, Duplicated(A, dR), Duplicated(transpose(A), transpose(dA)), Duplicated(B, dB))
end

@testset "DuplicatedReturn" begin
    moo(x) = fill(x, 10)

    @test_throws ErrorException autodiff(Reverse, moo, Active(2.1))
    fo, = autodiff(Forward, moo, Duplicated(2.1, 1.0))
    for i in 1:10
        @test 1.0 ≈ fo[i]
    end
end

@testset "GCPreserve" begin
    function f(x, y)
        GC.@preserve x y begin
            ccall(:memcpy, Cvoid,
                (Ptr{Float64},Ptr{Float64},Csize_t), x, y, 8)
        end
        nothing
    end
    autodiff(Reverse, f, Duplicated([1.0], [0.0]), Duplicated([1.0], [0.0]))
    autodiff(Forward, f, Duplicated([1.0], [0.0]), Duplicated([1.0], [0.0]))
end

@testset "GCPreserve2" begin
    function f!(a_out, a_in)
           a_out[1:end-1] .= a_in[2:end]
           return nothing
    end
    a_in = rand(4)
    a_out = a_in

    shadow_a_out = ones(4)
    shadow_a_in = shadow_a_out

    autodiff(Reverse, f!, Const, Duplicated(a_out, shadow_a_out), Duplicated(a_in, shadow_a_in))

    @test shadow_a_in ≈ Float64[0.0, 1.0, 1.0, 2.0]
    @test shadow_a_out ≈ Float64[0.0, 1.0, 1.0, 2.0]

    autodiff(Forward, f!, Const, Duplicated(a_out, shadow_a_out), Duplicated(a_in, shadow_a_in))

    @test shadow_a_in ≈ Float64[1.0, 1.0, 2.0, 2.0]
    @test shadow_a_out ≈ Float64[1.0, 1.0, 2.0, 2.0]
end

@testset "UndefVar" begin
    function f_undef(x, y)
        if x
            undefinedfnthowmagic()
        end
        y
    end
    @test 1.0 ≈ autodiff(Reverse, f_undef, false, Active(2.14))[1]
    @test_throws Base.UndefVarError autodiff(Reverse, f_undef, true, Active(2.14))

    @test 1.0 ≈ autodiff(Forward, f_undef, false, Duplicated(2.14, 1.0))[1]
    @test_throws Base.UndefVarError autodiff(Forward, f_undef, true, Duplicated(2.14, 1.0))
end

@testset "Return GC error" begin
	t = 0.0

	function tobedifferentiated(cond, a)::Float64
		if cond
			t + t
		else
			0.0
		end
	end

	@test 0.0 ≈ autodiff(Reverse, tobedifferentiated, true, Active(2.1))[1]
	@test 0.0 ≈ autodiff(Forward, tobedifferentiated, true, Duplicated(2.1, 1.0))[1]

	function tobedifferentiated2(cond, a)::Float64
		if cond
			a + t
		else
			0.0
		end
	end

	@test 1.0 ≈ autodiff(Reverse, tobedifferentiated2, true, Active(2.1))[1]
	@test 1.0 ≈ autodiff(Forward, tobedifferentiated2, true, Duplicated(2.1, 1.0))[1]

    @noinline function copy(dest, p1, cond)
        bc = convert(Broadcast.Broadcasted{Nothing}, Broadcast.instantiate(p1))

        if cond
            return nothing
        end

        bc2 = Broadcast.preprocess(dest, bc)
        @inbounds    dest[1] = bc2[1]

        nothing
    end

    function mer(F, F_H, cond)
        p1 = Base.broadcasted(Base.identity, F_H)
        copy(F, p1, cond)

        # Force an apply generic
        flush(stdout)
        nothing
    end

    L_H = Array{Float64, 1}(undef, 2)
    L = Array{Float64, 1}(undef, 2)

    F_H = [1.0, 0.0]
    F = [1.0, 0.0]

    autodiff(Reverse, mer, Duplicated(F, L), Duplicated(F_H, L_H), true)
    autodiff(Forward, mer, Duplicated(F, L), Duplicated(F_H, L_H), true)
end

@testset "Split GC" begin
    @noinline function bmat(x)
        data = [x]
        return data
    end

    function f(x::Float64)
        @inbounds return bmat(x)[1]
    end
    @test 1.0 ≈ autodiff(Reverse, f, Active(0.1))[1]
    @test 1.0 ≈ autodiff(Forward, f, Duplicated(0.1, 1.0))[1]
end

@testset "Array Copy" begin
	F = [2.0, 3.0]

	dF = [0.0, 0.0]

	function copytest(F)
		F2 = copy(F)
		@inbounds F[1] = 1.234
		@inbounds F[2] = 5.678
		@inbounds F2[1] * F2[2]
	end
	autodiff(Reverse, copytest, Duplicated(F, dF))
	@test F ≈ [1.234, 5.678]
	@test dF ≈ [3.0, 2.0]

    @test 31.0 ≈ autodiff(Forward, copytest, Duplicated([2.0, 3.0], [7.0, 5.0]))[1]
end

@testset "No inference" begin
    c = 5.0
    @test 5.0 ≈ autodiff(Reverse, (A,)->c * A, Active, Active(2.0))[1]
    @test 5.0 ≈ autodiff(Forward, (A,)->c * A, Duplicated(2.0, 1.0))[1]
end

@testset "Type-instable capture" begin
    L = Array{Float64, 1}(undef, 2)

    F = [1.0, 0.0]

    function main()
        t = 0.0

        function cap(m)
            t = m
        end

        @noinline function inner(F, cond)
            if cond
                genericcall(F)
            end
        end

        function tobedifferentiated(F, cond)
            inner(F, cond)
            # Force an apply generic
            -t
            nothing
        end
        autodiff(Reverse, tobedifferentiated, Duplicated(F, L), false)
        autodiff(Forward, tobedifferentiated, Duplicated(F, L), false)
    end

    main()
end

@testset "Arrays are double pointers" begin
    @noinline function func_scalar(X)
        return X
    end

    function timsteploop_scalar(FH1)
        G = Float64[FH1]
        k1 = @inbounds func_scalar(G[1])
        return k1
    end
    @test Enzyme.autodiff(Reverse, timsteploop_scalar, Active(2.0))[1] ≈ 1.0
    @test Enzyme.autodiff(Forward, timsteploop_scalar, Duplicated(2.0, 1.0))[1] ≈ 1.0

    @noinline function func(X)
        return @inbounds X[1]
    end
    function timsteploop(FH1)
        G = Float64[FH1]
        k1 = func(G)
        return k1
    end
    @test Enzyme.autodiff(Reverse, timsteploop, Active(2.0))[1] ≈ 1.0
    @test Enzyme.autodiff(Forward, timsteploop, Duplicated(2.0, 1.0))[1] ≈ 1.0
end

@testset "Type" begin
    function foo(in::Ptr{Cvoid}, out::Ptr{Cvoid})
        markType(Float64, in)
        ccall(:memcpy,Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), out, in, 8)
    end

    x = [2.0]
    y = [3.0]
    dx = [5.0]
    dy = [7.0]

    @test markType(x) === nothing
    @test markType(zeros(Float32, 64)) === nothing
    @test markType(view(zeros(64), 16:32)) === nothing

    GC.@preserve x y begin
        foo(Base.unsafe_convert(Ptr{Cvoid}, x), Base.unsafe_convert(Ptr{Cvoid}, y))
    end

    GC.@preserve x y dx dy begin
      autodiff(foo,
                Duplicated(Base.unsafe_convert(Ptr{Cvoid}, x), Base.unsafe_convert(Ptr{Cvoid}, dx)),
                Duplicated(Base.unsafe_convert(Ptr{Cvoid}, y), Base.unsafe_convert(Ptr{Cvoid}, dy)))
    end
end

@testset "BLAS" begin
    x = [2.0, 3.0]
    dx = [0.2,0.3]
    y = [5.0, 7.0]
    dy = [0.5,0.7]
    Enzyme.autodiff(Reverse, (x,y)->x' * y, Duplicated(x, dx), Duplicated(y, dy))
    @show x, dx, y, dy
    @test dx ≈ [5.2, 7.3]
    @test dy ≈ [2.5, 3.7]

    f_exc(x) = sum(x*x)
    y = [[1.0, 2.0] [3.0,4.0]]
    f_x = zero.(y)
    Enzyme.autodiff(Reverse, f_exc, Duplicated(y, f_x))
    @test f_x ≈ [7.0 9.0; 11.0 13.0]
end

@testset "Exception" begin

    f_no_derv(x) = ccall("extern doesnotexist", llvmcall, Float64, (Float64,), x)
    @test_throws Enzyme.Compiler.NoDerivativeException autodiff(Reverse, f_no_derv, Active, Active(0.5))

    f_union(cond, x) = cond ? x : 0
    g_union(cond, x) = f_union(cond,x)*x
    @test_throws Enzyme.Compiler.IllegalTypeAnalysisException autodiff(Reverse, g_union, Active, true, Active(1.0))

    # TODO: Add test for NoShadowException
end

@testset "Array push" begin

    function pusher(x, y)
        push!(x, y)
        x[1] + x[2]
    end

    x  = [2.3]
    dx = [0.0]
    @test 1.0 ≈ first(Enzyme.autodiff(Reverse, pusher, Duplicated(x, dx), Active(2.0)))
    @test x ≈ [2.3, 2.0]
    @test dx ≈ [1.0]
end

@testset "Batch Forward" begin
    square(x)=x*x
    bres = autodiff(Forward, square, BatchDuplicatedNoNeed, BatchDuplicated(3.0, (1.0, 2.0, 3.0)))
    @test length(bres) == 1
    @test length(bres[1]) == 3
    @test bres[1][1] ≈  6.0
    @test bres[1][2] ≈ 12.0
    @test bres[1][3] ≈ 18.0

    bres = autodiff(Forward, square, BatchDuplicatedNoNeed, BatchDuplicated(3.0 + 7.0im, (1.0+0im, 2.0+0im, 3.0+0im)))
    @test bres[1][1] ≈  6.0 + 14.0im
    @test bres[1][2] ≈ 12.0 + 28.0im
    @test bres[1][3] ≈ 18.0 + 42.0im

    squareidx(x)=x[1]*x[1]
    inp = Float32[3.0]

    # Shadow offset is not the same as primal so following doesn't work
    # d_inp = Float32[1.0, 2.0, 3.0]
    # autodiff(Forward, squareidx, BatchDuplicatedNoNeed, BatchDuplicated(view(inp, 1:1), (view(d_inp, 1:1), view(d_inp, 2:2), view(d_inp, 3:3))))

    d_inp = (Float32[1.0], Float32[2.0], Float32[3.0])
    bres = autodiff(Forward, squareidx, BatchDuplicatedNoNeed, BatchDuplicated(inp, d_inp))
    @test bres[1][1] ≈  6.0
    @test bres[1][2] ≈ 12.0
    @test bres[1][3] ≈ 18.0
end

@testset "Batch Reverse" begin
    function refbatchbwd(out, x)
        v = x[]
        out[1] = v
        out[2] = v*v
        out[3] = v*v*v
        nothing
    end

    dxs = (Ref(0.0), Ref(0.0), Ref(0.0))
    out = Float64[0,0,0]
    x = Ref(2.0)

    autodiff(Reverse, refbatchbwd, BatchDuplicated(out, Enzyme.onehot(out)), BatchDuplicated(x, dxs))
    @test dxs[1][] ≈  1.0
    @test dxs[2][] ≈  4.0
    @test dxs[3][] ≈ 12.0

    function batchbwd(out, v)
        out[1] = v
        out[2] = v*v
        out[3] = v*v*v
        nothing
    end

    bres = Enzyme.autodiff(Reverse, batchbwd, BatchDuplicated(out, Enzyme.onehot(out)), Active(2.0))
    @test length(bres) == 1
    @test length(bres[1]) == 3
    @test bres[1][1] ≈  1.0
    @test bres[1][2] ≈  4.0
    @test bres[1][3] ≈ 12.0
end

@testset "Jacobian" begin
    function inout(v)
       [v[2], v[1]*v[1], v[1]*v[1]*v[1]]
    end

    jac = Enzyme.jacobian(Reverse, inout, [2.0, 3.0], #=n_outs=# Val(3), Val(1))
    @test size(jac) == (3, 2)
    @test jac ≈ [ 0.0   1.0;
                  4.0   0.0;
                  12.0  0.0]

    jac = Enzyme.jacobian(Forward, inout, [2.0, 3.0], Val(1))
    @test size(jac) == (3, 2)
    @test jac ≈ [ 0.0   1.0;
                  4.0   0.0;
                  12.0  0.0]

    @test jac == Enzyme.jacobian(Forward, inout, [2.0, 3.0])
    @test jac == ForwardDiff.jacobian(inout, [2.0, 3.0])

    jac = Enzyme.jacobian(Reverse, inout, [2.0, 3.0], #=n_outs=# Val(3), Val(2))
    @test size(jac) == (3, 2)
    @test jac ≈ [ 0.0   1.0;
                  4.0   0.0;
                  12.0  0.0]

    jac = Enzyme.jacobian(Forward, inout, [2.0, 3.0], Val(2))
    @test size(jac) == (3, 2)
    @test jac ≈ [ 0.0   1.0;
                  4.0   0.0;
                  12.0  0.0]

    function f_test_1(A, x)
        u = A*x[2:end] .+ x[1]
        return u
    end

    function f_test_2(A, x)
        u = Vector{Float64}(undef, length(x)-1)
        u .= A*x[2:end] .+ x[1]
        return u
    end

    function f_test_3!(u, A, x)
        u .= A*x[2:end] .+ x[1]
    end

    J_r_1(A, x) = Enzyme.jacobian(Reverse, θ -> f_test_1(A, θ), x, Val(5))
    J_r_2(A, x) = Enzyme.jacobian(Reverse, θ -> f_test_2(A, θ), x, Val(5))
    J_r_3(u, A, x) = Enzyme.jacobian(Reverse, θ -> f_test_3!(u, A, θ), x, Val(5))

    J_f_1(A, x) = Enzyme.jacobian(Forward, θ -> f_test_1(A, θ), x)
    J_f_2(A, x) = Enzyme.jacobian(Forward, θ -> f_test_2(A, θ), x)
    J_f_3(u, A, x) = Enzyme.jacobian(Forward, θ -> f_test_3!(u, A, θ), x)

    x = ones(6)
    A = Matrix{Float64}(LinearAlgebra.I, 5, 5)
    u = Vector{Float64}(undef, 5)

    @test J_r_1(A, x) == [
        1.0  1.0  0.0  0.0  0.0  0.0;
        1.0  0.0  1.0  0.0  0.0  0.0;
        1.0  0.0  0.0  1.0  0.0  0.0;
        1.0  0.0  0.0  0.0  1.0  0.0;
        1.0  0.0  0.0  0.0  0.0  1.0;
    ]

    @test_broken J_r_2(A, x) == [
        1.0  1.0  0.0  0.0  0.0  0.0;
        1.0  0.0  1.0  0.0  0.0  0.0;
        1.0  0.0  0.0  1.0  0.0  0.0;
        1.0  0.0  0.0  0.0  1.0  0.0;
        1.0  0.0  0.0  0.0  0.0  1.0;
    ]

    # Function fails verification in test/CI
    # @test J_f_1(A, x) == [
    #     1.0  1.0  0.0  0.0  0.0  0.0;
    #     1.0  0.0  1.0  0.0  0.0  0.0;
    #     1.0  0.0  0.0  1.0  0.0  0.0;
    #     1.0  0.0  0.0  0.0  1.0  0.0;
    #     1.0  0.0  0.0  0.0  0.0  1.0;
    # ]
    # @test J_f_2(A, x) == [
    #     1.0  1.0  0.0  0.0  0.0  0.0;
    #     1.0  0.0  1.0  0.0  0.0  0.0;
    #     1.0  0.0  0.0  1.0  0.0  0.0;
    #     1.0  0.0  0.0  0.0  1.0  0.0;
    #     1.0  0.0  0.0  0.0  0.0  1.0;
    # ]

    # Bug on (augmented) forward pass deducing if
	# shadow value is used
    # @show J_r_3(u, A, x)
    # @show J_f_3(u, A, x)
end

@testset "Forward on Reverse" begin

	function speelpenning(y, x)
		ccall(:memmove, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
								  y, x, 2 * 8)
		return nothing
	end

	x = [0.5, 0.3]
	y = zeros(2)
    dx = ones(2)
    rx = zeros(2)
    drx = zeros(2)
    dy = zeros(2)
    ry = ones(2)
    dry = zeros(2)

    function foo(y, dy, x, dx)
        Enzyme.autodiff_deferred(speelpenning, Const, Duplicated(y, dy), Duplicated(x, dx))
        return nothing
    end

    autodiff(Forward, foo, Duplicated(x, dx), Duplicated(rx, drx), Duplicated(y, dy), Duplicated(ry, dry))
end

using CUDA
if CUDA.functional() && VERSION >= v"1.7.0"
    include("cuda.jl")
end

using Random

@testset "Random" begin
	f_rand(x) = x*rand()
	f_randn(x, N) = x*sum(randn(N))
    autodiff(f_rand, Active, Active(1.0))
    autodiff(f_randn, Active, Active(1.0), Const(64))
end


text:
using Statistics, StatsBase

function MNmodel(θ, rndseed=1234)
    n = 1000
    μ1, μ2, σ1, σ2, prob = θ
    d1=randn(n).*σ1 .+ μ1
    d2=randn(n).*(σ1+σ2) .+ (μ1 - μ2) # second component lower mean and higher variance
    ps=rand(n).<prob
    data=zeros(n)
    data[ps].=d1[ps]
    data[.!ps].=d2[.!ps]
    return data
end

function auxstat(data)
    r = 0.0 : 0.1 : 1.0
    sqrt(1000.).*vcat(mean(data), std(data), skewness(data), kurtosis(data),
        quantile.(Ref(data),r))
end

function auxstat(θ, reps)
    auxstat.([MNmodel(θ, rand(1:Int64(1e12))) for i = 1:reps])
end

function TrueParameters()
    [1.0, 1.0, 0.2, 1.8, 0.4] # first component N(1,0.2) second component N(0,2)
end

function PriorSupport()
    lb = [0.0, 0.0, 0.0, 0.0, 0.05] # there is always at least 5% prob for each component
    ub = [3.0, 3.0, 1.0, 3.0, 0.95]
    lb,ub
end

function PriorDraw()
    lb, ub = PriorSupport()
    θ = (ub-lb).*rand(size(lb,1)) + lb
end

function InSupport(θ)
    lb,ub = PriorSupport()
    all(θ .>= lb) & all(θ .<= ub)
end

# prior should be an array of distributions, one for each parameter
lb, ub = PriorSupport() # need these in Prior
macro Prior()
    return :( arraydist([Uniform(lb[i], ub[i]) for i = 1:size(lb,1)]) )
end

ext: jl

max_stars_repo_name: mcreel/SimulatedNeuralMoments

max_issues_repo_name: mcreel/SimulatedNeuralMoments

max_forks_repo_name: mcreel/SimulatedNeuralMoments


text:
# Turing.jl version of model at https://github.com/stan-dev/blob/master/misc/cluster/lda/lda.stan

# data {
#   int<lower=2> K;               // num topics
#   int<lower=2> V;               // num words
#   int<lower=1> M;               // num docs
#   int<lower=1> N;               // total word instances
#   int<lower=1,upper=V> w[N];    // word n
#   int<lower=1,upper=M> doc[N];  // doc ID for word n
#   vector<lower=0>[K] alpha;     // topic prior
#   vector<lower=0>[V] beta;      // word prior
# }
# parameters {
#   simplex[K] theta[M];   // topic dist for doc m
#   simplex[V] phi[K];     // word dist for topic k
# }
# model {
#   for (m in 1:M)
#     theta[m] ~ dirichlet(alpha);  // prior
#   for (k in 1:K)
#     phi[k] ~ dirichlet(beta);     // prior
#   for (n in 1:N) {
#     real gamma[K];
#     for (k in 1:K)
#       gamma[k] <- log(theta[doc[n],k]) + log(phi[k,w[n]]);
#     increment_log_prob(log_sum_exp(gamma));  // likelihood
#   }
# }

# V = ldastandata[1]["V"]
# K = ldastandata[1]["K"]
# M = ldastandata[1]["M"]
# theta = Vector{Vector{Real}}(undef, M)
# for m = 1:M
#   theta[m] = Vector{Real}(undef, K)
# end
# phi = Vector{Vector{Real}}(undef, K)
# for k = 1:K
#   phi[k] = Vector{Real}(undef, V)
# end

# phi_dot_theta = [Vector{Real}(undef, V) for m=1:M]

@model ldamodel(K, V, M, N, w, doc, beta, alpha) = begin
  theta = TArray{Vector{Real}, 1}(M)
  @simd for m = 1:M
    @inbounds theta[m] ~ Dirichlet(alpha)
  end

  phi = TArray{Vector{Real}, 1}(K)
  @simd for k = 1:K
    @inbounds phi[k] ~ Dirichlet(beta)
  end

  # z = tzeros(Int, N)
  # for n = 1:N
  #   z[n] ~ Categorical(theta[doc[n]])
  # end

  # phi_dot_theta = [log.([dot(map(p -> p[v], phi), theta[m]) for v = 1:V]) for m = 1:M]

  # @simd for m=1:M
  #   @inbounds phi_dot_theta[m] = log.([dot(map(p -> p[v], phi), theta[m]) for v = 1:V])
  # end

  #for n = 1:N
  #  # phi_dot_theta = [dot(map(p -> p[v], phi), theta[doc[n]]) for v = 1:V]
  #  # w[n] ~ Categorical(phi_dot_theta)
  #  Turing.acclogp!(vi, phi_dot_theta[doc[n]][w[n]])
  #end

  # lp = mapreduce(n->phi_dot_theta[doc[n]][w[n]], +, 1:N)

  # lp = phi_dot_theta[doc[1]][w[1]]
  # @simd for n = 2:N
  #   @inbounds lp += phi_dot_theta[doc[n]][w[n]]
  # end

  phi_dot_theta = log.(hcat(phi...) * hcat(theta...))
  _lp += mapreduce(n->phi_dot_theta[w[n], doc[n]], +, 1:N)

  # Turing.acclogp!(vi, lp)
end

@model ldamodel_vec(K, V, M, N, w, doc, beta, alpha) = begin
  theta = Matrix{Real}(undef, K, M)
  theta ~ [Dirichlet(alpha)]

  phi = Matrix{Real}(undef, V, K)
  phi ~ [Dirichlet(beta)]

  phi_dot_theta = log.(phi * theta)
  #for n = 1:N
  #  Turing.acclogp!(vi, phi_dot_theta[w[n], doc[n]])
  #end
  _lp += mapreduce(n->phi_dot_theta[w[n], doc[n]], +, 1:N)

end


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project
Pkg.add("JSON")

using JSON

string_json = """{
    "this": ["is", "json"],
    "numbers": [85, 16, 12.0],
    "and": [true, false, null]
}"""

json = JSON.parse(string_json)

json


text:
module Treesitter
using Reexport
export Parser, parse_filewith, parse_string, set_language, tree_sitter_c
export show, NodeLocation, Tree
include("LibTreesitter.jl")

@reexport using .LibTreesitter

tree_sitter_c() = ccall((:tree_sitter_c, LibTreesitter.libtreesitter_c), Ptr{TSLanguage}, ())

include("nodes.jl")
include("treecursor.jl")

include("tree.jl")
include("parser.jl")
include("debug.jl")

# Export all names starting with node, language, cursor, or tree.
foreach(names(@__MODULE__, all = true)) do s
    if startswith(string(s), r"node|language|cursor|tree"i)
        @eval export $s
    end
end

end # module Treesitter


text:
#=

This file is the entry-point for the Scattering module.

It includes this module's source files and exports the relevant keywords.

=#

module Scattering

using Parameters                # For default values in structs
using DocStringExtensions       # For simplifying docstring
using FastGaussQuadrature       # For fast Gauss-Legendre quadrature points
using JLD2                      # For saving and loading
using ProgressMeter             # For showing progress
using KernelAbstractions        # For heterogeneous (GPU+CPU) programming
using CUDA                      # For GPU programming
using Distributions             # Distributions from Julia
using ForwardDiff, DiffResults  # Automatic Differentiation tools
using LinearAlgebra             # For calculations
using StatsBase                 # Fit statistics for truncation
using StaticArrays              #

using ..Architectures: device

include("types.jl")                       # All types used in this module
include("mie_helper_functions.jl")        # Mie file-related functions
include("legendre_functions.jl")          # Recursions for associated Legendre Polynomials
include("truncate_phase.jl")              # Functions for truncation
include("compute_wigner_values.jl")       # Recursive Wigner 3j calculations
include("make_mie_model.jl")                   # Convenience functions to create Mie model
include("compute_NAI2.jl")                # Compute phase function w/ NAI2
include("compute_PCW.jl")                 # Compute phase function w/ PCW
include("phase_function_autodiff.jl")     # Auto-differentiation
include("show_utils.jl")                  # Pretty-print

# Export make functions/types
export make_mie_model, reconstruct_phase

# Export types
export NAI2, PCW, Aerosol, MieModel, Stokes_IQUV, Stokes_I, Stokes_IQU,
       δBGE, GreekCoefs, AerosolOptics, AbstractFourierDecompositionType

export compute_B, compute_ab, GreekCoefs, comp_ab, compute_mie_π_τ!,
       compute_wigner_values, save_wigner_values, load_wigner_values,
       compute_Sl, gausslegendre, compute_aerosol_optical_properties,
       compute_ref_aerosol_extinction,
       ConjugateTransposePairs, AbstractPolarizationType,
       AbstractAerosolType, AbstractAerosolType, MieModel,
       AbstractTruncationType, phase_function

end


text:
#=
    The File contains `DOK_Sparsik<T>` and related funcs implementation

    `DOK_Sparsik` implements `AbstractSparseMatrix` interface
=#
#------------------------------------------------------------------------------

"""
    DOK_Sparsik{T<:Field}

Provides basic 2-dimensional sparse vector interface, i.e Sparse Matrix.
Parametrized with the type of coefficient field `T`.

Supports these *fast* operations:

 `inner`, `scale`, `reduce`, `iterate`, matvec and matmat (`Base.prod`).

One can construct `BidimSparsik` from a dense matrix with the `from_dense` function.
"""
mutable struct DOK_Sparsik{T <: AbstractAlgebra.Field} <: AbstractSparseMatrix{T}
    # This implementation stores rows data regularly while constructing
    # columns data only when needed
    m::Int
    n::Int

    field::T

    nnz_rows::Vector{Int}
    nnz_cols::Vector{Int}

    rows::Dict{Int, Sparsik{T}}
    cols::Dict{Int, Sparsik{T}}
end

#------------------------------------------------------------------------------

Nemo.base_ring(v::DOK_Sparsik) = v.field

#------------------------------------------------------------------------------

# O(k)
function Base.deepcopy_internal(x::DOK_Sparsik, stackdict::IdDict)
    y = DOK_Sparsik(x.m, x.n, x.field,
            Base.deepcopy_internal(x.nnz_rows, stackdict),
            Base.deepcopy_internal(x.nnz_cols, stackdict),
            Base.deepcopy_internal(x.rows, stackdict),
            Base.deepcopy_internal(x.cols, stackdict))
    stackdict[x] = y
    return y
end

#------------------------------------------------------------------------------

# checks whether the `A.col` and `A.nnz_cols` fields are valid
# note that zero matrix object is always thorough
# O(1)
# Problems with data relevance when scaling/reducing inplace,
# have not broke everything, yet keep this feature in mind
function is_thorough(A::DOK_Sparsik)
    return length(A.nnz_cols) != 0 || length(A.nnz_rows) == 0
end

#------------------------------------------------------------------------------

function field(A::DOK_Sparsik)
    return A.field
end

# returns the dim of the ambient space
function dim(A::DOK_Sparsik)
    return size(A, 1) * size(A, 2)
end

# returns the order of the matrix
# `A` must be square
function order(A::DOK_Sparsik)
    return size(A, 1)
end

function issquare(A::DOK_Sparsik)
    return size(A, 1) == size(A, 2)
end

#------------------------------------------------------------------------------

# returns the tuple of `A` dimensions
# O(1)
function Base.size(A::DOK_Sparsik)
    return (A.m, A.n)
end

# returns the ith dimension of `A`
function Base.size(A::DOK_Sparsik, i::Int)
    return size(A)[i]
end

# returns number of nonzero entries of `A`
# O(k) where k is the number of nonzeroes in `A`
function Base.length(A::DOK_Sparsik)
    return sum(map(row -> length(row), values(A.rows)))
end

#------------------------------------------------------------------------------

# returns the density of `A`
# which equals the number of nonzeroes
# divied by the dim of the ambient space
# O(k) where k is the number of nonzeroes in `A`
function density(A::DOK_Sparsik)
    return length(A) / dim(A)
end

#------------------------------------------------------------------------------

# Export the matrix as a list of tripes [row_idx, col_idx, entry]
# for the nonzero entries
# In case `one_based` is set, add 1 to every index
function export_as_triples(A::DOK_Sparsik)
    result = []
    for i in A.nnz_rows
        for (j, x) in A.rows[i]
            push!(result, [i, j, x])
        end
    end
    return result
end

#------------------------------------------------------------------------------

# reconstructs fields `A.cols` and `A.nnz_cols` from rows data
# so that the resulting object is considered to be thorough
# let n = size(A, 2), k = length(A)
# O(n + k)
function reconstruct!(A::DOK_Sparsik)
    #if is_thorough(A)
    #    return A
    #end

    cols_data = [Int[] => Dict{Int, valtype(A)}()
                for _ in 1 : size(A, 2)]

    for i in A.nnz_rows
        for (j, x) in A.rows[i]
            push!(cols_data[j][1], i)
            cols_data[j][2][i] = x
        end
    end

    map(x-> unique!(sort!(x[1])), cols_data)

    for (j, col) in enumerate(cols_data)
        if !isempty(col[1])
            A.cols[j] = Sparsik(size(A, 1), A.field, col...)
            push!(A.nnz_cols, j)
        end
    end

    unique!(sort!(A.nnz_cols))

    return A
end

#------------------------------------------------------------------------------

# reconstructs fields `A.cols` and `A.nnz_cols` from rows data
# so that the resulting object is considered to be thorough
# let K = length(A)
#     k = length(A.nnz_rows)
# O(Klogk)
function reconstruct_2!(A::DOK_Sparsik)
    if is_thorough(A)
        return A
    end

    # k-way merge
    nnz_buffer = BinaryMinHeap(
                    [(first_nonzero(row), idx, 1)
                    for (idx, row) in A.rows]
                )

    # reconstructing ordered cols-nnz list
    while !isempty(nnz_buffer)
        (col_idx, row_idx, rel_idx) = pop!(nnz_buffer)
        row = A.rows[row_idx]

        if isempty(A.nnz_cols) || last(A.nnz_cols) != col_idx
            push!(A.nnz_cols, col_idx)
        end
        if rel_idx < length(row)
            rel_idx += 1
            push!(nnz_buffer, (row.nonzero[rel_idx], row_idx, rel_idx))
        end
    end

    # reconstructing cols-data dict
    for i in A.nnz_rows
        for (j, x) in A.rows[i]
            if !haskey(A.cols, j)
                A.cols[j] = zero_sparsik(size(A, 1), A.field)
            end
            push!(A.cols[j].nonzero, i)
            A.cols[j].data[i] = x
        end
    end

    return A
end

#------------------------------------------------------------------------------

# constructs `DOK_Sparsik` object from
#   `m`           - number of rows
#   `n`           - number of cols
#   `field`       - field of coefficients
#   `nnz_rows`    - ordered 1-dimensional Array of indices of nonzero rows
#   `rows`        - Dict of items (i => row) where i is
#                       the index of the nonzero row in the matrix
#                       and row - its Sparsik representation
#                   data must be converted to fmpq
#
# the resulting object is Not Thorough
# O(1)
function from_rows(m, n, field, nnz_rows, rows)
    return DOK_Sparsik(m, n, field,
                            nnz_rows, Int[],
                            rows, Dict{Int, Sparsik{typeof(field)}}())
end

#------------------------------------------------------------------------------

# constructs `DOK_Sparsik` object from
#   `m`           - number of rows
#   `n`           - number of cols
#   `nnz_coords`  - a collection-like object of items (i, j, x),
#                       where i and j mark
#                       the position of the nonzero element
#                       in the matrix and x - its value,
#                       the order is not implied
#   `field`       - field of coefficients
#
# the resulting object is Thorough
#
# O(klogk + R)
# where k is length(nnz_coords) and R is the reconstruction cost
function from_COO(m, n, nnz_coords, field)
    nnz_rows = Int[]
    rows = Dict{Int, Sparsik{typeof(field)}}()

    for (i, j, x) in nnz_coords
        if !haskey(rows, i)
            rows[i] = Sparsik(n, field, Int[], Dict{Int, elem_type(field)}())
        end
        rows[i].data[j] = field(x)
        push!(rows[i].nonzero, j)
        push!(nnz_rows, i)
    end

    unique!(sort!(nnz_rows))
    map(xs -> unique!(sort!(xs.nonzero)), values(rows))

    # and this can be speed up
    return reconstruct!(from_rows(m, n, field, nnz_rows, rows))
end

function from_COO(m, n, x::Dict, field)
    from_COO(m, n, ((k[1], k[2], x[k]) for k in keys(x)), field)
end

#------------------------------------------------------------------------------

# constructs `DOK_Sparsik` object from
#   `A`           - a collection-like object of 2 dimensions
#   `field`       - field of coefficients
#
# the resulting object is Thorough
#
# let (m, n) = size(A)
# O(mn + R)
# where R is reconstruction cost
function from_dense(A::AbstractMatrix, field)
    m, n = size(A, 1), size(A, 2)
    nnz_rows = Int[]
    rows = Dict{Int, Sparsik{typeof(field)}}()

    for i in 1 : m

        row = from_dense(A[i, :], field)

        if !iszero(row)
            push!(nnz_rows, i)
            rows[i] = row
        end
    end

    return reconstruct!(from_rows(m, n, field, nnz_rows, rows))
end

#------------------------------------------------------------------------------

# returns A + c * B
# mutates A, constructs A.rows inplace
# the resulting object is Not Thorough
#
# note that `length(A)` and `number of nonzeroes in A` are synonyms
# let k = length(A), r = length(B)
# O(k + r)
function reduce!(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}, c) where {T}
    FF = base_ring(A)
    nnz_rows = Int[]
    # equality of field types ensures that
    # valtype(A) == valtype(B)
    rows = Dict{Int, valtype(A)}()

    i, j = 1, 1
    A_nnz_rows = deepcopy(A.nnz_rows)
    B_nnz_rows = deepcopy(B.nnz_rows)
    empty!(A.nnz_rows)

    while i <= length(A_nnz_rows) || j <= length(B_nnz_rows)
        if i > length(A_nnz_rows)
            new_idx = B_nnz_rows[j]
            j += 1
        elseif j > length(B_nnz_rows)
            new_idx = A_nnz_rows[i]
            i += 1
        else
            new_idx = min(A_nnz_rows[i], B_nnz_rows[j])
            if A_nnz_rows[i] > B_nnz_rows[j]
                j += 1
            elseif A_nnz_rows[i] < B_nnz_rows[j]
                i += 1
            else
                i += 1
                j += 1
            end
        end

        if !haskey(A.rows, new_idx)
            A.rows[new_idx] = zero_sparsik(size(A, 2), A.field)
        end
        if haskey(B.rows, new_idx)
            reduce!(A.rows[new_idx], B.rows[new_idx], c)
        end

        if !iszero(A.rows[new_idx])
            push!(A.nnz_rows, new_idx)
        else
            delete!(A.rows, new_idx)
        end
    end

    return A
end

#------------------------------------------------------------------------------

# returns A + c * B
# the resulting object is Not Thorough
#
# let k = length(A), r = length(B)
# O(k + r)
function Base.reduce(A::DOK_Sparsik, B::DOK_Sparsik, c)
    return reduce!(deepcopy(A), B, c)
end

#------------------------------------------------------------------------------

# returns A × B
# the resulting object is Not Thorough
#
# let k = length(A), r = length(B)
# O(kr) if `B` Is Thorough
# O(kr + R) where R is reconstruction cost otherwise
function Base.prod(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T}
    if !is_thorough(B)
        reconstruct!(B)
    end

    FF = base_ring(A)

    nnz_rows = Int[]

    # we can create it only when needed
    # TODO: fix
    rows = Dict{Int, Sparsik{T}}()

    for i in A.nnz_rows

        # we can create it only when needed
        # TODO: fix
        row_vals = Dict{Int, valtype(A)}()
        row_indices = Int[]

        for j in B.nnz_cols
            product = inner(A.rows[i], B.cols[j])
            # TODO
            # global ZEROS
            # global TOTALS
            # TOTALS += 1
            # iszero(product) && (ZEROS += 1)

            if !iszero(product)
                push!(row_indices, j)
                row_vals[j] = product
            end
        end

        if !isempty(row_indices)
            push!(nnz_rows, i)
            rows[i] = Sparsik(size(A, 2), A.field, row_indices, row_vals)
        end
    end

    # we can return 0 instead of instantiating new object
    # as long as this method is not a part of interface
    return from_rows(
        size(A, 1), size(B, 2), FF,
        nnz_rows, rows
    )
end

function do_overlap(bv::UInt64, bu::UInt64)
    return (bv & bu) != 0
end

function prod_bloom(
    A::DOK_Sparsik{T}, B::DOK_Sparsik{T},
    bloomA::Dict{Int, UInt64}, bloomB::Dict{Int, UInt64}) where {T}

    FF = base_ring(A)
    nnz_rows = Int[]
    rows = Dict{Int, Sparsik{T}}()

    for i in A.nnz_rows
        row_vals = Dict{Int, valtype(A)}()
        row_indices = Int[]
        for j in B.nnz_cols
            if do_overlap(bloomA[i], bloomB[j])
                product = inner(A.rows[i], B.cols[j])
                if !iszero(product)
                    push!(row_indices, j)
                    row_vals[j] = product
                end
            end
        end

        if !isempty(row_indices)
            push!(nnz_rows, i)
            rows[i] = Sparsik(size(A, 2), A.field, row_indices, row_vals)
        end
    end

    return from_rows(
        size(A, 1), size(B, 2), FF,
        nnz_rows, rows
    )
end

#------------------------------------------------------------------------------

#
# returns A(v) = Av
# if k = length(A) and r = length(v)
# O(k)
function Base.prod(A::DOK_Sparsik, v::Sparsik)
    nonzero = Int[]
    field = base_ring(A)
    data = Dict{Int, elem_type(field)}()

    for idx in A.nnz_rows
        product = inner(A.rows[idx], v)
        if ! iszero(product)
            push!(nonzero, idx)
            data[idx] = product
        end
    end

    return Sparsik(size(A, 1), field, nonzero, data)
end

#------------------------------------------------------------------------------

# returns A * c modifying the `A` object itself
# the resulting object is Not Thorough
#
# let k = length(A)
# O(k)
function scale!(A::DOK_Sparsik, c)
    if iszero(c)
        return empty!(A)
    end

    for i in A.nnz_rows
        scale!(A.rows[i], c)
    end

    return A
end

#------------------------------------------------------------------------------

# returns A * c, does not modify `A`
# the resulting object is Not Thorough
#
# let k = length(A)
# O(k)
function scale(A::DOK_Sparsik, c)
    return scale!(deepcopy(A), c)
end

#------------------------------------------------------------------------------

# returns A[i, j] or zero if not present
# O(1)
function Base.get(A::DOK_Sparsik, i::Int, j::Int)
    if !haskey(A.rows, i)
        return zero(A.field)
    end
    return A.rows[i][j]
end

function Base.haskey(A::DOK_Sparsik, i::Int, j::Int)
    return haskey(A.rows, i) && haskey(A.rows[i].data, j)
end

function Base.haskey(A::DOK_Sparsik, i::Int)
    return haskey(A, to_cartesian(A, i)...)
end

#------------------------------------------------------------------------------

# erases all entries of `A`
# the resulting object is Thorough
# (since empty object is)
#
# let k = length(A)
# O(k)
function Base.empty!(A::DOK_Sparsik)
    empty!(A.nnz_rows)
    empty!(A.rows)
    empty!(A.nnz_cols)
    empty!(A.cols)
    return A
end

#------------------------------------------------------------------------------

# returns zero matrix of the dimensions (m, n)
function zero_sparsik(m, n, field)
    from_rows(m, n, field, Int[], Dict{Int, Sparsik{typeof(field)}}())
    #zero(field)
end

zero_sparsik(sz::Tuple{Int, Int}, field) = zero_sparsik(sz..., field)

# returns zero matrix of the dimensions of `A`
Base.zero(A::DOK_Sparsik) = zero_sparsik(size(A)..., A.field)

Base.iszero(A::DOK_Sparsik{T}) where {T} = length(A.nnz_rows) == 0

#------------------------------------------------------------------------------

==(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T} = (size(A) == size(B) &&
            A.field == B.field;
            A.rows == B.rows;
            A.nnz_rows == B.nnz_rows)
!=(A::DOK_Sparsik{T}, u::DOK_Sparsik{T}) where {T} = !(A == u)

+(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T} = reduce(A, B, one(A.field))
-(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T} = reduce(A, B, - one(A.field))
-(A::DOK_Sparsik) = scale(A, -one(A.field))

*(A::DOK_Sparsik, c::Number) = scale(A, c)
*(c::Number, A::DOK_Sparsik) = scale(A, c)
*(A::DOK_Sparsik, c::AbstractAlgebra.FieldElem)= scale(A, c)
*(c::AbstractAlgebra.FieldElem, A::DOK_Sparsik) = scale(A, c)
*(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T} = prod(A, B)

*(A::DOK_Sparsik{T}, v::Sparsik{T}) where {T} = prod(A, v)

#------------------------------------------------------------------------------

Base.repr(::MIME"text/plain", A::DOK_Sparsik) = "$(join(map(i -> "|$(join(map(j -> get(A, i, j), 1 : A.n), "\t"))|", 1 : A.m), "\n"))"
Base.show(io::IO, A::DOK_Sparsik) = print(io, repr(MIME("text/plain"), A))

#------------------------------------------------------------------------------

function to_cartesian(A::DOK_Sparsik, idx::Int)
    return to_cartesian(size(A), idx)
 end

 function to_plain(A, i::Int, j::Int)
     return to_plain(size(A), i, j)
 end

function to_cartesian(sz::Tuple{Int, Int}, idx::Int)
    cols = sz[2]
    if mod(idx, cols) == 0
        i = div(idx, cols)
    else
        i = div(idx, cols) + 1
    end
        j = idx - cols * (i - 1)
    return (i, j)
 end

function to_plain(sz::Tuple{Int, Int}, i::Int, j::Int)
    return (i - 1) * sz[1] + j
end

#------------------------------------------------------------------------------

# returns the idx-th element of m if m is treated as
# a vector stretched across the rows
# O(1)
Base.getindex(m::DOK_Sparsik, idx::Int) = m[to_cartesian(m, idx)...]

# returns the (i, j) element of m as a matrix
# O(1)
Base.getindex(m::DOK_Sparsik, i::Int, j::Int) = get(m, i, j)

function first_nonzero(m::DOK_Sparsik)
    if iszero(m)
        return -1
    end
    return (m.nnz_rows[1] - 1) * size(m, 2) + first_nonzero(m.rows[m.nnz_rows[1]])
end

#------------------------------------------------------------------------------

Base.valtype(A::DOK_Sparsik) = elem_type(A.field)
Base.eltype(A::DOK_Sparsik) = (Int, valtype(A))

#------------------------------------------------------------------------------

# iterates over A nonzeroes
# each returned item is of type eltype(A)
# linear in the nimber of nonzeroes in A
function Base.iterate(A::DOK_Sparsik, state=(1, 1))
    (i, j) = state
    if i > length(A.nnz_rows) || j > length(A.rows[A.nnz_rows[i]])
        (i, j) = (i + 1, 1)
    end
    if i > length(A.nnz_rows)
        return nothing
    else
        r_idx = A.nnz_rows[i]
        c_idx = A.rows[r_idx].nonzero[j]
        state = (i, j + 1)
        return ((to_plain(A, r_idx, c_idx), A[r_idx, c_idx]), state)
    end
end

#-----------------------------------------------------------------------------

# returns a new AbstractSparsik object consisting of elements
# from the given `A` each converted to the `field`
#
# O(k) if k is the number of nonzero in `A`
function modular_reduction(A::DOK_Sparsik, field)
    new_nnz_rows = Int[]
    new_rows = Dict{Int, Sparsik{typeof(field)}}()

    for i in A.nnz_rows
        y = modular_reduction(A.rows[i], field)
        if !iszero(y)
            push!(new_nnz_rows, i)
            new_rows[i] = y
        end
    end

    return from_rows(size(A)..., field, new_nnz_rows, new_rows)
end

#-----------------------------------------------------------------------------

# returns a new DOK_Sparsik object consisting of elements
# from the `A` each reconstructed from A.field to QQ
#
# roughly
# O(klog²(ch)) if k is the number of nonzero in `v`
# and ch is the field characteristic
function rational_reconstruction(A::DOK_Sparsik)
    new_nnz_rows = Int[]
    new_rows = Dict{Int, Sparsik{typeof(Nemo.QQ)}}()

    for i in A.nnz_rows
        y = rational_reconstruction(A.rows[i])
        if !iszero(y)
            push!(new_nnz_rows, i)
            new_rows[i] = y
        end
    end

    return from_rows(size(A)..., Nemo.QQ, new_nnz_rows, new_rows)
end

#-----------------------------------------------------------------------------

# returns standard inner product of A and B
# if the vectors A and B have k and r nonzeroes respectively
# O(min(k, r))
function inner(A::DOK_Sparsik{T}, B::DOK_Sparsik{T}) where {T}
    ans = zero(A.field)

    if length(A) > length(B)
        inner(B, A)
    end

    for (i, row) in A.rows
        if haskey(B.rows, i)
            ans += inner(row, B.rows[i])
        end
    end

    return ans
end

# ad-hoc inner
function inner(A::DOK_Sparsik{T}, c::AbstractAlgebra.FieldElem) where {T}
    if iszero(c)
        return zero(c)
    end

    error("not implemented!")
end

inner(c::AbstractAlgebra.FieldElem, A::DOK_Sparsik{T}) where {T} = inner(A, c)

#-----------------------------------------------------------------------------

# returns the Id matrix
# this one can be pretty slow
function Base.one(A::DOK_Sparsik{T}) where {T}
    field = base_ring(A)
    reconstruct!(from_rows(
        size(A)...,
        field,
        Array(1:order(A)),
        Dict{Int, Sparsik{T}}(
            i => unit_sparsik(order(A), i, field)
            for i in 1:order(A)
        )
    ))
end

#-----------------------------------------------------------------------------

function random_sparsik(sz::Tuple{Int, Int}, field; density=0.1)
    @assert 0 <= density <= 1

    # outer product
    ×(xs, ys) = (i, j) -> (xs[i], ys[j])

    underlying = (1:sz[1]) × (1:sz[2])
    n = prod(sz)
    λ = density * n

    nnz = Dict{Tuple{Int, Int}, elem_type(field)}()
    for _ in 1:λ
        while (idx = to_cartesian(sz, rand(1:n))) != nothing
            !haskey(nnz, idx) && break
        end
        nnz[idx...] = rand(field)
    end

    return from_COO(sz..., nnz, field)
end

#-----------------------------------------------------------------------------

# vectorizes the matrix
# O(k) if k is nnz count
function Base.vec(A::DOK_Sparsik)
    m, n = size(A)
    field = base_ring(A)
    nonzero = Int[]
    data = Dict{Int, elem_type(field)}()
    for idx in A.nnz_rows
        for (i, x) in A.rows[idx]
            data[i + (idx - 1) * n] = x
            push!(nonzero, i + (idx - 1) * n)
        end
    end

    return Sparsik(m*n, field, nonzero, data)
end

#-----------------------------------------------------------------------------

# restricts the given DOK_Sparsik object to the given `coords`
# if the sparsik is treated as a vector
# `coords` - a collection-like object of Ints
#
# O(length(coords))
function restrict(A::DOK_Sparsik, coords)
    nonzero = Int[]
    field = base_ring(A)
    data = Dict{Int, elem_type(field)}()
    for (i, idx) in enumerate(coords)
        if !iszero(A[idx])
            push!(nonzero, i)
            data[i] = A[idx]
        end
    end
    return Sparsik(length(coords), field, nonzero, data)
end

# the same as the above
# BUT !
# O(length(coords))
function restrict(A::DOK_Sparsik, coords, first_value)
    nonzero = Int[]
    field = base_ring(A)
    data = Dict{Int, elem_type(field)}()

    if ! iszero(first_value)
        push!(nonzero, 1)
        data[1] = first_value
    end

    for (i, idx) in enumerate(coords)
        if !iszero(A[idx])
            push!(nonzero, i + 1)
            data[i + 1] = A[idx]
        end
    end
    return Sparsik(length(coords) + 1, field, nonzero, data)
end

# -----------------------------------------------------------------------------

# well, it transposes the matrix
# O(R)
# where R is the reconstruction cost
function transpose!(A::DOK_Sparsik)
    reconstruct!(A)
    A.m, A.n = A.n, A.m
    A.nnz_rows, A.nnz_cols = A.nnz_cols, A.nnz_rows
    A.rows, A.cols = A.cols, A.rows
    return A
end

# -----------------------------------------------------------------------------

function to_dense(A::DOK_Sparsik)
    ans = fill(zero(base_ring(A)), size(A))
    for (i, x) in A
        ans[to_cartesian(A, i)...] = x
    end
    return ans
end

# -----------------------------------------------------------------------------

function unit_sparsik(sz::Tuple{Int, Int}, idx::Int, field)
    i, j = to_cartesian(sz, idx)
    return DOK_Sparsik(
        sz...,
        field,
        [ i ],
        Int[],
        Dict(i => unit_sparsik(sz[2], j, field)),
        Dict{Int, Sparsik{typeof(field)}}()
    )
end

# -----------------------------------------------------------------------------

function tr(A::DOK_Sparsik)
    if iszero(A)
        return zero(base_ring(A))
    end
    sum(A[i, i] for i in A.nnz_rows)
end


text:
# Julia wrapper for header: /usr/include/scip/scip_solvingstats.h
# Automatically generated using Clang.jl wrap_c

function SCIPgetNRuns(scip)
    ccall((:SCIPgetNRuns, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNReoptRuns(scip)
    ccall((:SCIPgetNReoptRuns, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPaddNNodes(scip, nnodes)
    ccall((:SCIPaddNNodes, libscip), Cvoid, (Ptr{SCIP_}, Clonglong), scip, nnodes)
end

function SCIPgetNNodes(scip)
    ccall((:SCIPgetNNodes, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNTotalNodes(scip)
    ccall((:SCIPgetNTotalNodes, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNFeasibleLeaves(scip)
    ccall((:SCIPgetNFeasibleLeaves, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNInfeasibleLeaves(scip)
    ccall((:SCIPgetNInfeasibleLeaves, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNObjlimLeaves(scip)
    ccall((:SCIPgetNObjlimLeaves, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDelayedCutoffs(scip)
    ccall((:SCIPgetNDelayedCutoffs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNLPs(scip)
    ccall((:SCIPgetNLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNLPIterations(scip)
    ccall((:SCIPgetNLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNNZs(scip)
    ccall((:SCIPgetNNZs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNRootLPIterations(scip)
    ccall((:SCIPgetNRootLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNRootFirstLPIterations(scip)
    ccall((:SCIPgetNRootFirstLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrimalLPs(scip)
    ccall((:SCIPgetNPrimalLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrimalLPIterations(scip)
    ccall((:SCIPgetNPrimalLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDualLPs(scip)
    ccall((:SCIPgetNDualLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDualLPIterations(scip)
    ccall((:SCIPgetNDualLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNBarrierLPs(scip)
    ccall((:SCIPgetNBarrierLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNBarrierLPIterations(scip)
    ccall((:SCIPgetNBarrierLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNResolveLPs(scip)
    ccall((:SCIPgetNResolveLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNResolveLPIterations(scip)
    ccall((:SCIPgetNResolveLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrimalResolveLPs(scip)
    ccall((:SCIPgetNPrimalResolveLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrimalResolveLPIterations(scip)
    ccall((:SCIPgetNPrimalResolveLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDualResolveLPs(scip)
    ccall((:SCIPgetNDualResolveLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDualResolveLPIterations(scip)
    ccall((:SCIPgetNDualResolveLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNNodeLPs(scip)
    ccall((:SCIPgetNNodeLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNNodeLPIterations(scip)
    ccall((:SCIPgetNNodeLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNNodeInitLPs(scip)
    ccall((:SCIPgetNNodeInitLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNNodeInitLPIterations(scip)
    ccall((:SCIPgetNNodeInitLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDivingLPs(scip)
    ccall((:SCIPgetNDivingLPs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNDivingLPIterations(scip)
    ccall((:SCIPgetNDivingLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNStrongbranchs(scip)
    ccall((:SCIPgetNStrongbranchs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNStrongbranchLPIterations(scip)
    ccall((:SCIPgetNStrongbranchLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNRootStrongbranchs(scip)
    ccall((:SCIPgetNRootStrongbranchs, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNRootStrongbranchLPIterations(scip)
    ccall((:SCIPgetNRootStrongbranchLPIterations, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNPriceRounds(scip)
    ccall((:SCIPgetNPriceRounds, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPricevars(scip)
    ccall((:SCIPgetNPricevars, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPricevarsFound(scip)
    ccall((:SCIPgetNPricevarsFound, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPricevarsApplied(scip)
    ccall((:SCIPgetNPricevarsApplied, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNSepaRounds(scip)
    ccall((:SCIPgetNSepaRounds, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNCutsFound(scip)
    ccall((:SCIPgetNCutsFound, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNCutsFoundRound(scip)
    ccall((:SCIPgetNCutsFoundRound, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNCutsApplied(scip)
    ccall((:SCIPgetNCutsApplied, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNConflictConssFound(scip)
    ccall((:SCIPgetNConflictConssFound, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNConflictConssFoundNode(scip)
    ccall((:SCIPgetNConflictConssFoundNode, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNConflictConssApplied(scip)
    ccall((:SCIPgetNConflictConssApplied, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetMaxDepth(scip)
    ccall((:SCIPgetMaxDepth, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetMaxTotalDepth(scip)
    ccall((:SCIPgetMaxTotalDepth, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNBacktracks(scip)
    ccall((:SCIPgetNBacktracks, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNActiveConss(scip)
    ccall((:SCIPgetNActiveConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNEnabledConss(scip)
    ccall((:SCIPgetNEnabledConss, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgDualbound(scip)
    ccall((:SCIPgetAvgDualbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgLowerbound(scip)
    ccall((:SCIPgetAvgLowerbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetDualbound(scip)
    ccall((:SCIPgetDualbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLowerbound(scip)
    ccall((:SCIPgetLowerbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetDualboundRoot(scip)
    ccall((:SCIPgetDualboundRoot, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetLowerboundRoot(scip)
    ccall((:SCIPgetLowerboundRoot, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetFirstLPDualboundRoot(scip)
    ccall((:SCIPgetFirstLPDualboundRoot, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetFirstLPLowerboundRoot(scip)
    ccall((:SCIPgetFirstLPLowerboundRoot, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetFirstPrimalBound(scip)
    ccall((:SCIPgetFirstPrimalBound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetPrimalbound(scip)
    ccall((:SCIPgetPrimalbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetUpperbound(scip)
    ccall((:SCIPgetUpperbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetCutoffbound(scip)
    ccall((:SCIPgetCutoffbound, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPupdateCutoffbound(scip, cutoffbound)
    ccall((:SCIPupdateCutoffbound, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cdouble), scip, cutoffbound)
end

function SCIPisPrimalboundSol(scip)
    ccall((:SCIPisPrimalboundSol, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetGap(scip)
    ccall((:SCIPgetGap, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetTransGap(scip)
    ccall((:SCIPgetTransGap, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetNSolsFound(scip)
    ccall((:SCIPgetNSolsFound, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNLimSolsFound(scip)
    ccall((:SCIPgetNLimSolsFound, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetNBestSolsFound(scip)
    ccall((:SCIPgetNBestSolsFound, libscip), Clonglong, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgPseudocost(scip, solvaldelta)
    ccall((:SCIPgetAvgPseudocost, libscip), Cdouble, (Ptr{SCIP_}, Cdouble), scip, solvaldelta)
end

function SCIPgetAvgPseudocostCurrentRun(scip, solvaldelta)
    ccall((:SCIPgetAvgPseudocostCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, Cdouble), scip, solvaldelta)
end

function SCIPgetAvgPseudocostCount(scip, dir)
    ccall((:SCIPgetAvgPseudocostCount, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetAvgPseudocostCountCurrentRun(scip, dir)
    ccall((:SCIPgetAvgPseudocostCountCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetPseudocostCount(scip, dir, onlycurrentrun)
    ccall((:SCIPgetPseudocostCount, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR, UInt32), scip, dir, onlycurrentrun)
end

function SCIPgetAvgPseudocostScore(scip)
    ccall((:SCIPgetAvgPseudocostScore, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetPseudocostVariance(scip, branchdir, onlycurrentrun)
    ccall((:SCIPgetPseudocostVariance, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR, UInt32), scip, branchdir, onlycurrentrun)
end

function SCIPgetAvgPseudocostScoreCurrentRun(scip)
    ccall((:SCIPgetAvgPseudocostScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgConflictScore(scip)
    ccall((:SCIPgetAvgConflictScore, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgConflictScoreCurrentRun(scip)
    ccall((:SCIPgetAvgConflictScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgConflictlengthScore(scip)
    ccall((:SCIPgetAvgConflictlengthScore, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgConflictlengthScoreCurrentRun(scip)
    ccall((:SCIPgetAvgConflictlengthScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgInferences(scip, dir)
    ccall((:SCIPgetAvgInferences, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetAvgInferencesCurrentRun(scip, dir)
    ccall((:SCIPgetAvgInferencesCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetAvgInferenceScore(scip)
    ccall((:SCIPgetAvgInferenceScore, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgInferenceScoreCurrentRun(scip)
    ccall((:SCIPgetAvgInferenceScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgCutoffs(scip, dir)
    ccall((:SCIPgetAvgCutoffs, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetAvgCutoffsCurrentRun(scip, dir)
    ccall((:SCIPgetAvgCutoffsCurrentRun, libscip), Cdouble, (Ptr{SCIP_}, SCIP_BRANCHDIR), scip, dir)
end

function SCIPgetAvgCutoffScore(scip)
    ccall((:SCIPgetAvgCutoffScore, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetAvgCutoffScoreCurrentRun(scip)
    ccall((:SCIPgetAvgCutoffScoreCurrentRun, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPgetDeterministicTime(scip)
    ccall((:SCIPgetDeterministicTime, libscip), Cdouble, (Ptr{SCIP_},), scip)
end

function SCIPprintOrigProblem(scip, file, extension, genericnames)
    ccall((:SCIPprintOrigProblem, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Cstring, UInt32), scip, file, extension, genericnames)
end

function SCIPprintTransProblem(scip, file, extension, genericnames)
    ccall((:SCIPprintTransProblem, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, Cstring, UInt32), scip, file, extension, genericnames)
end

function SCIPprintStatusStatistics(scip, file)
    ccall((:SCIPprintStatusStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintTimingStatistics(scip, file)
    ccall((:SCIPprintTimingStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintOrigProblemStatistics(scip, file)
    ccall((:SCIPprintOrigProblemStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintTransProblemStatistics(scip, file)
    ccall((:SCIPprintTransProblemStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintPresolverStatistics(scip, file)
    ccall((:SCIPprintPresolverStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintConstraintStatistics(scip, file)
    ccall((:SCIPprintConstraintStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintConstraintTimingStatistics(scip, file)
    ccall((:SCIPprintConstraintTimingStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintPropagatorStatistics(scip, file)
    ccall((:SCIPprintPropagatorStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintConflictStatistics(scip, file)
    ccall((:SCIPprintConflictStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintSeparatorStatistics(scip, file)
    ccall((:SCIPprintSeparatorStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintPricerStatistics(scip, file)
    ccall((:SCIPprintPricerStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintBranchruleStatistics(scip, file)
    ccall((:SCIPprintBranchruleStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintHeuristicStatistics(scip, file)
    ccall((:SCIPprintHeuristicStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintCompressionStatistics(scip, file)
    ccall((:SCIPprintCompressionStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintLPStatistics(scip, file)
    ccall((:SCIPprintLPStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintNLPStatistics(scip, file)
    ccall((:SCIPprintNLPStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintRelaxatorStatistics(scip, file)
    ccall((:SCIPprintRelaxatorStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintTreeStatistics(scip, file)
    ccall((:SCIPprintTreeStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintRootStatistics(scip, file)
    ccall((:SCIPprintRootStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintSolutionStatistics(scip, file)
    ccall((:SCIPprintSolutionStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintConcsolverStatistics(scip, file)
    ccall((:SCIPprintConcsolverStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintBendersStatistics(scip, file)
    ccall((:SCIPprintBendersStatistics, libscip), Cvoid, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintStatistics(scip, file)
    ccall((:SCIPprintStatistics, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintReoptStatistics(scip, file)
    ccall((:SCIPprintReoptStatistics, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintBranchingStatistics(scip, file)
    ccall((:SCIPprintBranchingStatistics, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}), scip, file)
end

function SCIPprintDisplayLine(scip, file, verblevel, endline)
    ccall((:SCIPprintDisplayLine, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{FILE}, SCIP_VERBLEVEL, UInt32), scip, file, verblevel, endline)
end

function SCIPgetNImplications(scip)
    ccall((:SCIPgetNImplications, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPwriteImplicationConflictGraph(scip, filename)
    ccall((:SCIPwriteImplicationConflictGraph, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring), scip, filename)
end

function SCIPstoreSolutionGap(scip)
    ccall((:SCIPstoreSolutionGap, libscip), Cvoid, (Ptr{SCIP_},), scip)
end


text:
# module TestIPython

if lowercase(get(ENV, "CI", "false")) == "true"
    let
        if VERSION < v"0.7.0-"
            setup_code = ""
        else
            setup_code = Base.load_path_setup_code()
        end
        path = joinpath(@__DIR__, "install_dependencies.jl")
        code = """
        $setup_code
        include("$(escape_string(path))")
        """
        run(`$(Base.julia_cmd()) -e $code`)
        # Run install_dependencies.jl in a separate process since it
        # may re-build PyCall.  In that case, we need to load
        # re-precompiled PyCall and IPython.
    end
end

include("preamble.jl")

IPython.envinfo()

ipy_opts = @time IPython._start_ipython(:ipython_options)
ipy_main = ipy_opts["user_ns"]["Main"]

@testset "Main" begin
    _setproperty!(ipy_main, :x, 17061)
    @test x == 17061
end

include("test_julia_repl.jl")
include("test_convenience.jl")

IPython.test_ipython_jl(inprocess=true)

# end  # module


text:
module TSne

# Numpy Math.sum => axis = 0 => sum down the columns. axis = 1 => sum along the rows
# Julia Base.sum => axis = 1 => sum down the columns. axis = 2 => sum along the rows

#
#  tsne.jl
#
#  This is a straight off Julia port of Laurens van der Maatens python implementation of tsne

export tsne, pca

using Distances

function Hbeta(D, beta = 1.0)
	#Compute the perplexity and the P-row for a specific value of the precision of a Gaussian distribution.
	P = exp(-copy(D) * beta)
	sumP = sum(P)
	H = log(sumP) + beta * sum(D .* P) / sumP
	P = P / sumP
	return (H, P)
end

function x2p(X, tol = 1e-5, perplexity = 30.0)
	#Performs a binary search to get P-values in such a way that each conditional Gaussian has the same perplexity.
	println("Computing pairwise distances...")
	(n, d) = size(X)
	sum_X = sum((X.^2),2)
	D = (-2 * (X * X') .+ sum_X)' .+ sum_X
	P = zeros(n, n)
	beta = ones(n, 1)
	logU = log(perplexity)

	# Loop over all datapoints
	range = [1:n]
	for i in 1:n

		# Print progress
		if mod(i, 500) == 0
			println("Computing P-values for point " * string(i) *  " of " * string(n) * "...")
        	end

		# Compute the Gaussian kernel and entropy for the current precision
		betamin = -Inf
		betamax =  Inf

		inds = range[range .!=i]
		Di = D[i, inds]
		(H, thisP) = Hbeta(Di, beta[i])

		# Evaluate whether the perplexity is within tolerance
		Hdiff = H - logU
		tries = 0
		while abs(Hdiff) > tol && tries < 50

			# If not, increase or decrease precision
			if Hdiff > 0
				betamin = beta[i]
				if betamax == Inf || betamax == -Inf
					beta[i] = beta[i] * 2
				else
					beta[i] = (beta[i] + betamax) / 2
				end
			else
				betamax = beta[i]
				if betamin == Inf || betamin == -Inf
					beta[i] = beta[i] / 2
				else
					beta[i] = (beta[i] + betamin) / 2
				end
			end

			# Recompute the values
			(H, thisP) = Hbeta(Di, beta[i])
			Hdiff = H - logU
			tries = tries + 1
		end
		# Set the final row of P
  		P[i, inds] = thisP

	end
	# Return final P-matrix
	println("Mean value of sigma: " * string(mean(sqrt(1 ./ beta))))
	return P
end

function pca(X, no_dims = 50)
	#Runs PCA on the NxD array X in order to reduce its dimensionality to no_dims dimensions.

	println("Preprocessing the data using PCA...")
	(n, d) = size(X)
	X = X - repmat(mean(X, 1), n, 1)
	C = (X' * X) ./ (size(X,1)-1)
	(l, M) = eig(C)
	sorder = sortperm(l,rev=true)
	M = M[:,sorder]
	ret_dims = no_dims > d ? d : no_dims
	Y = X * M[:,1:ret_dims]
	return Y
end

function grad!(dY, Y, P, n, no_dims, Q)
  fill!(dY, 0.0)
  pairwise!(Q, SqEuclidean(), Y')
  sum_Q = 0.0
  for col = 1:n, row = 1:n
    if col != row
      Q[row, col] = 1 / (1 + Q[row, col])
      sum_Q += Q[row, col]
    end
  end
  for col = 1:n, row = 1:n
    mult = (P[row, col] - (Q[row, col] / sum_Q)) * Q[row, col]
    for d = 1:no_dims
      dY[col, d] += (Y[col, d] - Y[row, d]) * mult
    end
  end
end

function cost!(P, Q, n)
  # destroys Q
  sum_Q = sum(Q)
  removeNA = (x) -> isnan(x) ? 0.0 : x
  for col = 1:n, row = 1:n
    v = P[row, col] * log(P[row, col] / max(Q[row, col] /  sum_Q, 1e-12))
    Q[row, col] = removeNA(v)
  end
  sum(Q)
end

function update!(Y, dY, iY, gains, momentum, min_gain, eta, n, no_dims)
  mean_Y = zeros(no_dims)
  for d = 1:no_dims, row = 1:n
    gains[row, d] = (gains[row, d] + 0.2) * ((dY[row, d] .> 0) != (iY[row, d] > 0)) +
                    (gains[row, d] * 0.8) * ((dY[row, d] .> 0) == (iY[row, d] > 0))
    gains[row, d] = gains[row, d] < min_gain ? min_gain : gains[row, d]
    iY[row, d] = momentum * iY[row, d] - eta * (gains[row, d] * dY[row, d])
    Y[row, d] = Y[row, d] + iY[row, d]
    mean_Y[d] += Y[row, d]
  end
  for d = 1:no_dims
    mean_Y[d] /= n
  end
  for d= 1:no_dims, row = 1:n
    Y[row, d] = Y[row, d] - mean_Y[d]
  end
end

function tsne(X, no_dims = 2, initial_dims = -1, max_iter = 1000, perplexity = 30.0)
	#Runs t-SNE on the dataset in the NxD array X to reduce its dimensionality to no_dims dimensions.
	# Diffrent from orginal, default is to not use PCA
	println("Initial X Shape is : " * string(size(X)))

	# Initialize variables
	if(initial_dims>0)
		X = pca(X, initial_dims)
	end
        # zero the mean
        X = X .- mean(X)
        # scale
        X = X / std(X)
	(n, d) = size(X)
	initial_momentum = 0.5
	final_momentum = 0.8
	eta = 500
	min_gain = 0.01
	Y = randn(n, no_dims) .* 0.0001
	dY = zeros(n, no_dims)
	iY = zeros(n, no_dims)
	gains = ones(n, no_dims)
        Q = zeros(n, n)

	# Compute P-values
	P = x2p(X, 1e-5, perplexity)
	P = P + P'
	P = P / sum(P)
	P = P * 12						# early exaggeration
	P = max(P, 1e-12)

	# Run iterations
	for iter in 1:max_iter
                grad!(dY, Y, P, n, no_dims, Q)

		# Perform the update
		if iter <= 250
			momentum = initial_momentum
		else
			momentum = final_momentum
		end
                update!(Y, dY, iY, gains, momentum, min_gain, eta, n, no_dims)

		# Compute current value of cost function
		if mod((iter + 1), 10) == 0
			println("Iteration ", (iter + 1), ": error is ", cost!(P, Q, n))
		end
		# Stop lying about P-values
		if iter == 250
			P = P / 12
		end
    	end
	# Return solution
	return Y
	end

end


text:
if !isdefined(@__MODULE__, :IoK8sApimachineryPkgApisMetaV1Time)
    const IoK8sApimachineryPkgApisMetaV1Time = String
else
    @warn("Skipping redefinition of IoK8sApimachineryPkgApisMetaV1Time to DateTime")
end


text:
using GR

using DelimitedFiles

function main()
    z = Z = readdlm("sans.dat")
    G = [ exp(-x^2 -y^2) for x in LinRange(-1.5, 1.5, 128), y in LinRange(-1.5, 1.5, 128) ]

    for t = 0:500
        surface(z, title="Time: $t s")
        z += 0.05 * Z .* G .* rand(Float64, (128, 128))
    end
end

main()


text:
# src/portfolios.jl
# Abstract Portfolios

export contains, quantity_contained, balance

"""```
contains(portfolio::AbstractPortfolio{K, N}, asset::K, quantity::N)
```
Check if the `portfolio` contains `quantity` amount of the given `asset`.
"""
function contains(portfolio::AbstractPortfolio{K,N}, asset::K, quantity::N)::Bool where {K,N}
    missing_api("contains", portfolio, asset, quantity)
end

"""```
contains(wallet::AbstractWallet{K, N}, asset::K, quantity::N)
```
Check if the `wallet` contains `quantity` amount of the given `asset`.
"""
function contains(wallet::AbstractWallet{K,N}, asset::K, quantity::N)::Bool where {K,N}
    return contains(portfolio(wallet), asset, quantity)
end

"""```
quantity_contained(portfolio::AbstractPortfolio{K}, asset::K)
```
Calculate the amount of `asset` contained in the `portfolio`. If the asset is not
in the wallet, a suitable null value should be returned.
"""
function quantity_contained(portfolio::AbstractPortfolio{K,N}, asset::K)::N where {K,N}
    missing_api("quantity_contained", portfolio, asset)
end

"""```
quantity_contained(wallet::AbstractWallet{K}, asset::K)
```
Calculate the amount of `asset` contained in the `wallet`. If the asset is not
in the wallet, a suitable null value should be returned.
"""
function quantity_contained(wallet::AbstractWallet{K,N}, asset::K)::N where {K,N}
    return quantity_contained(portfolio(wallet), asset)
end

"""```
balance(portfolio::AbstractPortfolio{K}, asset::K, base)
```
Calculate the total value of a given `asset` contained in the `portfolio` in the
given `base` asset (not necessarily contained in the `portfolio`).
"""
function balance(portfolio::AbstractPortfolio{K}, asset::K, base) where {K}
    missing_api("balance", portfolio, asset, base)
end

"""```
balance(wallet::AbstractWallet{K}, asset::K, base)
```
Calculate the total value of a given `asset` contained in the `wallet` in the
given `base` asset (not necessarily contained in the `wallet`).
"""
function balance(wallet::AbstractWallet{K}, asset::K, base) where {K}
    return balance(portfolio(wallet), asset, base)
end

"""```
balance(portfolio::AbstractPortfolio, base)
```
Calculate the total value contained in the `portfolio` in the given `base` asset
(not necessarily contained in the `portfolio`).
"""
function balance(portfolio::AbstractPortfolio, base)
    return sum((asset, _) -> balance(portfolio, asset, base), portfolio)
end

"""```
balance(wallet::AbstractWallet, base)
```
Calculate the total value contained in the `wallet` in the given `base` asset
(not necessarily contained in the `wallet`).
"""
function balance(wallet::AbstractWallet, base)
    return balance(portfolio(wallet), base)
end


text:
# Model
include_model("hopper3D")
model = free_time_model(model)

# Horizon
T = 51

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_uu[end] = 2.0 * h
_ul[end] = 0.5 * h
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp_init = MRP(UnitQuaternion(RotZ(0.0) * RotY(0.0) * RotX(0.0)))

z_h = 0.25
q1 = [0.0, 0.0, 0.5 + z_h, mrp_init.x, mrp_init.y, mrp_init.z, 0.25]
x1 = [q1; q1]
# qT = [0.0, 0.0, 0.5 + z_h, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
# xT = [qT; qT]

xl, xu = state_bounds(model, T,
    [model.qL; model.qL], [model.qU; model.qU],
    x1 = [q1; Inf * ones(nq)])

# Objective
Qq = Diagonal([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 1.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0, 1.0, 1.0, zeros(model.m - model.nu - 1)..., 0.0])

obj_tracking = quadratic_time_tracking_objective(
    [t < T ? 0.0 * Q : 0.0 * QT for t = 1:T],
    [R for t = 1:T-1],
    [[q1; q1] for t = 1:T],
    [zeros(model.m) for t = 1:T],
    1.0)

obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect(4:6))

obj_penalty = PenaltyObjective(1.0e5, model.m-1)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity])

# Constraints
include_constraints(["contact", "free_time", "loop"])
con_contact = contact_constraints(model_ft, T)
con_free_time = free_time_constraints(T)
con_loop = loop_constraints(model_ft, (1:model_ft.n), 1, T)
con = multiple_constraints([con_contact, con_free_time, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = [[q1; q1] for t = 1:T] #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [1.0e-3 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time z̄, info = solve(prob, copy(z0), tol = 1.0e-6, c_tol = 1.0e-6, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
tf, t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = mean(h̄)
@save joinpath(pwd(), "examples/trajectories/hopper3D_vertical_gait.jld2") z̄ x̄ ū h̄ q u γ b

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = ū[1][end])


text:
import StaticArrays: SVector, MVector
import DelayEmbeddings: Dataset

"""
    triangulate(pts)

Triangulate a set of points and, depending on the type of the input,
return an appropriate instance of a subtype of `AbstracTriangulationPartition`.

## Arguments
- **`pts`**. The points. May be an instance of `AbstractArray{T, 2`} (the largest dimension will be taken as the index), `Vector{Vector}`, `Vector{SVector}`, `Vector{MVector}` (each vector is a point) or a `Dataset`.
"""
function triangulate end

"""
    triangulate_m(pts)

Triangulate a set of points and, depending on the type of the input,
return an appropriate mutable instance of a subtype of `AbstractMutableTriangulationPartition`.

## Arguments
- **`pts`**. The points. May be an instance of `AbstractArray{T, 2`} (the largest dimension will be taken as the index), `Vector{Vector}`, `Vector{SVector}`, `Vector{MVector}` (each vector is a point) or a `Dataset`.
"""
function triangulate_m end

"""
    triangulate_f(pts)

Triangulate a set of points and, depending on the type of the input,
return an appropriate instance of a subtype of `AbstractTriangulationPartitionFull`,
which not only contains the points and the Delaunay triangulation, but also
the simplices and their radii, centroids, orientations and volumes.

## Arguments
- **`pts`**. The points. May be an instance of `AbstractArray{T, 2`} (the largest dimension will be taken as the index), `Vector{Vector}`, `Vector{SVector}`, `Vector{MVector}` (each vector is a point) or a `Dataset`.
"""
function triangulate_f end

"""
    triangulate_mf(pts)

Triangulate a set of points and, depending on the type of the input,
return an appropriate mutable instance of a subtype of `AbstractMutableTriangulationPartitionFull`,
which not only contains the points and the Delaunay triangulation, but also
the simplices and their radii, centroids, orientations and volumes.

## Arguments
- **`pts`**. The points. May be an instance of `AbstractArray{T, 2`} (the largest dimension will be taken as the index), `Vector{Vector}`, `Vector{SVector}`, `Vector{MVector}` (each vector is a point) or a `Dataset`.
"""
function triangulate_mf end

function triangulate(pts::Vector{Vector{T}}) where T
    DT = DelaunayTriangulation(pts)
    dim = length(pts[1])
    PointTriangulation{dim, T}(pts, DT)
end

function triangulate_m(pts::Vector{Vector{T}}) where T
    DT = DelaunayTriangulation(pts)
    dim = length(pts[1])
    MutablePointTriangulation{dim, T}(pts, DT)
end

function triangulate_f(pts::Vector{Vector{T}}) where T
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    PointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations, volumes, centroids)
end

function triangulate_mf(pts::Vector{Vector{T}}) where T
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MutablePointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate(pts::AbstractArray{T, 2}) where T
    DT = DelaunayTriangulation(pts)
    dim = minimum(size(pts))

    if size(pts, 1) > size(pts, 2)
        pts = transpose(pts)
    end
    MatrixTriangulation{dim, T}(pts, DT)
end

function triangulate_m(pts::AbstractArray{T, 2}) where T
    DT = DelaunayTriangulation(pts)
    dim = minimum(size(pts))

    if size(pts, 1) > size(pts, 2)
        pts = transpose(pts)
    end

    MutableMatrixTriangulation{dim, T}(pts, DT)
end

function triangulate_f(pts::AbstractArray{T, 2}) where T
    DT = DelaunayTriangulation(pts)

    if size(pts, 1) > size(pts, 2)
        pts = transpose(pts)
    end

    n_simplices = length(DT)
    dim = minimum(size(pts))
    simplices = [Simplex(pts[:, DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MatrixTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate_mf(pts::AbstractArray{T, 2}) where T
    DT = DelaunayTriangulation(pts)

    if size(pts, 1) > size(pts, 2)
        pts = transpose(pts)
    end

    n_simplices = length(DT)
    dim = minimum(size(pts))
    simplices = [Simplex(pts[:, DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MutableMatrixTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate(pts::Dataset{D, T}) where {D, T}
    DT = DelaunayTriangulation(pts)
    DatasetTriangulation{D, T}(pts, DT)
end

function triangulate_m(pts::Dataset{D, T}) where {D, T}
    DT = DelaunayTriangulation(pts)
    MutableDatasetTriangulation{D, T}(pts, DT)
end

function triangulate_f(pts::Dataset{D, T}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = minimum(size(pts))
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MatrixTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate_mf(pts::Dataset{D, T}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = minimum(size(pts))
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MutableMatrixTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate(pts::Vector{SVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    PointTriangulation{D, T}(pts, DT)
end

function triangulate_m(pts::Vector{SVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    MutablePointTriangulation{D, T}(pts, DT)
end

function triangulate_f(pts::Vector{SVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    PointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate_mf(pts::Vector{SVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MutablePointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate(pts::Vector{MVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    PointTriangulation{D, T}(pts, DT)
end

function triangulate_m(pts::Vector{MVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    MutablePointTriangulation{D, T}(pts, DT)
end

function triangulate_f(pts::Vector{MVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    PointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

function triangulate_mf(pts::Vector{MVector{D, T}}) where {D, T}
    DT = DelaunayTriangulation(pts)
    n_simplices = length(DT)
    dim = length(pts[1])
    simplices = [Simplex(pts[DT[i]]) for i = 1:n_simplices]
    radii = radius.(simplices)
    centroids = centroid.(simplices)
    orientations = orientation.(simplices)
    volumes = volume.(simplices)

    MutablePointTriangulationFull{dim, T}(pts, DT, simplices, radii, orientations,
        volumes, centroids)
end

export
triangulate,
triangulate_m,
triangulate_f,
triangulate_mf