text:
using StructArrays
using StructArrays: staticschema, iscompatible, _promote_typejoin, append!!
using OffsetArrays: OffsetArray
using StaticArrays
import Tables, PooledArrays, WeakRefStrings
using TypedTables: Table
using DataAPI: refarray, refvalue
using Adapt: adapt, Adapt
using Test

using Documenter: doctest
doctest(StructArrays)

@testset "index" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray((a = a, b = b))
    @test (@inferred t[2,2]) == (a = 4, b = 7)
    @test (@inferred t[2,1:2]) == StructArray((a = [3, 4], b = [6, 7]))
    @test_throws BoundsError t[3,3]
    @test (@inferred view(t, 2, 1:2)) == StructArray((a = view(a, 2, 1:2), b = view(b, 2, 1:2)))
end

@testset "components" begin
    t = StructArray(a = 1:10, b = rand(Bool, 10))
    @test StructArrays.propertynames(t) == (:a, :b)
    @test StructArrays.propertynames(StructArrays.components(t)) == (:a, :b)
end

@testset "utils" begin
    t = StructArray(rand(ComplexF64, 2, 2))
    T = staticschema(eltype(t))
    @test StructArrays.eltypes(T) == NamedTuple{(:re, :im), Tuple{Float64, Float64}}
    @test StructArrays.map_params(eltype, T) == NamedTuple{(:re, :im), Tuple{Float64, Float64}}
    @test StructArrays.map_params(eltype, StructArrays.astuple(T)) == Tuple{Float64, Float64}
    @test !iscompatible(typeof((1, 2)), typeof(([1],)))
    @test iscompatible(typeof((1, 2)), typeof(([1], [2])))
    @test !iscompatible(typeof((1, 2)), typeof(([1.1], [2])))
    @test iscompatible(typeof(()), typeof(()))
    @test _promote_typejoin(Tuple{Int, Missing}, Tuple{Int, Int}) == Tuple{Int, Union{Int, Missing}}
    @test _promote_typejoin(Pair{Int, Missing}, Pair{Int, Int}) == Pair{Int, Union{Int, Missing}}
    @test _promote_typejoin(NamedTuple{(:a, :b), Tuple{Int, Missing}}, NamedTuple{(:a, :b), Tuple{Int, Int}}) == NamedTuple{(:a, :b), Tuple{Int, Union{Int, Missing}}}
    @test _promote_typejoin(Tuple{}, Tuple{}) == Tuple{}
    @test _promote_typejoin(Tuple{Int}, Tuple{Int, Int}) == Tuple{Int, Vararg{Int, N} where N}

    @test StructArrays.astuple(Tuple{Int}) == Tuple{Int}
    @test StructArrays.strip_params(Tuple{Int}) == Tuple
    @test StructArrays.astuple(NamedTuple{(:a,), Tuple{Float64}}) == Tuple{Float64}
    @test StructArrays.strip_params(NamedTuple{(:a,), Tuple{Float64}}) == NamedTuple{(:a,)}
end

@testset "indexstyle" begin
    s = StructArray(a=rand(10,10), b=view(rand(100,100), 1:10, 1:10))
    T = typeof(s)
    @test IndexStyle(T) === IndexCartesian()
    @test StructArrays.index_type(s) == CartesianIndex{2}
    @test s[100] == s[10, 10] == (a=s.a[10,10], b=s.b[10,10])
    s[100] = (a=1, b=1)
    @test s[100] == s[10, 10] == (a=1, b=1)
    s[10, 10] = (a=0, b=0)
    @test s[100] == s[10, 10] == (a=0, b=0)
    @inferred IndexStyle(StructArray(a=rand(10,10), b=rand(10,10)))
    s = StructArray(a=rand(10,10), b=rand(10,10))
    T = typeof(s)
    @test StructArrays.index_type(s) == Int
    @inferred IndexStyle(s)
    @test s[100] == s[10, 10] == (a=s.a[10,10], b=s.b[10,10])
    s[100] = (a=1, b=1)
    @test s[100] == s[10, 10] == (a=1, b=1)
    s[10, 10] = (a=0, b=0)
    @test s[100] == s[10, 10] == (a=0, b=0)

    # inference for "many" types, both for linear ad Cartesian indexing
    @inferred StructArrays.index_type(ntuple(_ -> rand(5), 2))
    @inferred StructArrays.index_type(ntuple(_ -> rand(5, 5), 3))
    @inferred StructArrays.index_type(ntuple(_ -> rand(5, 5, 5), 4))

    @inferred StructArrays.index_type(ntuple(_ -> view(rand(5), 1:3), 2))
    @inferred StructArrays.index_type(ntuple(_ -> view(rand(5, 5), 1:3, 1:2), 3))
    @inferred StructArrays.index_type(ntuple(_ -> view(rand(5, 5, 5), 1:3, 1:2, 1:4), 4))

    @inferred StructArrays.index_type(ntuple(n -> n == 1 ? rand(5, 5) : view(rand(5, 5), 1:2, 1:3), 5))
    @inferred IndexStyle(StructArray(a=rand(10,10), b=view(rand(100,100), 1:10, 1:10)))
end

@testset "replace_storage" begin
    v = StructArray(a=rand(10), b = fill("string", 10))
    v_pooled = StructArrays.replace_storage(v) do c
        isbitstype(eltype(c)) ? c : convert(PooledArrays.PooledArray, c)
    end
    @test eltype(v) == eltype(v_pooled)
    @test all(v.a .== v_pooled.a)
    @test all(v.b .== v_pooled.b)
    @test !isa(v_pooled.a, PooledArrays.PooledArray)
    @test isa(v_pooled.b, PooledArrays.PooledArray)
end

@testset "roweq" begin
    a = ["a", "b", "a", "a"]
    b = PooledArrays.PooledArray(["x", "y", "z", "x"])
    s = StructArray((a, b))
    @test StructArrays.roweq(s, 1, 1)
    @test !StructArrays.roweq(s, 1, 2)
    @test !StructArrays.roweq(s, 1, 3)
    @test StructArrays.roweq(s, 1, 4)
    strs = WeakRefStrings.StringArray(["a", "a", "b"])
    @test StructArrays.roweq(strs, 1, 1)
    @test StructArrays.roweq(strs, 1, 2)
    @test !StructArrays.roweq(strs, 1, 3)
    @test !StructArrays.roweq(strs, 2, 3)
end

@testset "permute" begin
    a = WeakRefStrings.StringVector(["a", "b", "c"])
    b = PooledArrays.PooledArray([1, 2, 3])
    c = [:a, :b, :c]
    s = StructArray(a=a, b=b, c=c)
    permute!(s, [2, 3, 1])
    @test s.a == ["b", "c", "a"]
    @test s.b == [2, 3, 1]
    @test s.c == [:b, :c, :a]
    s = StructArray(a=[1, 2], b=["a", "b"])
    t = StructArray(a=[3, 4], b=["c", "d"])
    copyto!(s, t)
    @test s == t

    s = StructArray(a=[1, 2], b=["a", "b"])
    t = StructArray(a=[3, 4], b=["c", "d"])
    copyto!(s, 1, t, 1, 2)
    @test s == t

    a = WeakRefStrings.StringVector(["a", "b", "c"])
    b = PooledArrays.PooledArray(["1", "2", "3"])
    c = [:a, :b, :c]
    s = StructArray(a=a, b=b, c=c)
    ref = refarray(s)
    Base.permute!!(ref, sortperm(s))
    @test issorted(s)
end

@testset "sortperm" begin
    c = StructArray(a=[1,1,2,2], b=[1,2,3,3], c=["a","b","c","d"])
    d = StructArray(a=[1,1,2,2], b=[1,2,3,3], c=["a","b","c","d"])
    @test issorted(c)
    @test sortperm(c) == [1,2,3,4]
    permute!(c, [2,3,4,1])
    @test c == StructArray(a=[1,2,2,1], b=[2,3,3,1], c=["b","c","d","a"])
    @test sortperm(c) == [4,1,2,3]
    @test !issorted(c)
    @test sort(c) == d
    sort!(c)
    @test c == d

    c = StructArray(a=[1,1,2,2], b=[1,2,3,3], c=PooledArrays.PooledArray(["a","b","c","d"]))
    d = StructArray(a=[1,1,2,2], b=[1,2,3,3], c=PooledArrays.PooledArray(["a","b","c","d"]))
    @test issorted(c)
    @test sortperm(c) == [1,2,3,4]
    permute!(c, [2,3,4,1])
    @test c == StructArray(a=[1,2,2,1], b=[2,3,3,1], c=PooledArrays.PooledArray(["b","c","d","a"]))
    @test sortperm(c) == [4,1,2,3]
    @test !issorted(c)
    @test sort(c) == d
    sort!(c)
    @test c == d
end

@testset "iterators" begin
    c = [1, 2, 3, 1, 1]
    d = StructArrays.GroupPerm(c)
    @test parent(d) == c
    @test eltype(d) == UnitRange{Int}
    @test Base.IteratorEltype(d) == Base.HasEltype()
    @test sortperm(d) == sortperm(c)
    s = collect(d)
    @test s == [1:3, 4:4, 5:5]
    t = collect(StructArrays.finduniquesorted(c))
    @test first.(t) == [1, 2, 3]
    @test last.(t) == [[1, 4, 5], [2], [3]]
    u = collect(StructArrays.uniquesorted(c))
    @test u == [1, 2, 3]
end

@testset "similar" begin
    t = StructArray(a = rand(10), b = rand(Bool, 10))
    s = similar(t)
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (10,)
    t = StructArray(a = rand(10, 2), b = rand(Bool, 10, 2))
    s = similar(t, 3, 5)
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (3, 5)
    s = similar(t, (3, 5))
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (3, 5)
end

@testset "empty" begin
    s = StructVector(a = [1, 2, 3], b = ["a", "b", "c"])
    empty!(s)
    @test isempty(s.a)
    @test isempty(s.b)
end

@testset "constructor from existing array" begin
    v = rand(ComplexF64, 5, 3)
    t = @inferred StructArray(v)
    @test size(t) == (5, 3)
    @test t[2,2] == v[2,2]
    t2 = convert(StructArray, v)::StructArray
    @test t2 == t
    t3 = StructArray(t)::StructArray
    @test t3 == t
    @test convert(StructArray, t) == t

    v = rand(ComplexF64, 5)
    t = @inferred StructVector(v)
    @test t[2] == v[2]
    @test size(t) == (5,)
    @test t == convert(StructVector, v)
    @test t == convert(StructVector, t)
end

@testset "tuple case" begin
    s = StructArray(([1], ["test"],))
    @test s[1] == (1, "test")
    @test Base.getproperty(s, 1) == [1]
    @test Base.getproperty(s, 2) == ["test"]
    t = StructArray{Tuple{Int, Float64}}(([1], [1.2]))
    @test t[1] == (1, 1.2)

    t[1] = (2, 3)
    @test t[1] == (2, 3.0)
    push!(t, (1, 2))
    @test getproperty(t, 1) == [2, 1]
    @test getproperty(t, 2) == [3.0, 2.0]
    @test pop!(t) == (1, 2.0)
    @test getproperty(t, 1) == [2]
    @test getproperty(t, 2) == [3.0]
end

@testset "constructor from slices" begin
    if VERSION >= v"1.1"
        X = [1.0 2.0; 3.0 4.0]
        @test StructArray{Complex{Float64}}(X; dims=1) == [Complex(1.0,3.0), Complex(2.0,4.0)]
        @test StructArray{Complex{Float64}}(X; dims=2) == [Complex(1.0,2.0), Complex(3.0,4.0)]

        X = [1.0 2.0; 3.0 4.0; 5.0 6.0]
        @test StructArray{Tuple{Float64,Complex{Float64}}}(X; dims=1) == [(1.0,Complex(3.0,5.0)), (2.0, Complex(4.0,6.0))]
    end
end

struct A
    x::Int
    y::Int
    A(x) = new(x, x)
end

@testset "internal constructor" begin
    v = A.([1, 2, 3])
    s = StructArray(v)
    @test s[1] == A(1)
    @test s[2] == A(2)
    @test s[3] == A(3)
end

@testset "kwargs constructor" begin
    a = [1.2]
    b = [2.3]
    @test StructArray(a=a, b=b) == StructArray((a=a, b=b))
    @test StructArray{ComplexF64}(re=a, im=b) == StructArray{ComplexF64}((a, b))
    f1() = StructArray(a=[1.2], b=["test"])
    f2() = StructArray{Pair}(first=[1.2], second=["test"])
    t1 = @inferred f1()
    t2 = @inferred f2()
    @test t1 == StructArray((a=[1.2], b=["test"]))
    @test t2 == StructArray{Pair}(([1.2], ["test"]))
end

@testset "complex" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray{ComplexF64}((a, b))
    @test t[2,2] == ComplexF64(4, 7)
    @test t[2,1:2] == StructArray{ComplexF64}(([3, 4], [6, 7]))
    @test view(t, 2, 1:2) == StructArray{ComplexF64}((view(a, 2, 1:2), view(b, 2, 1:2)))
end

@testset "copy" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray{ComplexF64}((a, b))
    t2 = @inferred copy(t)
    @test t2[1,1] == 1.0 + im*4.0
    t2[1,1] = 2.0 + im*4.0
    # Test we actually did a copy
    @test t[1,1] == 1.0 + im*4.0
    # Test that `copy` works, even when the array type changes (e.g. views)
    s = rand(10, 2)
    v = StructArray{ComplexF64}((view(s, :, 1), view(s, :, 2)))
    v2 = copy(v)
    @test v2.re isa Vector
    @test v2.im isa Vector
    @test v2.re == copy(v.re)
    @test v2.im == copy(v.im)
end

@testset "undef initializer" begin
    t = @inferred StructArray{ComplexF64}(undef, 5, 5)
    @test eltype(t) == ComplexF64
    @test size(t) == (5,5)
    c = 2 + im
    t[1,1] = c
    @test t[1,1] == c
end

@testset "resize!" begin
    t = StructArray{Pair}(([3, 5], ["a", "b"]))
    resize!(t, 5)
    @test length(t) == 5
    p = 1 => "c"
    t[5] = p
    @test t[5] == p
end

@testset "concat" begin
    t = StructArray{Pair}(([3, 5], ["a", "b"]))
    push!(t, (2 => "d"))
    @test t == StructArray{Pair}(([3, 5, 2], ["a", "b", "d"]))
    @test pop!(t) == (2 => "d")
    push!(t, (2 => "c"))
    @test t == StructArray{Pair}(([3, 5, 2], ["a", "b", "c"]))
    append!(t, t)
    @test t == StructArray{Pair}(([3, 5, 2, 3, 5, 2], ["a", "b", "c", "a", "b", "c"]))
    t = StructArray{Pair}(([3, 5], ["a", "b"]))
    t2 = StructArray{Pair}(([1, 6], ["a", "b"]))
    @test cat(t, t2; dims=1)::StructArray == StructArray{Pair}(([3, 5, 1, 6], ["a", "b", "a", "b"])) == vcat(t, t2)
    @test vcat(t, t2) isa StructArray
    @test cat(t, t2; dims=2)::StructArray == StructArray{Pair}(([3 1; 5 6], ["a" "a"; "b" "b"])) == hcat(t, t2)
    @test hcat(t, t2) isa StructArray
end

f_infer() = StructArray{ComplexF64}((rand(2,2), rand(2,2)))
g_infer() = StructArray([(a=(b="1",), c=2)], unwrap = t -> t <: NamedTuple)
tup_infer() = StructArray([(1, 2), (3, 4)])
cols_infer() = StructArray(([1, 2], [1.2, 2.3]))

@testset "inferrability" begin
    @inferred f_infer()
    @inferred g_infer()
    @test g_infer().a.b == ["1"]
    s = @inferred tup_infer()
    @test StructArrays.components(s) == ([1, 3], [2, 4])
    @test s[1] == (1, 2)
    @test s[2] == (3, 4)
    @inferred cols_infer()
end

@testset "propertynames" begin
    a = StructArray{ComplexF64}((Float64[], Float64[]))
    @test sort(collect(propertynames(a))) == [:im, :re]
end

@testset "tables" begin
    s = StructArray([(a=1, b="test")])
    @test Tables.schema(s) == Tables.Schema((:a, :b), (Int, String))
    @test Tables.rows(s)[1].a == 1
    @test Tables.rows(s)[1].b == "test"
    @test Tables.columns(s).a == [1]
    @test Tables.columns(s).b == ["test"]
    @test Tables.istable(s)
    @test Tables.istable(typeof(s))
    @test Tables.rowaccess(s)
    @test Tables.rowaccess(typeof(s))
    @test Tables.columnaccess(s)
    @test Tables.columnaccess(typeof(s))
    @test Tables.getcolumn(Tables.columns(s), 1) == [1]
    @test Tables.getcolumn(Tables.columns(s), :a) == [1]
    @test Tables.getcolumn(Tables.columns(s), 2) == ["test"]
    @test Tables.getcolumn(Tables.columns(s), :b) == ["test"]
    @test append!(StructArray([1im]), [(re = 111, im = 222)]) ==
        StructArray([1im, 111 + 222im])
    @test append!(StructArray([1im]), (x for x in [(re = 111, im = 222)])) ==
        StructArray([1im, 111 + 222im])
    @test append!(StructArray([1im]), Table(re = [111], im = [222])) ==
        StructArray([1im, 111 + 222im])
    # Testing integer column "names":
    @test invoke(append!, Tuple{StructVector,Any}, StructArray(([0],)), StructArray(([1],))) ==
        StructArray(([0, 1],))
end

struct S
    x::Int
    y::Float64
    S(x) = new(x, x)
end

StructArrays.createinstance(::Type{<:S}, x, y) = S(x)

@testset "inner" begin
    v = StructArray{S}(([1], [1]))
    @test v[1] == S(1)
    @test v[1].y isa Float64
end

@testset "arrayof" begin
    v = StructArrays.arrayof(Missing, (2,))
    @test v isa Array{Missing, 1}
    v = StructArrays.arrayof(Int, (2,))
    @test v isa Array{Int, 1}
end

const initializer = StructArrays.ArrayInitializer(t -> t <: Union{Tuple, NamedTuple, Pair})
collect_structarray_rec(t) = collect_structarray(t, initializer = initializer)

@testset "collectnamedtuples" begin
    v = [(a = 1, b = 2), (a = 1, b = 3)]
    @test collect_structarray_rec(v) == StructArray((a = Int[1, 1], b = Int[2, 3]))

    # test inferrability with constant eltype
    itr = [(a = 1, b = 2), (a = 1, b = 2), (a = 1, b = 12)]
    el, st = iterate(itr)
    dest = initializer(typeof(el), (3,))
    dest[1] = el
    @inferred StructArrays.collect_to_structarray!(dest, itr, 2, st)

    v = [(a = 1, b = 2), (a = 1.2, b = 3)]
    @test collect_structarray_rec(v) == StructArray((a = [1, 1.2], b = Int[2, 3]))
    @test typeof(collect_structarray_rec(v)) == typeof(StructArray((a = Real[1, 1.2], b = Int[2, 3])))
    @test StructArray(v[i] for i in eachindex(v)) == StructArray((a = [1, 1.2], b = Int[2, 3]))

    s = StructArray(a = [1, 2], b  = [3, 4])
    @test collect_structarray(LazyRow(s, i) for i in eachindex(s)) == s
    @test collect_structarray_rec(LazyRow(s, i) for i in eachindex(s)) == s

    v = [(a = 1, b = 2), (a = 1.2, b = "3")]
    @test collect_structarray_rec(v) == StructArray((a = [1, 1.2], b = Any[2, "3"]))
    @test typeof(collect_structarray_rec(v)) == typeof(StructArray((a = Real[1, 1.2], b = Any[2, "3"])))

    v = [(a = 1, b = 2), (a = 1.2, b = 2), (a = 1, b = "3")]
    @test collect_structarray_rec(v) == StructArray((a = Real[1, 1.2, 1], b = Any[2, 2, "3"]))
    @test typeof(collect_structarray_rec(v)) == typeof(StructArray((a = Real[1, 1.2, 1], b = Any[2, 2, "3"])))

    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray((a = [2, 4, 6, 8], b = [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? (a = 1.2, b =i-1) : (a = i+1, b = i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr_real) == StructArray((a = Real[1.2, 4, 6, 8], b = [0, 2, 4, 6]))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray((a = Int[], b = Int[]))

    itr = (i for i in 0:-1)
    tuple_itr = ((a = i+1, b = i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray((a = Int[], b = Int[]))
end

@testset "collecttuples" begin
    v = [(1, 2), (1, 3)]
    @test collect_structarray_rec(v) == StructArray((Int[1, 1], Int[2, 3]))
    @inferred collect_structarray_rec(v)

    @test collect_structarray(v) == StructArray((Int[1, 1], Int[2, 3]))
    @inferred collect_structarray(v)

    v = [(1, 2), (1.2, 3)]
    @test collect_structarray_rec(v) == StructArray((Real[1, 1.2], Int[2, 3]))

    v = [(1, 2), (1.2, "3")]
    @test collect_structarray_rec(v) == StructArray((Real[1, 1.2], Any[2, "3"]))
    @test typeof(collect_structarray_rec(v)) == typeof(StructArray((Real[1, 1.2], Any[2, "3"])))

    v = [(1, 2), (1.2, 2), (1, "3")]
    @test collect_structarray_rec(v) == StructArray((Real[1, 1.2, 1], Any[2, 2, "3"]))
    # length unknown
    itr = Iterators.filter(isodd, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray(([2, 4, 6, 8], [0, 2, 4, 6]))
    tuple_itr_real = (i == 1 ? (1.2, i-1) : (i+1, i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr_real) == StructArray((Real[1.2, 4, 6, 8], [0, 2, 4, 6]))
    @test typeof(collect_structarray_rec(tuple_itr_real)) == typeof(StructArray((Real[1.2, 4, 6, 8], [0, 2, 4, 6])))

    # empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray((Int[], Int[]))

    itr = (i for i in 0:-1)
    tuple_itr = ((i+1, i-1) for i in itr)
    @test collect_structarray_rec(tuple_itr) == StructArray((Int[], Int[]))
end

@testset "collectscalars" begin
    v = (i for i in 1:3)
    @test collect_structarray_rec(v) == [1,2,3]
    @inferred collect_structarray_rec(v)

    v = (i == 1 ? 1.2 : i for i in 1:3)
    @test collect_structarray_rec(v) == collect(v)

    itr = Iterators.filter(isodd, 1:100)
    @test collect_structarray_rec(itr) == collect(itr)
    real_itr = (i == 1 ? 1.5 : i for i in itr)
    @test collect_structarray_rec(real_itr) == collect(real_itr)
    @test eltype(collect_structarray_rec(real_itr)) == Real

    #empty
    itr = Iterators.filter(t -> t > 10, 1:8)
    tuple_itr = (exp(i) for i in itr)
    @test collect_structarray_rec(tuple_itr) == Float64[]

    itr = (i for i in 0:-1)
    tuple_itr = (exp(i) for i in itr)
    @test collect_structarray_rec(tuple_itr) == Float64[]

    t = collect_structarray_rec((a = i,) for i in (1, missing, 3))
    @test StructArrays.components(t)[1] isa Array{Union{Int, Missing}}
    @test isequal(StructArrays.components(t)[1], [1, missing, 3])
end

pair_structarray((first, last)) = StructArray{Pair{eltype(first), eltype(last)}}((first, last))

@testset "collectpairs" begin
    v = (i=>i+1 for i in 1:3)
    @test collect_structarray_rec(v) == pair_structarray([1,2,3] => [2,3,4])
    @test eltype(collect_structarray_rec(v)) == Pair{Int, Int}

    v = (i == 1 ? (1.2 => i+1) : (i => i+1) for i in 1:3)
    @test collect_structarray_rec(v) == pair_structarray([1.2,2,3] => [2,3,4])
    @test eltype(collect_structarray_rec(v)) == Pair{Real, Int}

    v = ((a=i,) => (b="a$i",) for i in 1:3)
    @test collect_structarray_rec(v) == pair_structarray(StructArray((a = [1,2,3],)) => StructArray((b = ["a1","a2","a3"],)))
    @test eltype(collect_structarray_rec(v)) == Pair{NamedTuple{(:a,), Tuple{Int64}}, NamedTuple{(:b,), Tuple{String}}}

    v = (i == 1 ? (a="1",) => (b="a$i",) : (a=i,) => (b="a$i",) for i in 1:3)
    @test collect_structarray_rec(v) == pair_structarray(StructArray((a = ["1",2,3],)) => StructArray((b = ["a1","a2","a3"],)))
    @test eltype(collect_structarray_rec(v)) == Pair{NamedTuple{(:a,), Tuple{Any}}, NamedTuple{(:b,), Tuple{String}}}

    # empty
    v = ((a=i,) => (b="a$i",) for i in 0:-1)
    @test collect_structarray_rec(v) == pair_structarray(StructArray((a = Int[],)) => StructArray((b = String[],)))
    @test eltype(collect_structarray_rec(v)) == Pair{NamedTuple{(:a,), Tuple{Int}}, NamedTuple{(:b,), Tuple{String}}}

    v = Iterators.filter(t -> t.first.a == 4, ((a=i,) => (b="a$i",) for i in 1:3))
    @test collect_structarray_rec(v) == pair_structarray(StructArray((a = Int[],)) => StructArray((b = String[],)))
    @test eltype(collect_structarray_rec(v)) == Pair{NamedTuple{(:a,), Tuple{Int}}, NamedTuple{(:b,), Tuple{String}}}

    t = collect_structarray_rec((b = 1,) => (a = i,) for i in (2, missing, 3))
    s = pair_structarray(StructArray(b = [1,1,1]) => StructArray(a = [2, missing, 3]))
    @test s[1] == t[1]
    @test ismissing(t[2].second.a)
    @test s[3] == t[3]
end

@testset "collect2D" begin
    s = (l for l in [(a=i, b=j) for i in 1:3, j in 1:4])
    v = collect_structarray(s)
    @test size(v) == (3, 4)
    @test v.a == [i for i in 1:3, j in 1:4]
    @test v.b == [j for i in 1:3, j in 1:4]
end

@testset "collectoffset" begin
    s = OffsetArray([(a=1,) for i in 1:10], -3)
    sa = StructArray(s)
    @test sa isa StructArray
    @test axes(sa) == (-2:7,)
    @test sa.a == fill(1, -2:7)

    sa = StructArray(i for i in s)
    @test sa isa StructArray
    @test axes(sa) == (-2:7,)
    @test sa.a == fill(1, -2:7)

    zero_origin(T, d) = OffsetArray{T}(undef, map(r -> r .- 1, d))
    sa = collect_structarray(
        [(a = 1,), (a = 2,), (a = 3,)],
        initializer = StructArrays.StructArrayInitializer(t -> false, zero_origin),
    )
    @test sa isa StructArray
    @test collect(sa.a) == 1:3
    @test sa.a isa OffsetArray

    sa = collect_structarray(
        (x for x in [(a = 1,), (a = 2,), (a = 3,)] if true),
        initializer = StructArrays.StructArrayInitializer(t -> false, zero_origin),
    )
    @test sa isa StructArray
    @test collect(sa.a) == 1:3
    @test sa.a isa OffsetArray
end

@testset "collectstructarrays" begin
    s = StructArray(a = rand(10), b = rand(10))
    t = StructArray(a = rand(10), b = rand(10))
    v = StructArray([s, t])
    @test v.a[1] == s.a
    @test v.a[2] == t.a
    @test v.b[1] == s.b
    @test v.b[2] == t.b
    @test v[1] == s
    @test v[2] == t

    s = LazyRows(StructArray(a = rand(10), b = rand(10)))
    t = LazyRows(StructArray(a = rand(10), b = rand(10)))
    v = StructArray([s, t])
    @test v.a[1] == s.a
    @test v.a[2] == t.a
    @test v.b[1] == s.b
    @test v.b[2] == t.b
    @test v[1] == s
    @test v[2] == t
end

@testset "hasfields" begin
    @test StructArrays.hasfields(ComplexF64)
    @test !StructArrays.hasfields(Any)
    @test StructArrays.hasfields(Tuple{Union{Int, Missing}})
    @test StructArrays.hasfields(typeof((a=1,)))
    @test !StructArrays.hasfields(NamedTuple)
    @test !StructArrays.hasfields(Tuple{Int, Vararg{Int, N}} where {N})
    @test StructArrays.hasfields(Missing)
    @test !StructArrays.hasfields(Union{Tuple{Int}, Missing})
    @test StructArrays.hasfields(Nothing)
    @test !StructArrays.hasfields(Union{Tuple{Int}, Nothing})
end

@testset "reshape" begin
    s = StructArray(a=[1,2,3,4], b=["a","b","c","d"])
    rs = reshape(s, (2, 2))
    @test rs.a == [1 3; 2 4]
    @test rs.b == ["a" "c"; "b" "d"]
end

@testset "lazy" begin
    s = StructArray{ComplexF64}((rand(10, 10), view(rand(100, 100), 1:10, 1:10)))
    rows = LazyRows(s)
    @test propertynames(rows) == (:re, :im)
    @test propertynames(rows[1]) == (:re, :im)
    @test staticschema(eltype(rows)) == staticschema(ComplexF64)
    @test getproperty(rows, 1) isa Matrix{Float64}
    @test getproperty(rows, :re) isa Matrix{Float64}
    @test IndexStyle(rows) isa IndexCartesian
    @test IndexStyle(typeof(rows)) isa IndexCartesian
    @test all(t -> t.re >= 0, s)
    @test all(t -> t.re >= 0, rows)
    rows[13].re = -12
    rows[13].im = 0

    s = StructArray(rand(ComplexF64, 10, 10))
    rows = LazyRows(s)
    @test propertynames(rows) == (:re, :im)
    @test propertynames(rows[1]) == (:re, :im)
    @test staticschema(eltype(rows)) == staticschema(ComplexF64)
    @test getproperty(rows, 1) isa Matrix{Float64}
    @test getproperty(rows, :re) isa Matrix{Float64}
    @test IndexStyle(rows) isa IndexLinear
    @test IndexStyle(typeof(rows)) isa IndexLinear
    @test all(t -> t.re >= 0, s)
    @test all(t -> t.re >= 0, rows)
    rows[13].re = -12
    rows[13].im = 0
    @test !all(t -> t.re >= 0, s)
    @test !all(t -> t.re >= 0, rows)

    @test !all(t -> t.re >= 0, s)
    @test !all(t -> t.re >= 0, rows)
    io = IOBuffer()
    show(io, rows[13])
    str = String(take!(io))
    @test str == "LazyRow(re = -12.0, im = 0.0)"

    io = IOBuffer()
    Base.showarg(io, rows, true)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "LazyRows(::Array{Float64,2}, ::Array{Float64,2}) with eltype LazyRow{Complex{Float64}}"
    else
        @test str == "LazyRows(::Matrix{Float64}, ::Matrix{Float64}) with eltype LazyRow{ComplexF64}"
    end
    io = IOBuffer()
    Base.showarg(io, rows, false)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "LazyRows(::Array{Float64,2}, ::Array{Float64,2})"
    else
        @test str == "LazyRows(::Matrix{Float64}, ::Matrix{Float64})"
    end
    s = StructArray((rand(10, 10), rand(10, 10)))
    rows = LazyRows(s)
    @test IndexStyle(rows) isa IndexLinear
    @test IndexStyle(typeof(rows)) isa IndexLinear
    @test all(t -> Tables.getcolumn(t, 1) >= 0, s)
    @test all(t -> getproperty(t, 1) >= 0, rows)
    setproperty!(rows[13], 1, -12)
    setproperty!(rows[13], 2, 0)
    @test !all(t -> Tables.getcolumn(t, 1) >= 0, s)
    @test !all(t -> getproperty(t, 1) >= 0, rows)

    io = IOBuffer()
    show(io, rows[13])
    str = String(take!(io))
    @test str == "LazyRow(-12.0, 0.0)"

    io = IOBuffer()
    Base.showarg(io, rows, true)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "LazyRows(::Array{Float64,2}, ::Array{Float64,2}) with eltype LazyRow{Tuple{Float64,Float64}}"
    else
        @test str == "LazyRows(::Matrix{Float64}, ::Matrix{Float64}) with eltype LazyRow{Tuple{Float64, Float64}}"
    end
    io = IOBuffer()
    Base.showarg(io, rows, false)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "LazyRows(::Array{Float64,2}, ::Array{Float64,2})"
    else
        @test str == "LazyRows(::Matrix{Float64}, ::Matrix{Float64})"
    end
end

@testset "refarray" begin
    s = PooledArrays.PooledArray(["a", "b", "c", "c"])
    @test refarray(s) == UInt8.([1, 2, 3, 3])

    s = WeakRefStrings.StringArray(["a", "b"])
    @test refarray(s) isa WeakRefStrings.StringArray{WeakRefStrings.WeakRefString{UInt8}}
    @test all(isequal.(s, refarray(s)))
    s = WeakRefStrings.StringArray(["a", missing])
    @test refarray(s) isa WeakRefStrings.StringArray{Union{WeakRefStrings.WeakRefString{UInt8}, Missing}}
    @test all(isequal.(s, refarray(s)))
    a = WeakRefStrings.StringVector(["a", "b", "c"])
    b = PooledArrays.PooledArray(["1", "2", "3"])
    c = [:a, :b, :c]
    s = StructArray(a=a, b=b, c=c)
    ref = refarray(s)
    @test ref[1].a isa WeakRefStrings.WeakRefString{UInt8}
    @test ref[1].b isa Integer
    for i in 1:3
        @test refvalue(s, ref[i]) == s[i]
    end
end

@testset "show" begin
    s = StructArray(Complex{Int64}[1+im, 2-im])
    io = IOBuffer()
    Base.showarg(io, s, true)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "StructArray(::Array{Int64,1}, ::Array{Int64,1}) with eltype Complex{Int64}"
    else
        @test str == "StructArray(::Vector{Int64}, ::Vector{Int64}) with eltype Complex{Int64}"
    end
    io = IOBuffer()
    Base.showarg(io, s, false)
    str = String(take!(io))
    if VERSION < v"1.6-"
        @test str == "StructArray(::Array{Int64,1}, ::Array{Int64,1})"
    else
        @test str == "StructArray(::Vector{Int64}, ::Vector{Int64})"
    end
end

@testset "append!!" begin
    dest_examples = [
        ("mutate", StructVector(a = [1], b = [2])),
        ("widen", StructVector(a = [1], b = [nothing])),
    ]
    itr = [(a = 1, b = 2), (a = 1, b = 2), (a = 1, b = 12)]
    itr_examples = [
        ("HasLength", () -> itr),
        ("StructArray", () -> StructArray(itr)),
        ("SizeUnknown", () -> (x for x in itr if isodd(x.a))),
        # Broken due to https://github.com/JuliaArrays/StructArrays.jl/issues/100:
        # ("empty", (x for x in itr if false)),
        ("stateful", () -> Iterators.Stateful(itr)),
    ]
    @testset "$destlabel $itrlabel" for (destlabel, dest) in dest_examples,
                                        (itrlabel, makeitr) in itr_examples

        @test vcat(dest, StructVector(makeitr())) == append!!(copy(dest), makeitr())
    end
end

struct ArrayConverter end

Adapt.adapt_storage(::ArrayConverter, xs::AbstractArray) = convert(Array, xs)

@testset "adapt" begin
    s = StructArray(a = 1:10, b = StructArray(c = 1:10, d = 1:10))
    t = adapt(ArrayConverter(), s)
    @test propertynames(t) == (:a, :b)
    @test s == t
    @test t.a isa Array
    @test t.b.c isa Array
    @test t.b.d isa Array
end

struct MyArray{T,N} <: AbstractArray{T,N}
    A::Array{T,N}
end
MyArray{T}(::UndefInitializer, sz::Dims) where T = MyArray(Array{T}(undef, sz))
Base.IndexStyle(::Type{<:MyArray}) = IndexLinear()
Base.getindex(A::MyArray, i::Int) = A.A[i]
Base.setindex!(A::MyArray, val, i::Int) = A.A[i] = val
Base.size(A::MyArray) = Base.size(A.A)
Base.BroadcastStyle(::Type{<:MyArray}) = Broadcast.ArrayStyle{MyArray}()
Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{MyArray}}, ::Type{ElType}) where ElType =
    MyArray{ElType}(undef, size(bc))

@testset "broadcast" begin
    s = StructArray{ComplexF64}((rand(2,2), rand(2,2)))
    @test isa(@inferred(s .+ s), StructArray)
    @test (s .+ s).re == 2*s.re
    @test (s .+ s).im == 2*s.im
    @test isa(@inferred(broadcast(t->1, s)), Array)
    @test all(x->x==1, broadcast(t->1, s))
    @test isa(@inferred(s .+ 1), StructArray)
    @test s .+ 1 == StructArray{ComplexF64}((s.re .+ 1, s.im))
    r = rand(2,2)
    @test isa(@inferred(s .+ r), StructArray)
    @test s .+ r == StructArray{ComplexF64}((s.re .+ r, s.im))

    # used inside of broadcast but we also test it here explicitly
    @test isa(@inferred(Base.dataids(s)), NTuple{N, UInt} where {N})

    s = StructArray{ComplexF64}((MyArray(rand(2,2)), MyArray(rand(2,2))))
    @test_throws MethodError s .+ s
end

@testset "staticarrays" begin
    # test that staticschema returns the right things
    for StaticVectorType = [SVector, MVector, SizedVector]
        @test StructArrays.staticschema(StaticVectorType{2,Float64}) == Tuple{Float64,Float64}
    end

    # test broadcast + components for vectors
    for StaticVectorType = [SVector, MVector, SizedVector]
        x = @inferred StructArray([StaticVectorType{2}(Float64[i;i+1]) for i = 1:2])
        y = @inferred StructArray([StaticVectorType{2}(Float64[i+1;i+2]) for i = 1:2])
        @test StructArrays.components(x) == ([1.0,2.0], [2.0,3.0])
        @test x .+ y == StructArray([StaticVectorType{2}(Float64[2*i+1;2*i+3]) for i = 1:2])
    end
    # test broadcast + components for general arrays
    for StaticArrayType = [SArray, MArray, SizedArray]
        x = @inferred StructArray([StaticArrayType{Tuple{1,2}}(ones(1,2) .+ i) for i = 0:1])
        y = @inferred StructArray([StaticArrayType{Tuple{1,2}}(2*ones(1,2) .+ i) for i = 0:1])
        @test StructArrays.components(x) == ([1., 2.], [1., 2.])
        @test x .+ y == StructArray([StaticArrayType{Tuple{1,2}}(3*ones(1,2) .+ 2*i) for i = 0:1])
    end

    # test FieldVector constructor (see https://github.com/JuliaArrays/StructArrays.jl/issues/205)
    struct FlippedVec2D <: FieldVector{2,Float64}
        x::Float64
        y::Float64
    end
    # tuple constructors should respect the flipped ordering
    FlippedVec2D(t::Tuple) = FlippedVec2D(t[2], t[1])

    # define a custom getindex to test StructArrays.component(::FieldArray) behavior
    Base.getindex(a::FlippedVec2D, index::Int) = index==1 ? a.y : a.x
    Base.Tuple(a::FlippedVec2D) = (a.y, a.x)
    a = StructArray([FlippedVec2D(1.0,2.0)])
    @test a.x == [1.0]
    @test a.y == [2.0]
    @test a.x[1] == a[1].x

    # test custom indices and components
    @test typeof(StructArrays.components(a)) == NamedTuple{(:x, :y), NTuple{2, Vector{Float64}}}
    @test StructArrays.components(a) == (; x = [1.0], y = [2.0])

    # test type stability of creating views with "many" homogeneous components
    for n in 1:10
        u = StructArray(randn(SVector{n, Float64}) for _ in 1:10, _ in 1:5)
        @inferred view(u, :, 1)
        @inferred view(u, 1, :)
    end
end


text:
using Bessels
using Test
import SpecialFunctions

@time @testset "besseli" begin include("besseli_test.jl") end
@time @testset "besselk" begin include("besselk_test.jl") end
@time @testset "besselj" begin include("besselj_test.jl") end
@time @testset "bessely" begin include("bessely_test.jl") end

ext: jl


text:
################################
######     Fixations      ######
################################

# E[Dn]  = LT(E[Dn+] + E[Dn-] + E[Dns])

# Neutral fixation corrected by background selection
"""

 - fixNeut()

Expected neutral fixations rate reduce by B value.

```math
\\mathbb{E}[D_{s}] = (1 - p_{-} - p_{+}) B \\frac{1}{2N}
```
# Returns
 - `Float64`: expected rate of neutral fixations.

"""
function fixNeut(param::parameters)
	# Synonymous probabilty * (fixation probabilty corrected by BGS value)
	out::Float64 = 0.25*(1.0/(param.B*param.NN))
	return out
end

# Negative fixations corrected by background selection
"""

	fixNegB(ppos)

Expected fixation rate from negative DFE.

```math
\\mathbb{E}[D_{n-}] =  p_{-}\\left(2^-\\alpha\\beta^\\alpha\\left(-\\zeta[\\alpha,\\frac{2+\\beta}{2}] + \\zeta[\\alpha,1/2(2-\\frac{1}{N+\\beta})]\\right)\\right)
```

# Arguments
 - `ppos::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of fixations from negative DFE.

"""
function fixNegB(param::parameters,ppos::Float64)
	# Non-synonymous proportion * negative alleles probability * fixation probability from gamma distribution
	out::Float64 = 0.75*(1-ppos)*(2^(-param.al))*(param.B^(-param.al))*(param.be^param.al)*(-zeta(param.al,1.0+param.be/(2.0*param.B))+zeta(param.al,0.5*(2-1.0/(param.N*param.B)+param.be/param.B)))
	return out
end

# Positive fixations
"""
	pFix()

Expected positive fixation rate.

```math
\\mathbb{E}[D_{n+}] =  p_{+} \\cdot B \\cdot (1 - e^{(-2s)})
```

# Arguments
 - `ppos::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of positive fixation.

"""
function pFix(param::parameters,gam::Int64)

	# Fixation probability
	s::Float64    = gam/(param.NN)
	pfix::Float64 = (1.0-ℯ^(-2.0*s))/(1.0-ℯ^(-2.0*gam))

	# Correcting pFix for large s following Uricchio et al. 2014
	if s >= 0.1
		pfix = ℯ^(-(1.0+s))
		lim = 0
		while(lim < 200)
			pfix = ℯ^((1.0+s)*(pfix-1.0))
			lim  = lim + 1
		end
		pfix = 1 -pfix
	end

	return pfix
end

# Positive fixations after apply Φ. reduction of positive fixations due deleterious linkage given a value B of background selection
"""

	fixPosSim(gamma,ppos)

Expected positive fixations rate reduced due to the impact of background selection and linkage. The probabilty of fixation of positively selected alleles is reduced by a factor Φ across all deleterious linked sites [`Analytical.phiReduction`](@ref).

```math
\\mathbb{E}[D_{n+}] =  \\Phi \\cdot \\mathbb{E}[D_{n+}]
```

# Arguments
 - `ppos::Float64`: positive selected alleles probabilty.

# Returns
 - `Float64`: expected rate of positive fixations under background selection

"""
function fixPosSim(param::parameters,gamma::Int64,ppos::Float64)

	redPlus = phiReduction(param,gamma)

	# Non-synonymous * positive alleles probability * B reduction * fixation probility
	out::Float64 = 0.75*ppos*redPlus*pFix(param,gamma)
	return out
end


text:
#=
test_timeevolution:
- Julia version: 1.5.2
- Author: qcancam
- Date: 2021-01-03
=#

include("../src/MatrixProductStates.jl")

using .MatrixProductStates, SparseArrays, Arpack, PyCall, PyPlot

#np=pyimport("numpy")
#@pyimport matplotlib.pyplot as plt

g = 1.0; L = 5
T=8.0
N=1000.0
τ = T/N

function H_TFIM(g, L)
    H = Array{Complex{Float64}, 2}[]
    id = [1  0;
          0  1]
    σˣ = [0  1;
          1  0]
    σᶻ = [1  0;
          0 -1]
    for i in 1:L-2
        push!(H,-σᶻ⊗σᶻ-g*σˣ⊗ id)
    end
    #=L-1=#
    push!(H,-σᶻ⊗σᶻ-g*(σˣ⊗ id+ id⊗ σˣ))

    return H
end

function H_TFIM_MPO(g, L)
    id = [1  0;
          0  1]
    σˣ = [0  1;
          1  0]
    σᶻ = [1  0;
          0 -1]
    W_tnsr = zeros(Complex{Float64}, 3, 3, 2, 2)
    W_tnsr[1, 1, :, :] = id
    W_tnsr[2, 1, :, :] = -σᶻ*0
    W_tnsr[3, 1, :, :] = -σˣ
    W_tnsr[3, 2, :, :] = σᶻ*0
    W_tnsr[3, 3, :, :] = id

    return MPO(W_tnsr, L)
end

H = H_TFIM(g, L)
ψ = randn(MPS{L, Complex{Float64}}, 100, 2)

ϕ,time_list,expect_list= real_time_evolution(ψ, H[1],H[2],H[L-1],H_TFIM_MPO(g,L), T,N,50)
#ϕ= real_time_evolution(ψ, H[1],H[2],H[L-1],H_TFIM_MPO(g,L), T,N,50)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(time_list, expect_list, "-", c="b", ms=10)
ax.set_xlabel("t")
ax.set_ylabel("⟨-∑σˣ(t)⟩")
#ax.set_title("2D Line Plot Sample")
#ax.grid(true)
plt.show()


text:
global CRRao_rng = MersenneTwister()

"""
```julia
using StableRNGs
CRRao.set_rng(StableRNG(1234))
```
"""
function set_rng(x)
    global CRRao_rng = x
end


text:
using Distributed, SharedArrays, LinearAlgebra

using ForwardDiff
using Distributions

struct KSDTest <: KernelGoodnessOfFitTest
    kernel::Kernel
    m::Integer
    α::Float64
end

struct KSDResult <: KernelGoodnessOfFitResult
    stat::Float64

    "p-value of test"
    p_val
    "size of test"
    α # size of test
    "`:reject` or `:accept` to indicate rejection or acceptance of H₀"
    result::Symbol # :reject or :accept
end

u_p(k::Kernel, x::AbstractVector, y::AbstractVector, s_x::AbstractVector, s_y::AbstractVector) = begin
    kernel(k, x, y) * dot(s_x, s_y) +
        dot(s_x, k_dy(k, x, y)) +
        dot(k_dx(k, x, y), s_y) +
        tr(k_dxdy(k, x, y))
end

u_p(k::Kernel, p::Distribution, x::AbstractVector, y::AbstractVector) = begin
    s_x = gradlogpdf(p, x)
    s_y = gradlogpdf(p, y)

    u_p(k, x, y, s_x, s_y)
end

# compute unbiased estimate
Sᵤ_estimate(samples::AbstractArray, f) = begin
    n = size(samples, 2)
    res = 0.0

    for i = 1:n
        for j = 1:n
            if i == j
                continue
            end
            r = f(samples[:, i], samples[:, j])
            res += (r / (n * (n - 1)))
        end
    end

    res
end

"Performs a goodness-of-fit test based on the Kernelized Stein Discrepancy (KSD)."
perform(t::KSDTest, q, xs) = begin
    bootstrapped = SharedArray{Float64}(t.m);

    f = (x, y) -> u_p(t.kernel, q, x, y)

    @sync @distributed for b = 1:t.m
        bootstrapped[b] = bootstrap_degenerate(xs, f)

        if (b % 10000) == 0
            println("[$b / $m] Done!")
        end
    end

    # Unbiased estimate
    # println("Computing unbiased estimate")
    Ŝᵤ = Sᵤ_estimate(xs, f)

    # compute average times this occurs
    pval = reduce(+, bootstrapped .> Ŝᵤ) / t.m

    res = :accept
    if pval < t.α
        res = :reject;
    end

    # TODO: return some result containing information
    KSDResult(Ŝᵤ, pval, t.α, res)
end


text:
abstract QuExponential <: QuPropagatorMethod

@doc """
Exponential solver, using Expokit.expmv

Step Propagation using the exponential solver Expokit.expmv.

### Fields :

* options :: Dict

Dictionary to set the size of Krylov subspace and tolerance by using
keys as `:m` and `:tol`.
""" ->
immutable QuExpokit <: QuExponential
    options::Dict{Symbol, Any}
end

QuExpokit() = QuExpokit(Dict())

@doc """
Exponential solver, using ExpmV.expmv

Step Propagation using the exponential solver ExpmV.expmv.

### Fields :

* options :: Dict

Dictionary to set M, precision, shift, full_term by using
keys as `:M`, `:precision`, `:shift`, `:full_term`
""" ->
immutable QuExpmV <: QuExponential
    options::Dict{Symbol, Any}
end

QuExpmV() = QuExpmV(Dict())

function propagate(prob::QuExpokit, eq::QuEquation, t, current_t, current_qustate)
    dt = t - current_t
    dims = size(current_qustate)
    next_state = Expokit.expmv(dt, -im*coeffs(operator(eq, t)), coeffs(vec(current_qustate)), m = get(prob.options, :m, 30), tol = get(prob.options, :tol, 1e-7))
    CQST = QuBase.similar_type(current_qustate)
    return CQST(reshape(next_state, dims), bases(current_qustate))
end

function propagate(prob::QuExpmV, eq::QuEquation, t, current_t, current_qustate)
    dt = t - current_t
    dims = size(current_qustate)
    next_state = ExpmV.expmv(dt, -im*coeffs(operator(eq, t)), coeffs(vec(current_qustate)), M = get(prob.options, :M, []), prec = get(prob.options, :prec, "double"),
                            shift = get(prob.options, :shift, false), full_term = get(prob.options, :full_term, false))
    CQST = QuBase.similar_type(current_qustate)
    return CQST(reshape(next_state, dims), bases(current_qustate))
end

export QuExpokit,
      QuExpmV


text:
# ## Buckling column

# ## Model

module buckling_column   #jl

using ReachabilityAnalysis

@taylorize function buckling_column!(dx, x, params, t)
    dx[1] = x[2]
    dx[2] = 2 * x[1] - x[1]^3 - 0.2 * x[2] + 0.1
    return dx
end

function model(X0)                                  #jl
    S = @system(x' = buckling_column!(x), dim:2)    #jl
    return IVP(S, X0)                               #jl
end                                                 #jl

end  #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#####
##### Diffusive fluxes
#####

@inline diffusive_flux_x(i, j, k, grid, clock, κᶠᶜᶜ::Number, c) = κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, c)
@inline diffusive_flux_y(i, j, k, grid, clock, κᶜᶠᶜ::Number, c) = κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, c)
@inline diffusive_flux_z(i, j, k, grid, clock, κᶜᶜᶠ::Number, c) = κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, c)

# Diffusivities-as-functions

@inline diffusive_flux_x(i, j, k, grid, clock, κ::Function, c) =
        diffusive_flux_x(i, j, k, grid, clock, κ(xnode(Face, i, grid), ynode(Center, j, grid), znode(Center, k, grid), clock.time), c)

@inline diffusive_flux_y(i, j, k, grid, clock, κ::Function, c) =
        diffusive_flux_y(i, j, k, grid, clock, κ(xnode(Center, i, grid), ynode(Face, j, grid), znode(Center, k, grid), clock.time), c)

@inline diffusive_flux_z(i, j, k, grid, clock, κ::Function, c) =
        diffusive_flux_z(i, j, k, grid, clock, κ(xnode(Center, i, grid), ynode(Center, j, grid), znode(Face, k, grid), clock.time), c)

#####
##### Laplacian diffusion operator
#####

"""
    ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, κˣ, κʸ, κᶻ, c)

Calculates diffusion for a tracer c via

    1/V * [δxᶜᵃᵃ(κˣ * Ax * ∂xᶠᵃᵃ(c)) + δyᵃᶜᵃ(κʸ * Ay * ∂yᵃᶠᵃ(c)) + δzᵃᵃᶜ(κᶻ * Az * ∂zᵃᵃᶠ(c))]

which will end up at the location `ccc`.
"""
@inline function ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, clock, κˣ, κʸ, κᶻ, c)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, clock, κˣ, c) +
                                    δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, clock, κʸ, c) +
                                    δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, clock, κᶻ, c))
end

@inline ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, clock, κ, c) = ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, clock, κ, κ, κ, c)


text:
struct ExtendedJumpArray{T3<:Number, T1, T<:AbstractArray{T3,T1},T2} <: AbstractArray{T3,1}
  u::T
  jump_u::T2
end

Base.length(A::ExtendedJumpArray) = length(A.u) + length(A.jump_u)
Base.size(A::ExtendedJumpArray) = (length(A),)
@inline function Base.getindex(A::ExtendedJumpArray,i::Int)
  i <= length(A.u) ? A.u[i] : A.jump_u[i-length(A.u)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I::Int...)
  prod(I) <= length(A.u) ? A.u[I...] : A.jump_u[prod(I)-length(A.u)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I::CartesianIndex{1})
  A[I[1]]
end
@inline Base.setindex!(A::ExtendedJumpArray,v,I...) = (A[CartesianIndices(A.u,I...)] = v)
@inline Base.setindex!(A::ExtendedJumpArray,v,I::CartesianIndex{1}) = (A[I[1]] = v)
@inline function Base.setindex!(A::ExtendedJumpArray,v,i::Int)
  i <= length(A.u) ? (A.u[i] = v) : (A.jump_u[i-length(A.u)] = v)
end

Base.IndexStyle(::Type{<:ExtendedJumpArray}) = IndexLinear()
Base.similar(A::ExtendedJumpArray) = ExtendedJumpArray(similar(A.u),similar(A.jump_u))
Base.similar(A::ExtendedJumpArray,::Type{S}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))
Base.zero(A::ExtendedJumpArray) = fill!(similar(A),0)

# Required for non-diagonal noise
LinearAlgebra.mul!(c::ExtendedJumpArray,A::AbstractVecOrMat,u::AbstractVector) = mul!(c.u,A,u)

# Ignore axes
Base.similar(A::ExtendedJumpArray,::Type{S},axes::Tuple{Base.OneTo{Int}}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))

# Stiff ODE solver
function ArrayInterface.zeromatrix(A::ExtendedJumpArray)
  u = [vec(A.u);vec(A.jump_u)]
  u .* u' .* false
end
LinearAlgebra.ldiv!(A,b::ExtendedJumpArray) = LinearAlgebra.ldiv!(A,[vec(b.u);vec(b.jump_u)])

function recursivecopy!(dest::T, src::T) where T<:ExtendedJumpArray
  recursivecopy!(dest.u,src.u)
  recursivecopy!(dest.jump_u,src.jump_u)
end
Base.show(io::IO,A::ExtendedJumpArray) = show(io,A.u)
TreeViews.hastreeview(x::ExtendedJumpArray) = true
plot_indices(A::ExtendedJumpArray) = eachindex(A.u)

## broadcasting

struct ExtendedJumpArrayStyle{Style <: Broadcast.BroadcastStyle} <: Broadcast.AbstractArrayStyle{Any} end
ExtendedJumpArrayStyle(::S) where {S} = ExtendedJumpArrayStyle{S}()
ExtendedJumpArrayStyle(::S, ::Val{N}) where {S,N} = ExtendedJumpArrayStyle(S(Val(N)))
ExtendedJumpArrayStyle(::Val{N}) where N = ExtendedJumpArrayStyle{Broadcast.DefaultArrayStyle{N}}()

# promotion rules
@inline function Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle{AStyle}, ::ExtendedJumpArrayStyle{BStyle}) where {AStyle, BStyle}
    ExtendedJumpArrayStyle(Broadcast.BroadcastStyle(AStyle(), BStyle()))
end
Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle{Style}, ::Broadcast.DefaultArrayStyle{0}) where Style<:Broadcast.BroadcastStyle = ExtendedJumpArrayStyle{Style}()
Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle, ::Broadcast.DefaultArrayStyle{N}) where N = Broadcast.DefaultArrayStyle{N}()

combine_styles(args::Tuple{})         = Broadcast.DefaultArrayStyle{0}()
@inline combine_styles(args::Tuple{Any})      = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]))
@inline combine_styles(args::Tuple{Any, Any}) = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]), Broadcast.BroadcastStyle(args[2]))
@inline combine_styles(args::Tuple)   = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]), combine_styles(Base.tail(args)))

function Broadcast.BroadcastStyle(::Type{ExtendedJumpArray{T,S}}) where {T, S}
    ExtendedJumpArrayStyle(Broadcast.result_style(Broadcast.BroadcastStyle(T)))
end

@inline function Base.copy(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}) where Style
    ExtendedJumpArray(copy(unpack(bc, Val(:u))),copy(unpack(bc, Val(:jump_u))))
end

@inline function Base.copyto!(dest::ExtendedJumpArray, bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}) where Style
    copyto!(dest.u,unpack(bc, Val(:u)))
    copyto!(dest.jump_u,unpack(bc, Val(:jump_u)))
    dest
end

# drop axes because it is easier to recompute
@inline unpack(bc::Broadcast.Broadcasted{Style}, i) where Style = Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
@inline unpack(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}, i) where Style = Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
unpack(x,::Any) = x
unpack(x::ExtendedJumpArray, ::Val{:u}) = x.u
unpack(x::ExtendedJumpArray, ::Val{:jump_u}) = x.jump_u

@inline unpack_args(i, args::Tuple) = (unpack(args[1], i), unpack_args(i, Base.tail(args))...)
unpack_args(i, args::Tuple{Any}) = (unpack(args[1], i),)
unpack_args(::Any, args::Tuple{}) = ()


text:
using Printf

function main_pre(n_x,n_y,s_x,e_x,s_y,e_y)
    xy_or_r = 1          # 1:x,y 2:r,t
    xnum = n_x           #       2:rnum
    ynum = n_y           #       2:tnum
    lenx = abs(e_x-s_x)           #       2:inr
    leny = abs(e_y-s_y)           #       2:outr
    st_ang  = -1/2*pi       # 2:angle
    en_ang  = -3/2*pi          # 2:angle

    outdir = "grid_morton"
    make_dir(outdir)
    result = "result_morton"
    make_dir(result)
    result = "post_result_morton"
    make_dir(result)
    nodes, xnum_max, ynum_max = mk_gird(xnum,ynum,lenx,leny,xy_or_r,s_x,s_y,st_ang,en_ang,outdir)
end

function mk_gird(xnum,ynum,lenx,leny,xy_or_r,s_x,s_y,st_ang,en_ang,outdir)
    """
    nodes[i,j,k]
    i   : x,r方向の番号
    j   : y,theta方向の番号
    k=1 : 点のx座標
    k=2 : 点のy座標

    nodes[1,:]      : x方向境界
    nodes[:,1]      : y方向境界
    nodes[xnum+2,:] : x方向境界
    nodes[:,ynum+2] : y方向境界
    """
    xnum_max = xnum+1+2
    ynum_max = ynum+1+2
    nodes = zeros(xnum_max, ynum_max, 3)
    if xy_or_r == 1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                x = lenx/(xnum)*(i-2) + s_x
                y = leny/(ynum)*(j-2) + s_y
                nodes[i,j,1] = x
                nodes[i,j,2] = y
            end
        end
    elseif xy_or_r == 2
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                r = lenx + (leny-lenx)/ynum*(j-2)
                theta = st_ang - (st_ang-en_ang)/xnum*(i-2)
                x = r * cos(theta)
                y = r * sin(theta)
                nodes[i,j,1]=x
                nodes[i,j,2]=y
            end
        end
    end

    #=
    仮想セルの作成：現在あるセル境界線を延長して作成
    そのまま延長した際にセルが交差するのを防ぐため、仮想セルの大きさを100^-1倍にする.
    ベクトルで考えれば下記のような計算になる．たぶん
    また、[1,1]等の角にはダミー値が入っている.
    =#
    for j in 1:ynum_max
        for k in 1:2
            nodes[1,j,k]        =  1.01*nodes[2,j,k] - 0.01*nodes[3,j,k]
            nodes[xnum_max,j,k] =  1.01*nodes[xnum_max-1,j,k] - 0.01*nodes[xnum_max-2,j,k]
        end
    end
    for i in 1:xnum_max
        for k in 1:2
            nodes[i,1,k]        =  1.01*nodes[i,2,k] - 0.01*nodes[i,3,k]
            nodes[i,ynum_max,k] =  1.01*nodes[i,ynum_max-1,k] - 0.01*nodes[i,ynum_max-2,k]
        end
    end

    fff=outdir*"/nodes"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        for i in 1:xnum_max
            for j in 1:ynum_max
                x = @sprintf("%8.8e", nodes[i,j,1])
                y = @sprintf("%8.8e", nodes[i,j,2])
                write(f,string(i)*" "*string(j)*" "*x*" "*y*"\n")
            end
        end
    end
    println("write "*fff)

    # nodes_forvtk
    nodes_num = zeros(Int,xnum_max, ynum_max)
    fff=outdir*"/nodes_forvtk"
    open(fff,"w") do f
        write(f,"nodes: xnum, ynum , x, y\n")
        k=1
        for i in 2:xnum_max-1
            for j in 2:ynum_max-1
                x = @sprintf("%8.8e", nodes[i,j,1])
                y = @sprintf("%8.8e", nodes[i,j,2])
                z = @sprintf("%8.8e", nodes[i,j,3])
                write(f,string(k)*" "*x*" "*y*" "*z*"\n")
                nodes_num[i,j] = k
                k = k+1
            end
        end
    end
    println("write "*fff)

    fff=outdir*"/nodesnum"
    open(fff,"w") do f
        write(f,"nodesnum: xnum_max, ynum_max\n")
        write(f,string(xnum_max)*" "*string(ynum_max)*"\n")
    end

    ### element ###
    fff=outdir*"/element_forvtk"
    open(fff,"w") do f
        write(f,"elements:cell_xnum, lup,rup,ldown,rdown \n")
        k=1
        for i in 2:xnum_max-2
            for j in 2:ynum_max-2
                d1 = @sprintf("%1.0f", nodes_num[i,j])
                d2 = @sprintf("%1.0f", nodes_num[i,j+1])
                d3 = @sprintf("%1.0f", nodes_num[i+1,j+1])
                d4 = @sprintf("%1.0f", nodes_num[i+1,j])
                write(f,string(k)*" "*d1*" "*d2*" "*d3*" "*d4*"\n")
                k = k+1
            end
        end
    end
    println("write "*fff)

    return  nodes,xnum_max,ynum_max
end

function make_dir(outdir)
    k = 0
    try rm(outdir,recursive=true)
    catch
        mkdir(outdir)
        k = 1
    end

    if k == 0
        mkdir(outdir)
    end
end

max_stars_repo_name: hide-dog/Qtree

max_issues_repo_name: hide-dog/Qtree

max_forks_repo_name: hide-dog/Qtree


text:
using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using HDF5

import PyPlot.plt
using LaTeXStrings

cmap = plt[:get_cmap]("plasma")

import YAML
config = YAML.load(open("config.yaml"))

grd = config["grid"]
grid = Grid(grd)

# Load the parameters file from the configuration setting, and make this an object.
pars = convert_dict(config["parameters"], config["model"])
vel = pars.vel # [km/s]
npix = config["npix"] # number of pixels
species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

# The list of inclinations we want to create
nincl = 10
incls = linspace(140., 158., nincl)

nchan = 10

# Read the parameters from the config file
pars = convert_dict(config["parameters"], config["model"])

# For *this purpose only*, read in the flagged data in addition to the unflagged data
# so that we can export a model for these otherwise flagged visibilities
dvarr = DataVis(config["data_file"], true)
# Do this as we do in `mach_three.jl`
for dset in dvarr
    # Conjugation is necessary for the SMA and ALMA
    visibilities.conj!(dset) # Swap UV convention
end

# Stick with the baselines contained in the first channel.
dv = dvarr[1]

for incl in incls
    pars.incl = incl
    image = @sprintf("image_%.0f.out", incl)
    im = imread(image)
    skim = imToSky(im, pars.dpc)
    corrfun!(skim) # alpha = 1.0

    # Determine dRA and dDEC from the image and distance
    dRA = abs(skim.ra[2] - skim.ra[1])/2. # [arcsec] the half-size of a pixel
    println("dRA is ", dRA, " arcsec")

    mvarr = Array(DataVis, nchan)

    for i=1:nchan

        # FFT the appropriate image channel
        vis_fft = transform(skim, i)

        # Interpolate the `vis_fft` to the same locations as the DataSet
        mvis = ModelVis(dv, vis_fft)

        # Apply the phase correction here, since there are fewer data points
        phase_shift!(mvis, pars.mu_RA + dRA, pars.mu_DEC - dRA)

        dvis = visibilities.ModelVis2DataVis(mvis)
        mvarr[i] = dvis

        # Now swap the model and residual visibilities back to ALMA/SMA convetion
        visibilities.conj!(mvarr[i])

    end

    model_file = @sprintf("model_%.0f.hdf5", incl)

    visibilities.write(mvarr, model_file)

end

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
include("datasstructures/rsparsebitset.jl")

include("absolute.jl")
include("disjunctive.jl")
include("alldifferent.jl")
include("compacttable.jl")
include("equal.jl")
include("notequal.jl")
include("lessorequal.jl")
include("greaterorequal.jl")
include("interval.jl")
include("sumtozero.jl")
include("sumlessthan.jl")
include("sumgreaterthan.jl")
include("islessorequal.jl")
include("inset.jl")
include("reifiedinset.jl")
include("binaryor.jl")
include("isbinaryor.jl")
include("setdiffsingleton.jl")
include("element2d.jl")
include("element1d.jl")
include("binarymaximum.jl")
include("setequalconstant.jl")

"""
    addOnDomainChange!(x::AbstractIntVar, constraint::Constraint)

Make sure `constraint` will be propagated if `x`'s domain changes.
"""
function addOnDomainChange!(x::Union{IntVar, BoolVar, IntSetVar}, constraint::Constraint)
    if !(constraint in x.onDomainChange)
        push!(x.onDomainChange, constraint)
    end
end

addOnDomainChange!(x::Union{IntVarView, BoolVarView}, constraint::Constraint) = addOnDomainChange!(x.x, constraint)

"""
    addToPropagate!(toPropagate::Set{Constraint}, constraints::Array{Constraint})

Add the constraints to `toPropagate` only if they are active.
"""
function addToPropagate!(toPropagate::Set{<:Constraint}, constraints::Array{<:Constraint})
    for constraint in constraints
        if constraint.active.value
            push!(toPropagate, constraint)
        end
    end
end

"""
    triggerDomainChange!(toPropagate::Set{Constraint}, x::AbstractIntVar)

Add the constraints that have to be propagated when the domain of `x` changes to `toPropagate`.
"""
triggerDomainChange!(toPropagate::Set{<:Constraint}, x::Union{AbstractIntVar, AbstractBoolVar, IntSetVar}) = addToPropagate!(toPropagate, getOnDomainChange(x))
getOnDomainChange(x::Union{IntVar, BoolVar, IntSetVar}) = x.onDomainChange
getOnDomainChange(x::Union{IntVarView, BoolVarView}) = getOnDomainChange(x.x)


text:
function AbstractConstantThrustPropulsionModel(arg0::jdouble, arg1::jdouble, arg2::Vector3D, arg3::JString)
    return AbstractConstantThrustPropulsionModel((jdouble, jdouble, Vector3D, JString), arg0, arg1, arg2, arg3)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_acceleration(obj::ThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::FieldAttitude, arg2::Vector{RealFieldElement})
    return jcall(obj, "getAcceleration", FieldVector3D, (FieldSpacecraftState, FieldAttitude, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function get_acceleration(obj::ThrustPropulsionModel, arg0::SpacecraftState, arg1::Attitude, arg2::Vector{jdouble})
    return jcall(obj, "getAcceleration", Vector3D, (SpacecraftState, Attitude, Vector{jdouble}), arg0, arg1, arg2)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_direction(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getDirection", Vector3D, ())
end

function get_direction(obj::ThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getDirection", Vector3D, (SpacecraftState,), arg0)
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getFlowRate", jdouble, ())
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getFlowRate", RealFieldElement, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getFlowRate", jdouble, (SpacecraftState,), arg0)
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getFlowRate", jdouble, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel, arg0::Vector{RealFieldElement})
    return jcall(obj, "getFlowRate", RealFieldElement, (Vector{RealFieldElement},), arg0)
end

function get_flow_rate(obj::AbstractConstantThrustPropulsionModel, arg0::Vector{jdouble})
    return jcall(obj, "getFlowRate", jdouble, (Vector{jdouble},), arg0)
end

function get_isp(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getIsp", jdouble, ())
end

function get_isp(obj::ThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getIsp", jdouble, (SpacecraftState,), arg0)
end

function get_mass_derivatives(obj::ThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getMassDerivatives", RealFieldElement, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_mass_derivatives(obj::ThrustPropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getMassDerivatives", jdouble, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_name(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getName", JString, ())
end

function get_parameters_drivers(obj::PropulsionModel)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_thrust(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getThrust", jdouble, ())
end

function get_thrust(obj::ThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getThrust", jdouble, (SpacecraftState,), arg0)
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel)
    return jcall(obj, "getThrustVector", Vector3D, ())
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "getThrustVector", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel, arg0::SpacecraftState)
    return jcall(obj, "getThrustVector", Vector3D, (SpacecraftState,), arg0)
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "getThrustVector", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel, arg0::Vector{RealFieldElement})
    return jcall(obj, "getThrustVector", FieldVector3D, (Vector{RealFieldElement},), arg0)
end

function get_thrust_vector(obj::AbstractConstantThrustPropulsionModel, arg0::Vector{jdouble})
    return jcall(obj, "getThrustVector", Vector3D, (Vector{jdouble},), arg0)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::PropulsionModel, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
######### logistic regression on binary response  ###########

# simulate dataset
srand(1)
n = 1000
nbeta = 10 # number of predictors, including intercept
X = [ones(n) randn((n, nbeta-1))]
beta0 = randn((nbeta,))
Y = rand(n) .< ( 1 ./ (1. + exp(- X * beta0)))

# define model
ex = quote
	vars ~ Normal(0, 1.0)  # Normal prior, std 1.0 for predictors
	prob = 1 / (1. + exp(- X * vars))
	Y ~ Bernoulli(prob)
end

m = model(ex, vars=zeros(nbeta))
mg = model(ex, gradient=true, vars=zeros(nbeta))

name = "binomial 10x1000"

f1 = ()-> m.eval(m.init)
f2 = ()-> mg.evalallg(mg.init)
f3 = ()-> run(m * RWM(0.1), steps=100)

res = [	benchmark( f1, "loglik eval", name, 1000) ;
		benchmark( f2, "loglik and gradient eval", name, 1000) ;
		benchmark( f3, "100 RWM steps", name, 10) ]


text:
"""
FunctionBlocks
"""

export ConstantBlock

mutable struct ConstantBlock <: AbstractFunctionBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    value

    function ConstantBlock(value, ::Type{Tv} = Float64) where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}(), signal(value, Tv)) 0 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::ConstantBlock)
    println(io, "Constant($(blk.value))")
end

"""
to expr
"""

function _toexpr(blk::ConstantBlock)
    _toexpr(blk.value)
end


text:
module DataAugmentation

using ColorBlendModes
using CoordinateTransformations
using ImageDraw
using Images
using Images: Colorant, permuteddimsview
using ImageTransformations
using ImageTransformations: center, _center, box_extrapolation, warp!
using Interpolations
using LinearAlgebra: I
using Parameters
using Rotations
using Setfield
using StaticArrays
using Statistics

include("./base.jl")
include("./wrapper.jl")
include("./buffered.jl")
include("./sequence.jl")
include("./utils/draw.jl")
include("./items/arrayitem.jl")
include("./projective/base.jl")
include("./items/image.jl")
include("./items/keypoints.jl")
include("./items/mask.jl")
include("./projective/bounds.jl")
include("./projective/compose.jl")
include("./projective/crop.jl")
include("./projective/affine.jl")
include("./preprocessing.jl")

export Item,
    Transform,
    ArrayItem,
    MapElem,
    Identity,
    Sequence,
    Project,
    Image,
    Keypoints,
    Polygon,
    ToEltype,
    ImageToTensor,
    Normalize,
    NormalizeIntensity,
    ToBinary,
    AddChannel,
    MaskMulti,
    MaskBinary,
    BoundingBox,
    ScaleKeepAspect,
    ScaleRatio,
    itemdata,
    Crop,
    CenterCrop,
    RandomCrop,
    ScaleFixed,
    Rotate,
    RandomResizeCrop,
    CenterResizeCrop,
    Buffered,
    BufferedThreadsafe,
    OneHot,
    apply,
    Reflect,
    FlipX,
    FlipY,
    PinOrigin,
    apply!,
    PadDivisible,
    ResizePadDivisible,
    onehot,
    showitem

end # module


text:
module LossFunctions

import Base.*
using Base.Cartesian
using Markdown
using SparseArrays
using InteractiveUtils
using RecipesBase
using LearnBase

import LearnBase:
    value, value!,
    deriv, deriv2, deriv!,
    value_deriv,
    scaled,
    isminimizable,
    isdifferentiable,
    istwicedifferentiable,
    isconvex,
    isstrictlyconvex,
    isstronglyconvex,
    isnemitski,
    isunivfishercons,
    isfishercons,
    islipschitzcont,
    islocallylipschitzcont,
    islipschitzcont_deriv, # maybe overkill
    isclipable,
    ismarginbased,
    isclasscalibrated,
    isdistancebased,
    issymmetric

export

    Loss,
    SupervisedLoss,
    MarginLoss,
    DistanceLoss,
    value,
    value!,
    deriv,
    deriv!,
    deriv2,
    deriv2!,

    ZeroOneLoss,
    LogitMarginLoss,
    PerceptronLoss,
    HingeLoss,
    L1HingeLoss,
    L2HingeLoss,
    SmoothedL1HingeLoss,
    ModifiedHuberLoss,
    L2MarginLoss,
    ExpLoss,
    SigmoidLoss,
    DWDMarginLoss,

    LPDistLoss,
    L1DistLoss,
    L2DistLoss,
    PeriodicLoss,
    HuberLoss,
    EpsilonInsLoss,
    L1EpsilonInsLoss,
    L2EpsilonInsLoss,
    LogitDistLoss,
    QuantileLoss,
    PinballLoss,

    MisclassLoss,
    PoissonLoss,
    LogitProbLoss,
    CrossEntropyLoss,
    ZeroOneLoss,

    OrdinalMarginLoss,
    OrdinalHingeLoss,

    weightedloss,

    AggMode

include("devutils.jl")
include("aggregatemode.jl")

include("supervised/supervised.jl")
include("supervised/sparse.jl")
include("supervised/distance.jl")
include("supervised/margin.jl")
include("supervised/scaledloss.jl")
include("supervised/weightedbinary.jl")
include("supervised/other.jl")
include("supervised/ordinal.jl")
include("supervised/io.jl")

end # module


text:
using FrameworkPacote
using NLPModels
using Test

@testset "Básicos" begin
  nlp = ADNLPModel(
    x -> (x[1] - 1)^2 + (x[2] - 2)^2 / 4,
    zeros(2)
  )
  output = gradiente(nlp)
  @test isapprox(output.solution, [1.0; 2.0], rtol=1e-4)
  @test output.objective < 1e-4
  @test output.dual_feas < 1e-4
  @test output.status == :first_order

  output = gradiente(nlp, max_iter = 0)
  @test output.status == :max_iter
  output = gradiente(nlp, max_time = 0)
  @test output.status == :max_time
end


text:
module Layouts

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export layout, page_container, page

register_normal_element("q__layout", context = @__MODULE__)
register_normal_element("q__page__container", context = @__MODULE__)
register_normal_element("q__page", context = @__MODULE__)

function layout(args...; kwargs...)
  q__layout(args...; kw(kwargs)...)
end

function page_container(args...; kwargs...)
  q__page__container(args...; kw(kwargs)...)
end

function page(args...; kwargs...)
  q__page(args...; kw(kwargs)...)
end

end