text:
using Test, AxisKeys

@testset "basics" begin

    R = KeyedArray(rand(1:99, 3,4), (['a', 'b', 'c'], 10:10:40))
    @test axiskeys(R) == (['a', 'b', 'c'], 10:10:40)
    @test axiskeys(R,5) === Base.OneTo(1)

    @test R('c', 40) == R[3, 4]
    @test R(:, 30) == R[:, 3]

    @test R('b') == R[2,:]
    @test axiskeys(R(20)) == (['a', 'b', 'c'],)

    @test AxisKeys.setkey!(R, 123, 'a', 10) == 123
    @test R[1,1] == 123

    @test R('a', 10, :) isa SubArray{Int,0} # trailing colon
    R('a', 10, :) .= 321
    @test R[1,1] == 321

    @test (R[2,:] .= 1:4) isa SubArray # dotview
    @test R[2,:] == KeyedArray(1:4, (10:10:40,))

    @test R[1:end] == R[:] == vec(R.data)
    @test R[1:end] isa Vector  # issue #48
    @test R[1:2, 1, 1] == R.data[1:2, 1, 1]
    @test axiskeys(R[:, [0.9,0.1,0.9,0.1] .> 0.5],2) == [10,30]
    @test ndims(R[R .> 0.5]) == 1
    newaxis = [CartesianIndex{0}()]
    @test axiskeys(R[[1,3],newaxis,:]) == (['a', 'c'], Base.OneTo(1), 10:10:40)

    @test axiskeys(R[[3,1,1], :]) == (['c','a','a'], 10:10:40) # repeated
    @test axiskeys(R(['c','a','a'], :)) == (['c','a','a'], 10:10:40)
    @test axiskeys(R(['c'])) == (['c'], 10:10:40) # issue #35, don't drop length 1 dim.

    @test_throws Exception R(:nope) # ideally ArgumentError
    @test_throws Exception R('z')   # ideally BoundsError?
    @test_throws Exception R(99)
    @test_throws Exception R('c', 99)
    @test_throws BoundsError axiskeys(R,0)
    @test_throws Exception KeyedArray(rand(3,4), (['a', 'b'], 10:10:40)) # keys wrong length

    C = wrapdims(rand(10), 'a':'j')
    @test C('a':'c') == C[1:3]
    @test C(Base.Fix2(<=,'c')) == C[1:3]
    @test_skip axiskeys(C(Base.Fix2(<=,'c')),1) == 'a':'c' # ('a':'z')[1:3] isa StepRangeLen

    @test axiskeys(KeyedArray(C, collect('a':'j'))) === ('a':'j',) # re-constructor
    @test axiskeys(KeyedArray(R, ('a':'c', Base.OneTo(4)))) === ('a':'c', 10:10:40)

    D = wrapdims(rand(2,10,3) .+ (1:10)'./10, ["cat", "dog"], 0:10:90, nothing)
    @test D("cat") == D[1,:,:]
    @test_throws Exception D(10) # ambiguous
    @test_throws Exception D("cat", 10) # too few
    @test_throws Exception D("cat", 10, -10) # out of bounds
    @test_throws Exception D(10, "cat", 3) # wrong order

    E = wrapdims(rand(2,3,4))
    @test axiskeys(E) == axes(E)

    F = wrapdims(rand(5), 'a':'z') # range needs adjusting
    @test axiskeys(F,1) == 'a':'e'
    @test_throws Exception wrapdims(rand(5), ['a','b','c'])
    @test_throws Exception KeyedArray(rand(5), ['a','b','c'])

    G = KeyedArray([[1,2], [3,4], [5,6]], 'α':'γ')  # array of arrays, issue 83
    @test G[end] == [5,6]
    @test axiskeys(G') == (1:1, 'α':'γ')
    @test begin G('α',:) .= Ref([7,8]); G[1] == [7,8] end  # trailing colon still OK here
end
@testset "selectors" begin

    V = wrapdims(rand(Int8, 11), 0:0.1:1)

    @test V(==(0.1)) == V[2:2]
    @test V(Near(0.12)) == V(0.1) == V[2]
    @test V(Interval(0.1, 0.3)) == V[2:4]

    @test V(Index[1]) == V[1]
    @test V(Index[2:3]) == V[2:3]
    @test V(Index[end]) == V[end]
    # @static if VERSION >= v"1.4"
    #     @test V(Index[begin]) == V[1] # syntax error on 1.0
    # end

    V2 = wrapdims(rand(Int8, 5), [1,2,3,2,1])
    @test V2(==(2)) == V2[[2,4]]
    @test V2(==(2.0)) == V2[[2,4]]
    @test V2(Near(2.3)) == V2[2]
    @test V2(Interval(0.5, 1.5)) == V2[[1,5]]

    R = KeyedArray(rand(1:99, 3,4), (['a', 'b', 'c'], 10:10:40))
    @test R(==('a')) == R[1:1, :]
    @test R(Near(23)) == R[:, 2]
    @test R(Near(23.5)) == R[:, 2] # promotes to Real & then matches
    @test R(Interval(17,23)) == R[:, 2:2]
    @test R(Base.Fix2(<=,23)) == R[:, 1:2]
    @test axiskeys(R(Base.Fix2(<=,23)), 2) isa AbstractRange

    @test_throws BoundsError V(Index[99])

    # faster Near using searchsortedfirst
    V3 = wrapdims(parent(V), collect(axiskeys(V,1)))
    xs = 0.3:0.001:0.5
    @test [V(Near(x)) for x in xs] == [V3(Near(x)) for x in xs]
    # ... and with decreasing keys:
    V4 = wrapdims(rand(Int8,10), 1:-0.1:0.1)
    V5 = wrapdims(parent(V4), collect(axiskeys(V4,1)))
    @test_skip [V4[Near(x)] for x in xs] == [V5[Near(x)] for x in xs]  # not broken on 1.3

end
@testset "reverse selectors" begin # https://github.com/mcabbott/AxisKeys.jl/pull/5

    V = wrapdims(rand(Int8, 11), 0:0.1:1)

    @test V[Key(0.1)] == V(0.1) == V[2]
    @test V[==(0.1)] == V[2:2]
    @test V[Interval(0.1, 0.3)] == V[2:4]

    newaxis = [CartesianIndex{0}()]
    @test axiskeys(V[newaxis, Key(0.1), newaxis]) === (Base.OneTo(1), Base.OneTo(1))

    A = wrapdims(rand(1:99, 5,7,3), alp='a':'e', num=21:27, sym=[:x, :y, :z])

    r5 = rand(5) .> 0.5
    @test A[r5, ==(21), Key(:x)] == A[r5, 1:1, 1]
    @test axiskeys(A[CartesianIndex(1,1), ==(:z), 1,1,1]) == ([:z],)
    @test dimnames(A[:,newaxis,Interval(25,99),Key(:y)]) == (:alp, :_, :num)

end
@testset "names" begin

    # constructor
    @test dimnames(wrapdims(rand(3), :a)) == (:a,)
    @test dimnames(wrapdims(rand(3), b=1:3), 1) == :b
    @test axiskeys(wrapdims(rand(3), b=1:3)) == (1:3,)
    @test axiskeys(wrapdims(rand(3), b=1:3), :b) == 1:3

    # @test namedaxiskeys(wrapdims(rand(3), b=1:3)) === (b = 1:3,)
    # @test namedaxes(wrapdims(rand(3), b=1:3)) === (b = Base.OneTo(3),)

    # internal functions of NamedDims
    @test NamedDims.dimnames(rand(2,2)) === (:_, :_)
    @test NamedDims.dim((:a, :b, :c), :b) == 2
    @test NamedDims.order_named_inds(Val((:a, :b, :c)); a=1, c=2:3) === (1, Colon(), 2:3)
    @test 0 == @allocated NamedDims.order_named_inds(Val((:a, :b, :c)); a=1, c=2:3)

    # indexing etc, of commutative wrappers
    data = rand(1:99, 3,4)
    N1 = wrapdims(data, obs = ['a', 'b', 'c'], iter = 10:10:40)
    N2 = NamedDimsArray(KeyedArray(data, axiskeys(N1)), dimnames(N1))
    N3 = KeyedArray(NamedDimsArray(data, dimnames(N1)), axiskeys(N1))

    @testset "with $(typeof(N).name) outside" for N in [N2, N3]
        @test axiskeys(N) == (['a', 'b', 'c'], 10:10:40)
        @test axiskeys(N, :iter) == 10:10:40
        @test dimnames(N) == (:obs, :iter)
        @test axes(N, :obs) == 1:3
        @test size(N, :obs) == 3

        @test propertynames(N) == (:obs, :iter)
        @test N.obs == ['a', 'b', 'c']

        @test N(obs='a', iter=40) == N[obs=1, iter=4]
        @test N(obs='a') == N('a') == N[1,:] == N[obs=1]
        @test N(obs='a') == N('a') == view(N, 1,:) == view(N, obs=1)

        @test dimnames(N(obs='a')) == (:iter,)
        @test axiskeys(N(obs='b')) == (10:10:40,)

        @test_throws Exception N(obs=55)  # ideally ArgumentError
        @test_throws Exception N(obs='z') # ideally BoundsError
    end

    @testset "named_axiskeys" begin
        arr = wrapdims(randn(2,3), a=1:2, b='a':'c')
        @inferred named_axiskeys(arr)
        @test named_axiskeys(arr) === (a=1:2, b='a':'c')
        @test Tuple(named_axiskeys(arr)) === axiskeys(arr)

        @test named_axiskeys(rename(arr, :a=>:x)) == (x=1:2, b='a':'c')
        @test named_axiskeys(rename(arr, (:y, :z))) == (y=1:2, z='a':'c')

        nonames = KeyedArray(randn(2,3), (1:2, 'a':'c'))
        @test_throws ErrorException named_axiskeys(nonames)

        v = wrapdims(randn(3), x=1:3)
        @test_throws ErrorException named_axiskeys(v .+ v') # duplicate names
    end

end
@testset "broadcasting" begin
    using Base: OneTo
    @testset "keys" begin

        V = wrapdims(rand(Int8, 11), 0:0.1:1)
        @test axiskeys(V .+ 10) == axiskeys(V)
        @test axiskeys(V .+ V') == (axiskeys(V,1), axiskeys(V,1))
        @test axiskeys(V .+ rand(3)') === (axiskeys(V,1), OneTo(3))
        @test axiskeys(V .+ rand(11)) == (axiskeys(V,1),)

        M = wrapdims(rand(Int8, 3,11), [:a, :b, :c], 0:0.1:1)
        @test axiskeys(M .+ V') == axiskeys(M)

        W = wrapdims(rand(11), 0.1:0.1:1.1)
        @test axiskeys(V .+ W') == (axiskeys(V,1), axiskeys(W,1))
        @test_throws Exception axiskeys(V .+ W)

    end
    @testset "with names" begin

        vec_x = wrapdims(ones(2), :x)
        vec_y = wrapdims(ones(2), y='α':'β')
        mat_r = wrapdims(ones(2,2), 11:12, 'α':'β')
        mat_y = wrapdims(ones(2,2), :_, :y)

        @test axiskeys(vec_x .+ mat_r .+ mat_y) == axiskeys(mat_r)
        @test dimnames(vec_x .+ mat_r .+ mat_y) == (:x, :y)

        @test axiskeys(sqrt.(vec_x .+ vec_y') ./ mat_r) == axiskeys(mat_r)
        @test dimnames(sqrt.(vec_x .+ vec_y') ./ mat_r) == (:x, :y)

        yy = vec_y .+ mat_y
        @test axiskeys(yy' .+ mat_r) == (11:12, 'α':'β')

        @test_throws Exception vec_x .+ vec_y

    end
    @testset "in-place" begin

        v1 = wrapdims(ones(2), ["a", "b"])
        v2 = wrapdims(ones(2), :μ)
        v3 = wrapdims(ones(2), [11, 22])
        z = zeros(2,2)

        @test dimnames(v1 .= v1 .+ v2) == (:μ,)
        @test v1[1] == 2

        @test dimnames(v2 .= v3 .+ 5) == (:μ,)
        @test v2[1] == 6

        @test axiskeys(z .= v1 .+ v2') == (["a", "b"], Base.OneTo(2))

        @test_throws Exception v3 .= v1 .+ v2
        @test_throws Exception zeros(2) .= v1 .+ v3

    end
    @testset "unify rules" begin

        using AxisKeys: who_wins
        @test who_wins(1:2, [1,2]) === 1:2
        @test who_wins(1:2, 1.0:2.0) === 1:2
        @test who_wins(OneTo(2), [3,4]) == [3,4]
        @test who_wins(1:2, [3,4]) === nothing

        using AxisKeys: unify_keys, unify_longest
        @test unify_keys((OneTo(2), 1:2), (3:4, [1,2])) === (3:4, 1:2)
        @test unify_longest((OneTo(2), 1:2), (3:4, [1,2], [0,1])) == (3:4, 1:2, [0,1])
        @test_throws Exception unify_keys((1:2,), ([3,4],))

    end
end
@testset "bitarray" begin

    x = wrapdims(rand(3,4), a=11:13, b=21:24.0)
    b = rand(12) .> 0.5
    @test length(x[b]) == sum(b)
    m = rand(3,4) .> 0.5 # BitArray{2}
    @test size(x[m]) == (sum(m),)

    # indexing a view, https://github.com/JuliaArrays/AxisArrays.jl/issues/179
    v = view(x, :, 1:2)
    b = rand(6) .> 0.5 # BitArray{1}
    @test length(v[b]) == sum(b)
    @test length(view(v, b)) == sum(b)
    m = rand(3,2) .> 0.5 # BitArray{2}
    @test size(v[m]) == (sum(m),)
    @test size(view(v, m)) == (sum(m),)

end
@testset "namedtuple" begin

    @test keys(NamedTuple(wrapdims([1,2,3], z='a':'c'))) == (:a, :b, :c)

    @test axiskeys(KeyedArray((a=1, b=2, c=3))) == ([:a, :b, :c],)
    @test axiskeys(wrapdims((a=1, b=2, c=3))) == ([:a, :b, :c],)
    @test dimnames(wrapdims((a=1, b=1+2, c=33), :z)) == (:z,)

end
@testset "mutation" begin

    V = wrapdims([3,5,7,11], μ=10:10:40)
    @test axiskeys(push!(V, 13)) == (10:10:50,)

    @test pop!(V) == 13
    @test axiskeys(V) == (10:10:40,)

    V2 = KeyedArray(rand(3), Base.OneTo(3))
    @test axiskeys(push!(V2, 0)) === (Base.OneTo(4),)
    @test axiskeys(append!(V2, [7,7])) === (Base.OneTo(6),)

    AxisKeys.nameouter() || # fails with nda(ka(...))
        @test axiskeys(append!(V2, V),1) == [1, 2, 3, 4, 5, 6, 10, 20, 30, 40]

    W = wrapdims([1,2,3], ["a", "b", "c"])
    push!(W, d=4)
    push!(W, "e" => 5)
    @test axiskeys(W,1) == ["a", "b", "c", "d", "e"]

end


text:
facts("spread constructors work") do

    context("long vertical") do

        p1 = LongPut(Ticker("SPX"), 1.50, 2000.00)
        p2 = ShortPut(Ticker("SPX"), 1.00, 2010.00)
        p3 = ShortPut(Ticker("SPX"), 2.00, 1990.00)

        c1 = LongCall(Ticker("SPX"), 1.50, 2000.00)
        c2 = ShortCall(Ticker("SPX"), 1.00, 2010.00)
        c3 = ShortCall(Ticker("SPX"), 2.00, 1990.00)

        @pending typeof(p1 + p2) => Spread{ShortPut, LongPut}
        @pending typeof(p1 + p3) => Spread{LongPut, ShortPut}
        @pending typeof(p1 + p2) => LongPutVertical
        @pending typeof(p1 + p3) => ShortPutVertical
    end

    context("short vertical") do

        p1 = LongPut(Ticker("SPX"), 1.50, 2000.00)
        p2 = ShortPut(Ticker("SPX"), 1.00, 2010.00)
        p3 = ShortPut(Ticker("SPX"), 2.00, 1990.00)

        c1 = LongCall(Ticker("SPX"), 1.50, 2000.00)
        c2 = ShortCall(Ticker("SPX"), 1.00, 2010.00)
        c3 = ShortCall(Ticker("SPX"), 2.00, 1990.00)

        spv = ShortVertical(p1, p2)
        spv = ShortVertical(p1, p2)
        lpv  = LongVertical(p1, p2)

        @pending typeof(p1 + p2) => Spread{ShortPut, LongPut}
        @pending typeof(p1 + p3) => Spread{LongPut, ShortPut}
        @pending typeof(p1 + p2) => LongPutVertical
        @pending typeof(p1 + p3) => ShortPutVertical
    end

    context("covered call") do

        stk  = Stock(Ticker("AAPL"), 100.00)
        stk  = LongStock(Ticker("AAPL"), 100.00)
        sstk = Stock(Ticker("AAPL"), -100.00)
        sstk = ShortStock(Ticker("AAPL"), 100.00)

        c1   = ShortCall(Ticker("AAPL"), 1.00, 120.00)
        p1   = ShortPut(Ticker("AAPL"), 1.00, 120.00)

        cc   = CoveredCall(stk, c1)

        @fact cc.risk => 99.00

        @pending typeof(stk + c1)  => Spread(Stock, ShortCall)
        @pending typeof(sstk + p1) => Spread(Stock, ShortPut)
        @pending typeof(stk + c1)  => CoveredCall
        @pending typeof(sstk + p1) => CoveredPut
    end

    context("butterfly") do

        p1 = LongPut(Ticker("SPX"), 2000.00, Date(2015,3,19))
        p2 = ShortPut(Ticker("SPX"), 2010.00, Date(2015,3,19))
        p3 = ShortPut(Ticker("SPX"), 1990.00, Date(2015,3,19))

        @pending typeof(p1 + 2p2 + p3) => Spread(ShortPutVertical, LongPutVertical)
        @pending typeof(p1 + 2p2 + p3) => LongPutButterfly
    end
end

ext: jl


text:
module FeatureDescriptors

export Descriptor
export categorical_keys, label, quantity_key, sources

"""
    Descriptor

The `Descriptor`s interface lets users define the properties and relationships between the
features being used in their models or other data-intensive workflows.

`Descriptor` supertypes define some general properties that can be overloaded for specific
subtypes as required, such as the tables to load or the quantitative variable to use.

Users may define their own Descriptor` type system for their own purposes but must extend
the following methods:
- [`sources`](@ref): The sources for the data needed to build the `Desriptor` feature.
- [`quantity_key`](@ref): The quantitative variable for the `Descriptor`.
- [`categorical_keys`](@ref): The categorical variables for the `Descriptor`.
"""
abstract type Descriptor end

"""
    sources(::Descriptor) -> Vector{String}

Returns the data sources for the [`Descriptor`](@ref), e.g. the names or paths to the tables
containing the required data: `["temperature.csv", "locations.csv"]`.
"""
function sources end

"""
    quantity_key(::Descriptor) -> Symbol

Returns the quantitative variable for the [`Descriptor`](@ref), e.g. `:temperature`.
"""
function quantity_key end

"""
    categorical_keys(::Descriptor) -> Vector

Returns the categorical variables for the [`Descriptor`](@ref), e.g. `[:colour, :shape]`.
Returns an empty vector if no categorical variables exist.
"""
function categorical_keys end

"""
    label(D::Descriptor) -> Symbol

Returns a symbol representation of the [`Descriptor`](@ref).
"""
function label(D::Type{<:Descriptor})
    s = replace(string(D), string(parentmodule(D)) => "")  # remove module name
    s = replace(lowercase(s), r"[^\w]" => "")  # remove any weird characters
    return Symbol(s)
end

# Allow sorting descriptors lexographically
# NOTE: isless(hash(A), hash(B)) would be faster, but:
# 1. we shouldn't need to sort that many descriptors
# 2. this is more human readable/understandable
Base.isless(A::Type{<:Descriptor}, B::Type{<:Descriptor}) = isless(string(A), string(B))

include("test_utils.jl")

end

max_stars_repo_name: invenia/FeatureDescriptors

max_issues_repo_name: invenia/FeatureDescriptors

max_forks_repo_name: invenia/FeatureDescriptors


text:
# ------------------------------------------------------------------------------
# operator fields
# ------------------------------------------------------------------------------

"""
```julia
OperatorField(
    basis,
    evaluationmodes
)
```

Finite Element operator input or output, with a basis and evaluation mode

# Arguments:
- `basis`:           finite element basis for the field
- `evaluationmodes`: array of basis evaluation modes,
                         note that quadrature weights must be listed in a
                         separate operator field

# Returns:
- Finite element operator field object

# Example:
```jldoctest
# basis
basis = TensorH1LagrangeBasis(4, 3, 2, 1);

# quadrature weights field, input only
weightsfield = OperatorField(basis, [EvaluationMode.quadratureweights], "quadrature weights");

# verify
println(weightsfield)

# input or output field
inputfield = OperatorField(basis, [
    EvaluationMode.interpolation,
    EvaluationMode.gradient,
    ],
    "gradient of weak form input"
);

# verify
println(inputfield)

# output
operator field:
  name:
    quadrature weights
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 3
    numbercomponents: 2
    dimension: 1
  evaluation mode:
    quadratureweights
operator field:
  name:
    gradient of weak form input
  tensor product basis:
    numbernodes1d: 4
    numberquadraturepoints1d: 3
    numbercomponents: 2
    dimension: 1
  evaluation modes:
    interpolation
    gradient
```
"""
struct OperatorField
    # data
    basis::AbstractBasis
    evaluationmodes::AbstractArray{EvaluationMode.EvalMode}
    name::String

    # inner constructor
    OperatorField(
        basis::AbstractBasis,
        evaluationmodes::AbstractArray{EvaluationMode.EvalMode},
    ) = (
        # validity checking
        if length(evaluationmodes) > 1 &&
           EvaluationMode.quadratureweights in evaluationmodes
            error("quadrature weights must be a separate operator field") # COV_EXCL_LINE
        end;

        # constructor
        new(basis, evaluationmodes)
    )

    # inner constructor
    OperatorField(
        basis::AbstractBasis,
        evaluationmodes::AbstractArray{EvaluationMode.EvalMode},
        name::String,
    ) = (
        # validity checking
        if length(evaluationmodes) > 1 &&
           EvaluationMode.quadratureweights in evaluationmodes
            error("quadrature weights must be a separate operator field") # COV_EXCL_LINE
        end;

        # constructor
        new(basis, evaluationmodes, name)
    )
end

# printing
# COV_EXCL_START
function Base.show(io::IO, field::OperatorField)
    print(io, "operator field:\n")

    # name
    if isdefined(field, :name)
        print(io, "  name:\n    ", field.name, "\n")
    end

    # basis
    print(io, "  ", field.basis)

    # evaluation modes
    if length(field.evaluationmodes) == 1
        print(io, "\n  evaluation mode:")
    else
        print(io, "\n  evaluation modes:")
    end
    for mode in field.evaluationmodes
        print(io, "\n    ", mode)
    end
end
# COV_EXCL_STOP

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


text:
import Zygote, ForwardDiff
using Random, Test
using YaoBlocks, YaoArrayRegister
using ChainRulesCore: Tangent

@testset "recursive_create_tangent" begin
    c = chain(put(5, 2 => chain(Rx(1.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    tc = YaoBlocks.AD.recursive_create_tangent(c)
    @test tc isa Tangent
end

@testset "construtors" begin
    @test Zygote.gradient(x -> x.list[1].blocks[1].theta, sum([chain(1, Rz(0.3))]))[1] == (n=nothing,
        list = NamedTuple{
            (:n, :blocks,),
            Tuple{Nothing, Vector{NamedTuple{(:block, :theta),Tuple{Nothing,Float64}}}},
        }[(n=nothing, blocks = [(block = nothing, theta = 1.0)],)],
    )
    @test_broken Zygote.gradient(
        x -> getfield(getfield(x, :content), :theta),
        Daggered(Rx(0.5)),
    )[1] == (content = (block = nothing, theta = 1.0),)
end

@testset "rules" begin
    h = put(5, 3 => Z) + put(5, 2 => X)
    c = chain(put(5, 2 => chain(Rx(1.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    r = rand_state(5)
    g0 = reinterpret(
        ComplexF64,
        ForwardDiff.gradient(
            x -> real(
                expect(h, ArrayReg([Complex(x[2i-1], x[2i]) for i = 1:length(x)÷2])),
            ),
            reinterpret(Float64, r.state),
        ),
    )
    @test Zygote.gradient(x -> real(expect(h, ArrayReg(x))), r.state)[1] ≈ g0
    @test Zygote.gradient(
        x -> real(expect(h, ArrayReg(reshape(statevec(x), :, 1)))),
        r,
    )[1].state ≈ g0
    @test Zygote.gradient(
        x -> real(expect(h, ArrayReg(reshape(state(x), :, 1)))),
        r,
    )[1].state ≈ g0
    @test Zygote.gradient(x -> real(expect(h, copy(x))), r)[1].state ≈ g0
    @test Zygote.gradient(x -> real(expect(h, parent(x'))), r)[1].state ≈ g0

    g1 = reinterpret(
        ComplexF64,
        ForwardDiff.gradient(
            x -> real(sum(abs2, [Complex(x[2i-1], x[2i]) for i = 1:length(x)÷2])),
            reinterpret(Float64, r.state),
        ),
    )
    @test Zygote.gradient(x -> real(sum(abs2, state(x'))), r)[1].state ≈ g1
    @test Zygote.gradient(x -> real(sum(abs2, statevec(x'))), r)[1].state ≈ g1
    # zygote does not work if `sin` is not here,
    # because it gives an adjoint of different type as the output matrix type.
    @test AD.extract_circuit_gradients!(
        Zygote.gradient(x -> real(sum(sin, Matrix(x))), c)[1].blocks,
        Float64[],
    ) ≈ ForwardDiff.gradient(x -> real(sum(sin, Matrix(dispatch(c, x)))), parameters(c))
end

@testset "adwith zygote" begin
    c = chain(put(5, 2 => chain(Rx(0.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    dispatch!(c, :random)

    function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        reg = apply(copy(reg), circuit)
        st = state(reg)
        sum(real(st .* st))
    end

    # apply
    reg0 = zero_state(5)
    params = rand!(parameters(c))
    paramsδ = Zygote.gradient(params -> loss(reg0, dispatch(c, params)), params)[1]
    regδ = Zygote.gradient(reg -> loss(reg, c), reg0)[1]
    fparamsδ = ForwardDiff.gradient(
        params -> loss(
            ArrayReg(Matrix{Complex{eltype(params)}}(reg0.state)),
            dispatch(c, params),
        ),
        params,
    )
    fregδ = ForwardDiff.gradient(
        x -> loss(
            ArrayReg([Complex(x[2i-1], x[2i]) for i = 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(parameters(c))),
        ),
        reinterpret(Float64, reg0.state),
    )
    @test fregδ ≈ reinterpret(Float64, regδ.state)
    @test fparamsδ ≈ paramsδ

    # expect and fidelity
    c = chain(
        put(5, 5 => Rx(1.5)),
        put(5, 1 => Rx(0.4)),
        put(5, 4 => Rx(0.2)),
        put(5, 2 => chain(Rx(0.4), Rx(0.5))),
        cnot(5, 3, 1),
        put(5, 3 => Rx(-0.5)),
    )
    h = chain(repeat(5, X, 1:5))
    reg = rand_state(5)
    function loss2(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        return 5 *
               real(expect(h, copy(reg) => circuit) + fidelity(reg, apply(reg, circuit)))
    end
    params = rand!(parameters(c))
    fδc = ForwardDiff.gradient(
        params -> loss2(
            ArrayReg(Matrix{Complex{eltype(params)}}(reg.state)),
            dispatch(c, params),
        ),
        params,
    )
    δr, δc = Zygote.gradient((reg, params) -> loss2(reg, dispatch(c, params)), reg, params)
    @test δc ≈ fδc

    fregδ = ForwardDiff.gradient(
        x -> loss2(
            ArrayReg([Complex(x[2i-1], x[2i]) for i = 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(params)),
        ),
        reinterpret(Float64, reg.state),
    )
    @test fregδ ≈ reinterpret(Float64, δr.state)

    # operator fidelity
    c = chain(
        put(5, 5 => Rx(1.5)),
        put(5, 1 => Rx(0.4)),
        put(5, 4 => Rx(0.2)),
        put(5, 2 => chain(Rx(0.4), Rx(0.5))),
        cnot(5, 3, 1),
        put(5, 3 => Rx(-0.5)),
    )
    h = chain(repeat(5, X, 1:5))
    function loss3(circuit::AbstractBlock{N}, h) where {N}
        return operator_fidelity(circuit, h)
    end
    params = rand!(parameters(c))
    fδc = ForwardDiff.gradient(params -> loss3(dispatch(c, params), h), params)
    δc, = Zygote.gradient(p -> loss3(dispatch(c, p), h), params)
    @test δc ≈ fδc

    # NOTE: operator back propagation in expect is not implemented!
    # to differentiate operators, we need to use the expensive `mat_back` function.
end

@testset "add block" begin
    H = sum([chain(5, put(k => Z)) for k = 1:5])
    c = chain(put(5, 2 => chain(Rx(0.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    dispatch!(c, :random)
    function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        reg = apply(copy(reg), circuit)
        st = state(reg)
        reg2 = apply(copy(reg), H)
        st2 = state(reg2)
        sum(real(st .* st2))
    end

    reg0 = zero_state(5)
    params = rand!(parameters(c))
    paramsδ = Zygote.gradient(params -> loss(reg0, dispatch(c, params)), params)[1]
    fparamsδ = ForwardDiff.gradient(
        params -> loss(
            ArrayReg(Matrix{Complex{eltype(params)}}(reg0.state)),
            dispatch(c, params),
        ),
        params,
    )
    @test fparamsδ ≈ paramsδ

    regδ = Zygote.gradient(reg -> loss(reg, c), reg0)[1]
    fregδ = ForwardDiff.gradient(
        x -> loss(
            ArrayReg([Complex(x[2i-1], x[2i]) for i = 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(parameters(c))),
        ),
        reinterpret(Float64, reg0.state),
    )
    @test fregδ ≈ reinterpret(Float64, regδ.state)
end

@testset "apply hami1" begin
    h = 0.2 * put(5, 3 => Z) + 0.5 * put(5, 2 => X) + 0.4 * repeat(5, Z, 1:3)
    reg0 = rand_state(5)
    loss(reg, params) = fidelity(reg0, apply(reg, dispatch(h, params)))

    reg = rand_state(5)
    params = randn(3)
    _adapt(::Type, reg) = reg
    _adapt(::Type{T}, reg::AbstractArrayReg) where T <:ForwardDiff.Dual = similar(reg, _adapt.(T, reg.state))
    _adapt(::Type{T}, x::Complex) where T <: ForwardDiff.Dual = T(x.re) + im * T(x.im)
    t1 = ForwardDiff.gradient(state->loss(ArrayReg{2}(state[1:32] .+ im .* state[33:64]), params), [real.(reg.state)..., imag.(reg.state)...])
    t2 = ForwardDiff.gradient(params->loss(_adapt(eltype(params), reg), params), params)
    g1, g2 = Zygote.gradient(loss, reg, params)
    @test g2 ≈ t2
    @test [real.(g1.state)..., imag.(g1.state)...] ≈ t1
end

@testset "apply hami2" begin
    N = 6
    cost = sum([put(N, i=>Z) for i in 1:N])
    theta = rand(96)
    c = chain(rand([[control(N, i, mod1(i+1,N)=>Rx(0.0)) for i=1:N]..., [put(N, i=>Ry(0.0)) for i=1:N]..., [put(N, i=>Rx(0.0)) for i=1:N]...], 96))
    new_var = dispatch(c, :random);
    function new_loss(theta::AbstractVector{T}) where T
        new_circ = dispatch(new_var, theta)
        diff_evalue = expect(cost, zero_state(Complex{T}, N) => new_circ)
        return 3*real.(diff_evalue)
    end

    g1 = Zygote.gradient(_theta->new_loss(_theta), theta)[1]
    g2 = ForwardDiff.gradient(_theta->new_loss(_theta), theta)
    @test g1 ≈ g2

    function new_loss2(reg::ArrayReg{T}) where T
        diff_evalue = expect(cost, reg => new_var)
        return 3*real.(diff_evalue)
    end
    regin = rand_state(N)
    g1 = Zygote.gradient(_theta->new_loss2(_theta), regin)[1]
    params = reinterpret(real(eltype(regin.state)), regin.state)
    g2 = ForwardDiff.gradient(p->new_loss2(arrayreg(reinterpret(complex(eltype(p)), p))), params)
    @test statevec(g1) ≈ reinterpret(ComplexF64, g2)
end

@testset "fix #348" begin
    n = 3
    h = sum([kron(n, i => Z, i+1=>Z) for i = 1:n-1])
    function wrap(U::AbstractBlock, θ::Vector{T}) where {T}
        n = nqubits(U)
        reg = apply(zero_state(Complex{T}, n), dispatch(U, θ))
        real(expect(h, reg))
    end

    hb = sum([put(n, i => X) for i = 1:n])
    U2 = chain(n, [time_evolve(hb, 0.0)])

    g1 = (wrap(U2, [0.618 + 1e-5]) - wrap(U2, [0.618-1e-5]))/2e-5
    g2 = Zygote.gradient(p -> wrap(U2, p), [0.618])[1][]
    @test isapprox(g1, g2, rtol=1e-3)
end


text:
# Advent of Code 2020 - Day 21

test_str = "mxmxvkd kfcds sqjhc nhms (contains dairy, fish)
trh fvjkl sbzzf mxmxvkd (contains dairy)
sqjhc fvjkl (contains soy)
sqjhc mxmxvkd sbzzf (contains fish)"

test_list = split(test_str, '\n')

list = readlines("day21_input.txt")

function count_non_allergen_appearance(list)
    allergen_assesment = Dict{String,Set{String}}()
    ingredient_count = Dict{String,Int}()
    for i in 1:length(list)
        m = match(r"(^.*)(?: \(contains )(.*)\)$", list[i])
        ingredients = Set(String.(split(m[1], " ")))
        allergens = Set(String.(split(m[2], ", ")))
        for allergen in allergens
            if allergen in keys(allergen_assesment)
                intersect!(allergen_assesment[allergen], ingredients)
            else
                push!(allergen_assesment, allergen => copy(ingredients)) # if we don't use copy we end up referencing the same object for different keys and changing one affects the other.
            end
        end
        for ingredient in ingredients
            if ingredient in keys(ingredient_count)
                ingredient_count[ingredient] += 1
            else
                push!(ingredient_count, ingredient => 1)
            end
        end
    end

    while maximum(length.(values(allergen_assesment))) > 1
        for (k,v) in allergen_assesment
            if length(v) == 1
                for j in keys(allergen_assesment)
                    if j != k
                        setdiff!(allergen_assesment[j], v)
                    end
                end
            end
        end
    end

    total = sum([v for (k,v) in ingredient_count if k ∉ union(values(allergen_assesment)...)])

    dangerous_ingredients = ""
    for allergen in sort(collect(keys(allergen_assesment)))
        dangerous_ingredients *= unique(allergen_assesment[allergen])[1]*","
    end
    dangerous_ingredients_list = dangerous_ingredients[1:end-1]

    return total, dangerous_ingredients_list
end

@show count_non_allergen_appearance(test_list) == (5, "mxmxvkd,sqjhc,fvjkl")
@show count_non_allergen_appearance(list) == (2584, "fqhpsl,zxncg,clzpsl,zbbnj,jkgbvlxh,dzqc,ppj,glzb")


text:
using GeoMakie, AbstractPlotting.MakieLayout, Makie

lons = LinRange(-179.5, 179.5, 360)
lats = LinRange(-89.5, 89.5, 180)

field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

source = LonLat()
dest   = WinkelTripel()
cs     = Proj4.CRS2CRS(source, dest)

xs, ys = xygrid(lons, lats)

xmin, xmax = extrema(xs)
ymin, ymax = extrema(ys)

aspect_ratio = (ymax - ymin) / (xmax - xmin)

scene, layout = layoutscene(
);

layout[1, 1] = lsc = LScene(scene; scenekw = (show_axis = false, scale_plot = false, transform_func = cs))

splot = surface!(lsc, xs, ys; color = field, shading = false, show_axis = false)

geoaxis!(lsc, Rect2D(-180, -90, 360, 180))

coastlines!(lsc)

layout[1, 2] = LColorbar(scene, splot; label = "Arbitrary data", width = 30)

colsize!(layout, 1, Relative(1))
rowsize!(layout, 1, Aspect(1, aspect_ratio))

layout[0, :] = LText(scene, "MakieLayout is cool!"; textsize = 40)

scene |> save("geo.png")
scene |> save("geo.pdf")


text:
using OMEinsum
using OMEinsumContractionOrders
using Random

function generate_vertical_rules(;Nv=2,χ=800)
	eincode = EinCode(((1,2,3,4,5),# T1
	(6,7,8,9,10),#T2 (dag)
	(2,6,11,12), #swapgate(nl,nu)
	(8,10,14,13), #swapgate(nf,nu)
	(3,13,15,16), #swapgate(nf,nr)
	(4,16,17,18), #swapgate(nl,nu)
	(17,19,20,21,22), #T4
	(23,24,25,26,27), #T3 (dag)
	(21,27,28,29),#swapgate(nl,nu)
	(20,19,30,31),#swapgate(nf,nr)
	(25,31,32,33),#swapgate(nf,nr)
	(23,33,9,34), #swapgate(nl,nu)
	(35,12,1,36), # ACu: E3
	(36,18,5,37), # FRu: E8
	(37,29,22,38), # FRo: E4
	(39,26,28,38),# ACd: E6
	(40,24,34,39), # FLo: E1
	(35,7,11,40) # FLu: E7
	),
	(15,30,14,32) #hamiltonian (ij di dj)
	)

	size_dict = [2^Nv for i = 1:40]
	size_dict[[3;8;14;15;30;32;20;25]] .= 4
	size_dict[35:40] .= χ
	sd = Dict(i=>size_dict[i] for i = 1:40)

	# for seed =40:100
	seed = 60
	Random.seed!(60)
	optcode = optimize_tree(eincode,sd; sc_target=28, βs=0.1:0.1:10, ntrials=2, niters=100, sc_weight=3.0)

	print("Vertical Contraction Complexity(seed=$(seed))",OMEinsum.timespace_complexity(optcode,sd),"\n")
	# You would better try some times to make it optimal (By the for-end iteration...)
	# end
	return optcode
end
const VERTICAL_RULES = generate_vertical_rules()

function generate_horizontal_rules(;Nv=2,χ=800)
    eincode = EinCode(((1,2,3,4,5),# T1
    (3,4,21,22),#swapgate(nf,nu)

	(6,7,8,9,10),#T2 (dag)
    (8,22,23,27),#swapgate(nf,nu)
    (10,40,17,27), #swapgate(nl,nu)

    (2,6,29,30), #swapgate(nl,nu)

    (16,17,18,19,20),# T4(dag)
    (18,26,25,16),#swapgate(nf,nu)

    (13,28,24,26),#swapgate(nf,nu)
    (12,28,5,33),#swapgate(nl,nu)
    (11,12,13,14,15),# T3

    (20,32,31,14),#swapgate(nl,nu)

    (39,7,29,38), #E1 FLo
    (39,30,1,34), #E2 ACu
    (34,33,11,35), #E3 ARu
    (35,31,15,36), #E4 FRo
    (37,19,32,36), #E5 ARd
    (38,9,40,37), #E6 ACd
	),
	(21,24,23,25) #hamiltonian (ij di dj)
	)

    size_dict = [2^Nv for i = 1:40]
	size_dict[[3;8;13;18;21;23;24;25]] .= 4
	size_dict[34:39] .= χ
	sd = Dict(i=>size_dict[i] for i = 1:40)

	# for seed = 1:100
	seed = 4
	Random.seed!(seed)
	optcode = optimize_tree(eincode,sd; sc_target=28, βs=0.1:0.1:10, ntrials=3, niters=100, sc_weight=4.0)
	print("Horizontal Contraction Complexity(seed=$(seed))",OMEinsum.timespace_complexity(optcode,sd),"\n")
	# end

	return optcode
end
const HORIZONTAL_RULES = generate_horizontal_rules()


text:
using GeometricBase
using GeometricBase: fromarray, test_interface
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solutions: createHDF5
using GeometricProblems.HarmonicOscillator
using Test
import HDF5

nt    = 10
Δt    = .1
ni    = 5

t0    = 0.
x0    = rand(2)
y0    = rand(2)
q0    = rand(1)
p0    = rand(1)
z0    = rand(3)
λ0    = rand(1)
μ0    = rand(2)

t1    = 1.
x1    = [rand(2) for j in 1:ni]
y1    = [rand(2) for j in 1:ni]
q1    = [rand(1) for j in 1:ni]
p1    = [rand(1) for j in 1:ni]
z1    = [rand(3) for j in 1:ni]
λ1    = [rand(1) for j in 1:ni]
μ1    = [rand(2) for j in 1:ni]

t2    = t1 + (t1-t0)
x2    = [rand(2) for j in 1:ni]
y2    = [rand(2) for j in 1:ni]
q2    = [rand(1) for j in 1:ni]
p2    = [rand(1) for j in 1:ni]
z2    = [rand(3) for j in 1:ni]
λ2    = [rand(1) for j in 1:ni]
μ2    = [rand(2) for j in 1:ni]

tx    = zero(x0)
ty    = zero(y0)
tq    = zero(q0)
tp    = zero(p0)
tz    = zero(z0)
tλ    = zero(λ0)
tμ    = zero(μ0)

xs    = [rand(2) for i in 1:nt]
ys    = [rand(2) for i in 1:nt]
qs    = [rand(1) for i in 1:nt]
ps    = [rand(1) for i in 1:nt]
zs    = [rand(3) for i in 1:nt]
λs    = [rand(1) for i in 1:nt]
μs    = [rand(2) for i in 1:nt]

Xs    = [rand(2) for i in 1:nt, j in 1:ni]
Ys    = [rand(2) for i in 1:nt, j in 1:ni]
Qs    = [rand(1) for i in 1:nt, j in 1:ni]
Ps    = [rand(1) for i in 1:nt, j in 1:ni]
Zs    = [rand(3) for i in 1:nt, j in 1:ni]
Λs    = [rand(1) for i in 1:nt, j in 1:ni]
Ms    = [rand(2) for i in 1:nt, j in 1:ni]

ode   = harmonic_oscillator_ode()
dae   = harmonic_oscillator_dae()
pode  = harmonic_oscillator_pode()
pdae  = harmonic_oscillator_pdae()

x100  = rand(ndims(ode),  100)
y100  = rand(ndims(ode),  100)
z100  = rand(ndims(dae),  100)
q100  = rand(ndims(pode), 100)
p100  = rand(ndims(pode), 100)
λ100  = rand(dae.m,  100)
μ100  = rand(pdae.m, 100)

h5file = "test.hdf5"

@testset "$(rpad("ODE Solution",80))" begin
    asol = AtomicSolution(ode)

    # test constructors and general functionality
    sol = Solution(ode, Δt, nt)
    @test typeof(sol) <: SolutionODE

    test_interface(sol)

    sol0 = Solution(similar(ode, x0), Δt, nt)
    @test typeof(sol0) <: SolutionODE

    sol1 = Solution(similar(ode, x1), Δt, nt)
    @test typeof(sol1) <: SolutionODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0)
    get_initial_conditions!(sol, tx, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test asol.t == t0
    @test asol.q == x0

    δt, δx = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0

    set_initial_conditions!(sol1, similar(ode, t1, x2))
    for j in 1:ni
        get_initial_conditions!(sol1, tx, j)
        @test tx == x2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(ode, x0), Δt, nt)
    sol2 = Solution(similar(ode, x0), Δt, nt)
    for i in 1:nt
        tx .= xs[i]
        asol.q .= xs[i]
        set_solution!(sol1, tx, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == xs
    @test sol2.q[1:nt] == xs

    sol1 = Solution(similar(ode, x1), Δt, nt)
    sol2 = Solution(similar(ode, x1), Δt, nt)
    for i in 1:nt
        for j in 1:ni
            tx .= Xs[i,j]
            asol.q .= Xs[i,j]
            set_solution!(sol1, tx, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt,:] == Xs
    @test sol2.q[1:nt,:] == Xs

    # test general hdf5 functions
    sol1 = Solution(ode, Δt, nt)
    h5 = createHDF5(sol1, h5file)
    @test typeof(h5) == HDF5.File
    close(h5)
    @test isfile(h5file)

    sol1 = Solution(ode, Δt, nt)
    h5 = createHDF5(sol1, h5file; overwrite=false)
    @test typeof(h5) == HDF5.File
    close(h5)
    @test isfile(h5file)
    rm(h5file)

    sol1 = Solution(ode, Δt, nt)
    write_to_hdf5(sol1, h5file)
    @test isfile(h5file)
    rm(h5file)

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_ode(x0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_ode(x1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(ode, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(ode, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(ode, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(ode, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == vcat(ode.q₀, [x100[:,j] for j in axes(x100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(ode, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == vcat(ode.q₀, [x100[:,j] for j in axes(x100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("PODE Solution",80))" begin
    asol = AtomicSolution(pode)

    # test constructors and general functionality
    sol = Solution(pode, Δt, nt)
    @test typeof(sol) <: SolutionPODE

    test_interface(sol)

    sol0 = Solution(similar(pode, q0, p0), Δt, nt)
    @test typeof(sol0) <: SolutionPODE

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    @test typeof(sol1) <: SolutionPODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, q0, p0)
    get_initial_conditions!(sol, tq, tp, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tq == q0
    @test tp == p0
    @test asol.t == t0
    @test asol.q == q0
    @test asol.p == p0

    δt, δq, δp = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δq == q0
    @test δp == p0

    set_initial_conditions!(sol1, similar(pode, t1, q2, p2))
    for j in 1:ni
        get_initial_conditions!(sol1, tq, tp, j)
        @test tq == q2[j]
        @test tp == p2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(pode, q0, p0), Δt, nt)
    sol2 = Solution(similar(pode, q0, p0), Δt, nt)
    for i in 1:nt
        tq .= qs[i]
        tp .= ps[i]
        asol.q .= qs[i]
        asol.p .= ps[i]
        set_solution!(sol1, tq, tp, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == qs
    @test sol1.p[1:nt] == ps
    @test sol2.q[1:nt] == qs
    @test sol2.p[1:nt] == ps

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    sol2 = Solution(similar(pode, q1, p1), Δt, nt)
    for i in 1:nt
        for j in 1:ni
            tq .= Qs[i,j]
            tp .= Ps[i,j]
            asol.q .= Qs[i,j]
            asol.p .= Ps[i,j]
            set_solution!(sol1, tq, tp, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt,:] == Qs
    @test sol1.p[1:nt,:] == Ps
    @test sol2.q[1:nt,:] == Qs
    @test sol2.p[1:nt,:] == Ps

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_pode(q0, p0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pode(q1, p1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(pode, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(pode, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pode, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(pode, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == vcat(pode.q₀, [q100[:,j] for j in axes(q100, 2)])
    @test sol2.p.d == vcat(pode.p₀, [p100[:,j] for j in axes(p100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pode, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == vcat(pode.q₀, [q100[:,j] for j in axes(q100, 2)[2:2:end]])
    @test sol2.p.d == vcat(pode.p₀, [p100[:,j] for j in axes(p100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("DAE Solution",80))" begin
    asol = AtomicSolution(dae)

    # test constructors and general functionality
    sol = Solution(dae, Δt, nt)
    @test typeof(sol) <: SolutionDAE

    test_interface(sol)

    sol0 = Solution(similar(dae, z0, λ0), Δt, nt)
    @test typeof(sol0) <: SolutionDAE

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    @test typeof(sol1) <: SolutionDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, z0, λ0)
    get_initial_conditions!(sol, tz, tλ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tz == z0
    @test tλ == λ0
    @test asol.t == t0
    @test asol.q == z0
    @test asol.λ == λ0

    δt, δz, δλ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δz == z0
    @test δλ == λ0

    set_initial_conditions!(sol1, similar(dae, t1, z2, λ2))
    for j in 1:ni
        get_initial_conditions!(sol1, tz, tλ, j)
        @test tz == z2[j]
        @test tλ == λ2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(dae, z0, λ0), Δt, nt)
    sol2 = Solution(similar(dae, z0, λ0), Δt, nt)
    for i in 1:nt
        tz .= zs[i]
        tλ .= λs[i]
        asol.q .= zs[i]
        asol.λ .= λs[i]
        set_solution!(sol1, tz, tλ, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[1:nt] == zs
    @test sol1.λ[1:nt] == λs
    @test sol2.q[1:nt] == zs
    @test sol2.λ[1:nt] == λs

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    sol2 = Solution(similar(dae, z1, λ1), Δt, nt)
    for i in 1:nt
        for j in 1:ni
            tz .= Zs[i,j]
            tλ .= Λs[i,j]
            asol.q .= Zs[i,j]
            asol.λ .= Λs[i,j]
            set_solution!(sol1, tz, tλ, i, j)
            set_solution!(sol2, asol, i, j)
        end
    end
    @test sol1.q[1:nt,:] == Zs
    @test sol1.λ[1:nt,:] == Λs
    @test sol2.q[1:nt,:] == Zs
    @test sol2.λ[1:nt,:] == Λs

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_dae(z0, λ0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_dae(z1, λ1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(dae, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(dae, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(dae, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(dae, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, z100[:,oj], λ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == vcat(dae.q₀, [z100[:,j] for j in axes(z100, 2)])
    @test sol2.λ.d == vcat(dae.λ₀, [λ100[:,j] for j in axes(λ100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(dae, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, z100[:,oj], λ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == vcat(dae.q₀, [z100[:,j] for j in axes(z100, 2)[2:2:end]])
    @test sol2.λ.d == vcat(dae.λ₀, [λ100[:,j] for j in axes(λ100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("PDAE Solution",80))" begin
    asol = AtomicSolution(pdae)

    # test constructors and general functionality
    sol = Solution(pdae, Δt, nt)
    @test typeof(sol) <: SolutionPDAE

    test_interface(sol)

    sol0 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    @test typeof(sol0) <: SolutionPDAE

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    @test typeof(sol1) <: SolutionPDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0, y0, μ0)
    get_initial_conditions!(sol, tx, ty, tμ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test ty == y0
    @test tμ == μ0
    @test asol.t == t0
    @test asol.q == x0
    @test asol.p == y0
    @test asol.λ == μ0

    tx .= 0
    ty .= 0
    get_initial_conditions!(sol, tx, ty, 1)
    @test tx == x0
    @test ty == y0

    δt, δx, δy, δμ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0
    @test δy == y0
    @test δμ == μ0

    set_initial_conditions!(sol1, similar(pdae, t1, x2, y2, μ2))
    for j in 1:ni
        get_initial_conditions!(sol1, tx, ty, tμ, j)
        @test tx == x2[j]
        @test ty == y2[j]
        @test tμ == μ2[j]
    end

    # test set/get solution
    sol1 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol2 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol3 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    for i in 1:nt
        tx .= xs[i]
        ty .= ys[i]
        tμ .= μs[i]
        asol.q .= xs[i]
        asol.p .= ys[i]
        asol.λ .= μs[i]
        set_solution!(sol1, tx, ty, tμ, i)
        set_solution!(sol2, tx, ty, i)
        set_solution!(sol3, asol, i)
    end
    @test sol1.q[1:nt] == xs
    @test sol1.p[1:nt] == ys
    @test sol1.λ[1:nt] == μs
    @test sol2.q[1:nt] == xs
    @test sol2.p[1:nt] == ys
    @test sol3.q[1:nt] == xs
    @test sol3.p[1:nt] == ys
    @test sol3.λ[1:nt] == μs

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol2 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol3 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    for i in 1:nt
        for j in 1:ni
            tx .= Xs[i,j]
            ty .= Ys[i,j]
            tμ .= Ms[i,j]
            asol.q .= Xs[i,j]
            asol.p .= Ys[i,j]
            asol.λ .= Ms[i,j]
            set_solution!(sol1, tx, ty, tμ, i, j)
            set_solution!(sol2, tx, ty, i, j)
            set_solution!(sol3, asol, i, j)
        end
    end
    @test sol1.q[1:nt,:] == Xs
    @test sol1.p[1:nt,:] == Ys
    @test sol1.λ[1:nt,:] == Ms
    @test sol2.q[1:nt,:] == Xs
    @test sol2.p[1:nt,:] == Ys
    @test sol3.q[1:nt,:] == Xs
    @test sol3.p[1:nt,:] == Ys
    @test sol3.λ[1:nt,:] == Ms

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_pdae(x0, y0, μ0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pdae(x1, y1, μ1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(pdae, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(pdae, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pdae, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(pdae, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, x100[:,oj], y100[:,oj], μ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == vcat(pdae.q₀, [x100[:,j] for j in axes(x100, 2)])
    @test sol2.p.d == vcat(pdae.p₀, [y100[:,j] for j in axes(y100, 2)])
    @test sol2.λ.d == vcat(pdae.λ₀, [μ100[:,j] for j in axes(μ100, 2)])
    # @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)
    # @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)
    # @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pdae, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, x100[:,oj], y100[:,oj], μ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == vcat(pdae.q₀, [x100[:,j] for j in axes(x100, 2)[2:2:end]])
    @test sol2.p.d == vcat(pdae.p₀, [y100[:,j] for j in axes(y100, 2)[2:2:end]])
    @test sol2.λ.d == vcat(pdae.λ₀, [μ100[:,j] for j in axes(μ100, 2)[2:2:end]])
    # @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)[:,1:2:end]
    # @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)[:,1:2:end]
    # @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end


text:
const irregular = Dict("one" => "first", "two" => "second", "three" => "third", "five" => "fifth",
                                "eight" => "eighth", "nine" => "ninth", "twelve" => "twelfth")
const suffix = "th"
const ysuffix = "ieth"


text:
# This script uses the functions in sandpile.jl
include("Sandpile.jl")
using .Sandpile
using Plots

x, y = 50, 50; # Dimensions of grid
N = 100_000; # Number of grains to be added
fᵪ = 4; # Critical value for sandpile model

for samples in [1, 10, 200]
    for init in ["zero", "random"]

            s = [run_sandpile!(sandpile_init(x, y, "zero"), N) for i in 1:samples]

        # Get counts of each slide size s
        unique_s = sort(unique(reduce(vcat, s)))
        counts = [count(==(element), reduce(vcat, s)) for element in unique_s]
        yticks = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
        xticks = [1e0, 1e1, 1e2, 1e3]
        plot(
            unique_s,
            counts ./ sum(counts);
            xaxis=:log,
            yaxis=:log,
            leg=false,
            yticks=yticks,
            xticks=xticks,
            size=(600, 600),
            dpi=300,
        )
        plot!(
            [1e0, 1e5],
            [1.25e-1, 1.25e-6];
            xaxis=:log,
            yaxis=:log,
            xlabel="s",
            ylabel="D(s)",
            xlim=(1e0, 1.5e3),
            ylim=(1e-5, 1.5e-1),
            linestyle=:dash,
            dpi=300,
        )
        savefig("../results/sandpile_dim$(x)x$(y)_N$(N)_init$(init)_$(samples)samples.png")

    end
end

max_stars_repo_name: romeroraa/sandpile-julia

max_issues_repo_name: romeroraa/sandpile-julia

max_forks_repo_name: romeroraa/sandpile-julia


text:
data12 = """
hl-WP
vl-fo
vl-WW
WP-start
vl-QW
fo-wy
WW-dz
dz-hl
fo-end
VH-fo
ps-vl
FN-dz
WP-ps
ps-start
WW-hl
end-QW
start-vl
WP-fo
end-FN
hl-QW
WP-dz
QW-fo
QW-dz
ps-dz
"""

""" recursive depth first search """
function allpaths(dverts, startvertex, currentpath, part1=true)
    startvertex == "end" && return [currentpath]
    paths = Vector{String}[]
    for nextvertex in dverts[startvertex]
        nextvertex == "start" && continue
        if isuppercase(first(nextvertex)) || !(nextvertex in currentpath)
            append!(
                paths, allpaths(dverts, nextvertex, vcat(currentpath, [nextvertex]), part1)
            )
        elseif !part1 && islowercase(first(nextvertex))
            append!(paths, allpaths(dverts, nextvertex, vcat(currentpath, [nextvertex])))
        end
    end
    return paths
end

function day12()
    part = [0, 0]
    day12lines = filter(!isempty, strip.(split(data12, "\n")))
    vertices = Dict{String,Vector{String}}()
    for line in day12lines
        v1, v2 = split(line, "-")
        push!(get!(vertices, v1, String[]), v2)
        push!(get!(vertices, v2, String[]), v1)
    end

    part[1] = length(allpaths(vertices, "start", ["start"]))
    part[2] = length(allpaths(vertices, "start", ["start"], false))
    return part
end

part = day12()
println("Part 1:", part[1])
println("Part 2:", part[2])


text:
# ---
# title: 1219. Path with Maximum Gold
# id: problem1219
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/path-with-maximum-gold/description/>
# hidden: true
# ---
#
# In a gold mine `grid` of size `m * n`, each cell in this mine has an integer
# representing the amount of gold in that cell, `0` if it is empty.
#
# Return the maximum amount of gold you can collect under the conditions:
#
#   * Every time you are located in a cell you will collect all the gold in that cell.
#   * From your position you can walk one step to the left, right, up or down.
#   * You can't visit the same cell more than once.
#   * Never visit a cell with `0` gold.
#   * You can start and stop collecting gold from  **any** position in the grid that has some gold.
#
#
#
# **Example 1:**
#
#
#
#     Input: grid = [[0,6,0],[5,8,7],[0,9,0]]
#     Output: 24
#     Explanation:
#     [[0,6,0],
#      [5,8,7],
#      [0,9,0]]
#     Path to get the maximum gold, 9 -> 8 -> 7.
#
#
# **Example 2:**
#
#
#
#     Input: grid = [[1,0,7],[2,0,6],[3,4,5],[0,3,0],[9,0,20]]
#     Output: 28
#     Explanation:
#     [[1,0,7],
#      [2,0,6],
#      [3,4,5],
#      [0,3,0],
#      [9,0,20]]
#     Path to get the maximum gold, 1 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7.
#
#
#
#
# **Constraints:**
#
#   * `1 <= grid.length, grid[i].length <= 15`
#   * `0 <= grid[i][j] <= 100`
#   * There are at most **25  **cells containing gold.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using HDF5
using LsqFit
using Glob

function derive_footprint_fits(fname::String)

    nc = h5open(fname, "r")

    # Find out what footprints we have
    fp = nc["Metadata/FootprintId"][:]
    un_fp = unique(nc["Metadata/FootprintId"][:])
    sort!(un_fp)

    SIF_ucert = nc["SIF_Uncertainty_740nm"][:]
    cont = nc["Science/continuum_radiance_757nm"][:]
    qual = nc["Quality_Flag"][:]

    output = Dict()

    @. model(x, p) = p[1] + p[2] * x

    for this_fp in un_fp
        idx = findall((qual .== 0) .& (fp .== this_fp))
        fit = curve_fit(model, cont[idx], (SIF_ucert[idx]).^2,
                        [0.0, 1.0],
                        lower=[0.0, 0.0], upper=[Inf, Inf])

        output[this_fp] = coef(fit)
    end

    return output

end


text:
using Parquet
using Test

function test_decode(file, parcompat=joinpath(dirname(@__FILE__), "parquet-compatibility"))
    p = ParFile(joinpath(parcompat, file))
    println("loaded $file")
    @test isa(p.meta, Parquet.FileMetaData)

    rgs = rowgroups(p)
    @test length(rgs) > 0
    println("\tfound $(length(rgs)) row groups")

    for rg in rgs
        ccs = columns(p, rg)
        println("\treading row group with $(length(ccs)) column chunks")

        for cc in ccs
            pgs = pages(p, cc)
            println("\t\treading column chunk with $(length(pgs)) pages, $(colname(cc))")

            for pg in pgs
                bpg = bytes(pg)
                println("\t\t\tread page with $(length(bpg)) bytes")
                if (pg.hdr._type == Parquet.PageType.DATA_PAGE) || (pg.hdr._type == Parquet.PageType.DATA_PAGE_V2)
                    #println("\t\t\treading data page")
                    vals, defn_levels, repn_levels = values(p, pg)
                    println("\t\t\tread data page, $(length(defn_levels)) defn_levels, $(length(repn_levels)) repn_levels, $(length(vals)) values")
                elseif pg.hdr._type == Parquet.PageType.DICTIONARY_PAGE
                    #println("\t\t\treading dictionary page")
                    vals = values(p, pg)[1]
                    println("\t\t\tread dictionary page, $(length(vals)) values")
                elseif pg.hdr._type == Parquet.PageType.INDEX_PAGE
                    println("\t\t\tnot reading index page yet")
                else
                    error("unknown page type $(pg.hdr._type)")
                end
            end

            println("\t\tre-reading column chunk for values. total $(length(pgs)) pages")
            vals, defn_levels, repn_levels = values(p, cc)
            println("\t\t\tread $(length(vals)) values, $(length(defn_levels)) defn levels, $(length(repn_levels)) repn levels")
        end
    end

    println("\tsuccess")
end

function test_decode_all_pages()
    for encformat in ("SNAPPY", "GZIP", "NONE")
        for fname in ("nation", "customer")
            test_decode("parquet-testdata/impala/1.1.1-$encformat/$fname.impala.parquet")
        end
    end
end

test_decode_all_pages()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "PulseAudio"
version = v"13.0.0-dev"

# Collection of sources required to build PulseAudio
sources = [
    "http://anongit.freedesktop.org/git/pulseaudio/pulseaudio.git" =>
    "a5f25af0434bd703c502c230f5a5ff0238b6b8a3",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd pulseaudio/
./autogen.sh LDFLAGS=-L$prefix/lib CFLAGS=-I$prefix/include --prefix=$prefix --host=$target --without-caps --disable-manpages
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpulse", :libpulse),
    ExecutableProduct(prefix, "pasuspender", :pasuspender),
    ExecutableProduct(prefix, "pacmd", :pacmd),
    ExecutableProduct(prefix, "pactl", :pactl),
    ExecutableProduct(prefix, "pulseaudio", :pulseaudio),
    LibraryProduct(prefix, "libpulse-simple", Symbol("libpulse-simple")),
    ExecutableProduct(prefix, "pacat", :pacat)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/yukota/libsndfileBuilder/releases/download/v1.0.28-0/build_libsndfile.v1.0.28.jl",
    "https://github.com/yukota/libatomic_opsBuilder/releases/download/v7.6.8-0/build_libatomic_ops.v7.6.8.jl",
    "https://github.com/yukota/SpeexDSPBuilder/releases/download/v1.2.0-rc3-0/build_SpeexDSP.v1.2.0-rc3.jl",
    "https://github.com/yukota/LibtoolBuilder/releases/download/v2.4.6-1/build_Libtool.v2.4.6.jl",
    "https://github.com/yukota/json-cBuilder/releases/download/v0.13.1-0/build_json-c.v0.13.1.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Gettext-v0.19.8-0/build_Gettext.v0.19.8.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Libiconv-v1.15-0/build_Libiconv.v1.15.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: yukota/PulseAudioBuilder

max_issues_repo_name: yukota/PulseAudioBuilder

max_forks_repo_name: yukota/PulseAudioBuilder


text:
using IteratorInterfaceExtensions
using Test

struct MyType
end

@testset "IteratorInterfaceExtensions" begin

@test isiterable(MyType()) == false
@test isiterable([1,2,3]) == true

@test_throws ErrorException getiterator(MyType())

@test [1,2,3] == getiterator([1,2,3])

@test IteratorSize2([1,2,3]) == Base.HasShape{1}()

end


text:
# Figure out which AD backend to test
const AD = get(ENV, "AD", "All")

function test_ad(f, x, broken = (); rtol = 1e-6, atol = 1e-6)
    finitediff = FiniteDifferences.grad(central_fdm(5, 1), f, x)[1]

    if AD == "All" || AD == "Tracker"
        if :Tracker in broken
            @test_broken Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        else
            @test Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "ForwardDiff"
        if :ForwardDiff in broken
            @test_broken ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "Zygote"
        if :Zygote in broken
            @test_broken Zygote.gradient(f, x)[1] ≈ finitediff rtol=rtol atol=atol
        else
            @test Zygote.gradient(f, x)[1] ≈ finitediff rtol=rtol atol=atol
        end
    end

    if AD == "All" || AD == "ReverseDiff"
        if :ReverseDiff in broken
            @test_broken ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    return
end


text:
function filter_seqs_sites(seqs, filter_fun)
    filtered_sites = Vector{Vector{Char}}()
    observed_alleles = Set(Char[])
    for pos in 1:length(seqs[1])
        column = [seq[pos] for seq in seqs]
        if filter_fun(column)
            push!(filtered_sites, column)
        end
    end
    return filtered_sites
end

function filter_fasta_sites(in_fasta, filter_fun)
    names = String[]
    seqs = String[]
    read_fasta!(in_fasta, names, seqs);

    for seq in seqs
        @assert length(seq) == length(seqs[1])
    end

    filtered_sites = filter_seqs_sites(seqs, filter_fun)
    new_seqs = [Char[] for s in seqs]
    for site in filtered_sites
        for (i, allele) in enumerate(site)
            push!(new_seqs[i], allele)
        end
    end

    return [
        FASTA.Record(name, LongCharSeq(join(seq)))
        for (name, seq) in zip(names, new_seqs)
    ]
end

function get_variant_only_seqs(in_fasta)
    return filter_fasta_sites(in_fasta, is_variant_site)
end

function is_variant_site(column)
    observed_alleles = Set(Char[])
    fill_observed_alleles!(observed_alleles, column);
    return (length(observed_alleles) >= 2)
end

function fill_observed_alleles!(observed_alleles, column)
    for a in column
        if occursin(a, valid_alleles)
            push!(observed_alleles, a);
        elseif !occursin(a, skip_chars)
            throw(DomainError(a, string("$a not a valid allele")));
        end
    end
end

function variants_only_fasta(in_fasta, out_fasta)
    variant_only_seqs = get_variant_only_seqs(in_fasta);
    w = FASTA.Writer(out_fasta)
    for record in variant_only_seqs
        write(w, record)
    end
end

function is_core_site(column)
    for a in column
        if a=='N' || a=='-'
            return false
        end
    end
    return true
end

function get_core_only_seqs(in_fasta)
    return filter_fasta_sites(in_fasta, is_core_site)
end

export filter_fasta_sites
export variants_only_fasta
export get_variant_only_seqs
export get_core_only_seqs


text:
__precompile__()
module GLPKMathProgInterface

export
    GLPKInterfaceLP,
    GLPKInterfaceMIP,
    GLPKSolverLP,
    GLPKSolverMIP

include("GLPKInterfaceBase.jl")
include("GLPKInterfaceLP.jl")
include("GLPKInterfaceMIP.jl")

const GLPKSolverLP = GLPKInterfaceLP.GLPKSolverLP
const GLPKSolverMIP = GLPKInterfaceMIP.GLPKSolverMIP

# Enables GLPK to act as a conic solver
import MathProgBase
MathProgBase.ConicModel(s::Union{GLPKSolverLP,GLPKSolverMIP}) = MathProgBase.LPQPtoConicBridge(MathProgBase.LinearQuadraticModel(s))
MathProgBase.supportedcones(::Union{GLPKSolverLP,GLPKSolverMIP}) = [:Free,:Zero,:NonNeg,:NonPos]

end


text:
using BioMedQuery.UMLS
using BioMedQuery.DBUtils

user= ""
psswd = ""
try
    user = ENV["UMLS_USER"]
    psswd = ENV["UMLS_PSSWD"]
catch
    println("UMLS tests require the following enviroment variables:")
    println("UMLS_USER")
    println("UMLS_PSSWD")
end

credentials = Credentials(user, psswd)
term = "obesity"
query = Dict("string"=>term, "searchType"=>"exact" )

@testset "Testing UMLS" begin
    tgt = get_tgt(credentials)
    all_results= search_umls(tgt, query)
    @test length(all_results[1]["result"]["results"]) == 2
    cui = best_match_cui(all_results)
    @test cui == "C0028754"
    sm = get_semantic_type(tgt, cui)
    @test sm[1] == "Disease or Syndrome"
end

@testset "Populate Semantic Network" begin

    host="localhost"
    username="root"
    password=""
    umls_sn_dbname="umls_sn_test"
    overwrite_db=true

    db_config = Dict(:host=>host,
                     :dbname=>umls_sn_dbname,
                     :username=>username,
                     :pswd=>password,
                     :overwrite=>overwrite_db)
    umls_sn_db = populate_net_mysql(db_config)

    #check collection of tables
    tables_query = BioMedQuery.DBUtils.select_all_tables(umls_sn_db)
    tables = ["SRDEF","SRFIL","SRFLD","SRSTR","SRSTRE1","SRSTRE2"]
    @test sort(tables) == sort(tables_query)

    db_query(umls_sn_db, "DROP DATABASE IF EXISTS $umls_sn_dbname;")

end


text:
module DungBase

using Dates, UUIDs, StaticArrays
import IntervalSets: width, (..), AbstractInterval, leftendpoint
using StatsBase, StructArrays, Dierckx, AngleBetweenVectors, LinearAlgebra

export VideoFile, WholeVideo, FragmentedVideo, DisjointVideo, AbstractTimeLine, AbstractPeriod, Instantaneous, Prolonged, Temporal, Board, Calibration, POI, Metadata, Run, Experiment, Track, Common, Point, PointCollection

export start, duration, stop, files, filenames, point, pointcollection, homing, searching, searchcenter, turningpoint

stop(x) = start(x) + duration(x)

include("videos.jl")

include("temporals.jl")

include("calibrations.jl")

include("pois.jl")

include("track.jl")

# include(joinpath(@__DIR__, "experimental_setup.jl"))

struct Metadata
    setup::Dict{Symbol, Any}
    comment::String
    date::DateTime
end

# Base.getproperty(x::Metadata, level::Symbol) = get(x.setup, level, missing)
#=getcoord(x::Run, poi::Symbol) = get(x.pois, poi, missing)
_getlevel(x::Metadata, level::Symbol) = get(x, level, missing)
getlevel(x::Run, level::Symbol) = _getlevel(x.metadata, level)=#

struct Run{T, M}
    data::T # pois::Dict{Symbol, T} or Track
    metadata::M
end

#=struct Run
    setup::Dict{Symbol, Any}
    pois::Dict{Symbol, POI}
    comment::String
    date::DateTime
end=#

struct Experiment
    runs::Vector{Run}
    description::String
end

end # module


text:
# -----------------------------------------------------------------------------
# Implementation of Atomistic MolecularDynamicsSimulator API
# -----------------------------------------------------------------------------

"""
    NBSimulator <: MolecularDynamicsSimulator

A wrapper around NBodySimulator to implement the Atomistic API.

**Type parameter descriptions**
- `S` the type of ODE solver
- `T<:Unitful.Time` the type for the time fields
- `C` the type for the thermostat field

**Field descriptions**
- `Δt::T` the time between timesteps, assumed to be in atomic units
- `steps::Integer` the number of timesteps for the simulation
- `t₀::T` the starting time of the simulation, assumed to be in atomic units;
    defaults to 0
- `thermostat::C` the thermostat for the simulation;
    many options are defined by `NBodySimulator`, but a user could also define a custom thermostat;
    defaults to the `NullThermostat`
- `simulator::S` the algorithm to be used for the ODE;
    defaults to VelocityVerlet
"""
struct NBSimulator{S<:OrdinaryDiffEqAlgorithm,T<:Unitful.Time,C<:Thermostat} <: MolecularDynamicsSimulator
    Δt::T
    steps::Int
    t₀::T
    thermostat::C
    simulator::S
end
function NBSimulator(Δt::T, steps::Int; t₀::Real=zero(T), thermostat::C=NullThermostat(), simulator::S=NBodySimulator.VelocityVerlet()) where {S<:OrdinaryDiffEqAlgorithm,T<:Real,C<:Thermostat}
    Δt, t₀ = promote(Δt * TIME_UNIT, t₀ * TIME_UNIT)
    NBSimulator(Δt, steps, t₀, thermostat, simulator)
end
function NBSimulator(Δt::T, steps::Integer; t₀::Unitful.Time=zero(T), thermostat::C=NullThermostat(), simulator::S=NBodySimulator.VelocityVerlet()) where {S<:OrdinaryDiffEqAlgorithm,T<:Unitful.Time,C<:Thermostat}
    Δt, t₀ = promote(Δt, t₀)
    NBSimulator(Δt, steps, t₀, thermostat, simulator)
end

# Extract the tuple of start_time, end_time from the simulator
time_range(simulator::NBSimulator) = Float64.(austrip.((simulator.t₀, simulator.t₀ + simulator.steps * simulator.Δt)))

function simulate(system::AbstractSystem{3}, simulator::NBSimulator, potential::AbstractPotential)
    wrapper = InteratomicPotentialParameters(potential)
    potentials = Dict{Symbol,PotentialParameters}(:custom => wrapper)
    nb_system = PotentialNBodySystem{ElementMassBody}(ElementMassBody.(system), potentials)
    simulation = NBodySimulation(nb_system, time_range(simulator), nbs_boundary_conditions(system), simulator.thermostat, austrip(u"k"))
    result = run_simulation(simulation, simulator.simulator, dt=austrip(simulator.Δt))
    NBSResult(result, wrapper.energy_cache)
end

# -----------------------------------------------------------------------------
# Integration with InteratomicPotentials
# -----------------------------------------------------------------------------

# Internal struct that represents a potential in the format accepted by NBodySimulator
# Wraps the underlying InteratomicPotential with a cache of the forces for the current timestep
# Also stores the potential energy from every timestep to be referenceable after the simulation
struct InteratomicPotentialParameters{P<:AbstractPotential} <: PotentialParameters
    potential::P
    timestep_cache::Ref{Float64}
    force_cache::Ref{Vector{SVector{3,Float64}}}
    energy_cache::Vector{ENERGY_TYPE}
    InteratomicPotentialParameters(potential::AbstractPotential) = new{typeof(potential)}(potential, Ref{Float64}(), Ref{Vector{SVector{3,Float64}}}(), Vector{ENERGY_TYPE}())
end

function NBodySimulator.get_accelerating_function(parameters::InteratomicPotentialParameters, simulation::NBodySimulation)
    bodies = simulation.system.bodies
    boundary_conditions = simulation.boundary_conditions
    (dv, u, v, t, i) -> begin
        if !isassigned(parameters.timestep_cache) || t != parameters.timestep_cache[]
            # TODO: avoid copying data at every timestep
            particles = [AtomsBase.Atom(b, SVector{3}(r), SVector{3}(v), boundary_conditions) for (b, r, v) ∈ zip(bodies, eachcol(u), eachcol(v))]
            system = FlexibleSystem(particles, get_bounding_box(boundary_conditions), get_boundary_conditions(boundary_conditions))
            eandf = energy_and_force(system, parameters.potential)
            parameters.timestep_cache[] = t
            parameters.force_cache[] = [austrip.(f) for f ∈ eandf.f]
            push!(parameters.energy_cache, uconvert(ENERGY_UNIT, eandf.e))
        end
        dv .+= parameters.force_cache[][i] ./ bodies[i].m
    end
end


text:
using NEAT

paramsDict = NEAT.loadConfig(joinpath(dirname(@__FILE__),"XOR_config.txt"))
config = NEAT.Config(paramsDict)
g = NEAT.Global(config)

# set node gene type

# XOR-2
INPUTS = [0 0; 0 1; 1 0; 1 1]
OUTPUTS = [0, 1, 1, 0]

function eval_fitness(population)
    for chromo in population
        net = NEAT.createPhenotype(chromo)
        error = 0.0
        for i = 1:size(INPUTS,1)
#             NEAT.flush!(net) # not strictly necessary in feedforward nets
            output = NEAT.activate(net.nntype, net, [INPUTS[i,1],INPUTS[i,2]])
            error += (output[1] - OUTPUTS[i])^2
        end
        chromo.fitness = max(0, 1-sqrt(error/length(OUTPUTS)))
    end
end

println("Beginning XOR test...")
g.cf.nn_activation = :sigm
p = NEAT.Population(g)
p.evaluate = eval_fitness
NEAT.epoch(g, p, 300, false)

winner = p.best_fitness[end]
println(winner)
println("Number of evaluations: $(p.generation)   winner.id=$(winner.id)")

# Let's check if it's really solved the problem
println("Best network output:")
correct = falses(length(OUTPUTS))
brain = NEAT.createPhenotype(winner)
println("example\t\tcorrect\t\tpredicted")
for i = 1:size(INPUTS,1)
    output = NEAT.activate(brain.nntype, brain, [INPUTS[i,1],INPUTS[i,2]])
    @printf("%5d\t\t%1.5f \t%1.5f\n", i ,OUTPUTS[i], round(output[1],0))
    if OUTPUTS[i] == round(output[1],2) correct[i] = true end
end

# # Visualize the winner network (requires PyDot)
# #visualize.draw_net(winner) # best chromosome

# # Plots the evolution of the best/average fitness (requires Biggles)
# #visualize.plot_stats(pop.stats)
# # Visualizes speciation
# #visualize.plot_species(pop.species_log)

# # saves the winner
# #file = open('winner_chromosome', 'w')
# #pickle.dump(winner, file)
# #file.close()

max_stars_repo_name: Andy-P/NEAT

max_issues_repo_name: Andy-P/NEAT

max_forks_repo_name: Andy-P/NEAT


text:
#Demo 2: Reflection coefficients
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import Graphics: width, height

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)
μ0 = 4pi*1e-7 #F/m
ϵ0 = 8.854e-12 #H/m

#Input data
#-------------------------------------------------------------------------------
fmax = 10e9
fstep = 50e6
f = collect(0:fstep:fmax)
ℓ = 50e-3 #Meters
ZL = Float64[1e6, 377, 60.0]
_colors = [blue, red, green]

#==Equations
===============================================================================#
Γ(Z; Zref::Real=50.0) = (Z - Zref) ./ (Z + Zref)
#ZC: Characteristic impedance
#ZL: Load impendance (termination)
function Zline(ℓ::Real, f::Vector, ZL::Number; ZC::Number=50.0, α::Real=0, μ::Real=μ0, ϵ::Real=ϵ0)
	j = im
	β = f*(2pi*sqrt(μ*ϵ))
	γ = α+j*β
	tanh_γℓ = tanh(γ*ℓ)
	return ZC*(ZL+ZC*tanh_γℓ)./(ZC+ZL*tanh_γℓ)
end
function Γline(ℓ::Real, f::Vector, ZL::Number; ZC::Number=50.0, Zref::Number=50.0, α::Real=0, μ::Real=μ0, ϵ::Real=ϵ0)
	return Γ(Zline(ℓ, f, ZL; ZC=ZC, α=α, μ=μ, ϵ=ϵ), Zref=Zref)
end

#Calculations
#-------------------------------------------------------------------------------
Γload = []
for ZLi in ZL
	_Γ = Γline(ℓ, f, ZLi, ZC=40)
	push!(Γload, _Γ)
end

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Transmission Line Example")
mplot.ncolumns = 2

plot_linf = InspectDR.Plot2D()
	plot_linf.axes = InspectDR.axes(:lin, :dB20)
	plot_linf.ext_full = InspectDR.PExtents2D(ymax=5)
#	plot_linf.layout.legend.enabled=true
plot_logf = InspectDR.Plot2D()
	plot_logf.axes = InspectDR.axes(:log10, :dB20)
	plot_logf.ext_full = InspectDR.PExtents2D(xmin=10e6,ymax=5)
	plot_logf.layout.grid = grid(vmajor=true, vminor=true, hmajor=true)
plot_ysmith = InspectDR.Plot2D()
	plot_ysmith.axes = InspectDR.axes(:smith, :Y)
	plot_ysmith.ext_full = InspectDR.PExtents2D(xmin=-1.2,xmax=1.2,ymin=-1.2,ymax=1.2)
	plot_ysmith.layout.legend.enabled=true
plot_smith = InspectDR.Plot2D()
	plot_smith.axes = InspectDR.axes(:smith, :Z, ref=50)
	plot_smith.ext_full = InspectDR.PExtents2D(xmin=-1.2,xmax=1.2,ymin=-1.2,ymax=1.2)
	plot_smith.layout.legend.enabled=true

for plot in [plot_linf, plot_logf]
	a = plot.annotation
	a.title = "Reflection Coefficient (Γ)"
	a.xlabel = "Frequency (Hz)"
	a.ylabel = "Magnitude (dB)"
end

a = plot_ysmith.annotation
	a.title = "Y-Smith Chart"
	a.xlabel = "Real(Γ)"
	a.ylabel = "Imaginary(Γ)"

a = plot_smith.annotation
	a.title = "Z-Smith Chart"
	a.xlabel = "Real(Γ)"
	a.ylabel = "Imaginary(Γ)"

#Select which plots to actually display:
plotlist = [plot_linf, plot_logf, plot_ysmith, plot_smith]
#plotlist = [plot_smith]

for plot in plotlist
	for i in 1:length(Γload)
		wfrm = add(plot, f, Γload[i], id="ZL=$(ZL[i])")
		wfrm.line = line(color=_colors[i], width=2)
	end

	add(mplot, plot)
end

gplot = display(InspectDR.GtkDisplay(), mplot)

#==Save multi-plot to file
===============================================================================#

maximize_square = true
if maximize_square
	#Target plot size to get square Smith plots without gaps:
	lyt = plot_smith.layout
		lyt.wdata = 500
		lyt.hdata = 500
	bb = InspectDR.plotbounds(lyt, plot_smith.axes) #Required
		mplot.wplot = width(bb)
		mplot.hplot = height(bb)
end

InspectDR.write_png("export_multiplot.png", mplot)
InspectDR.write_svg("export_multiplot.svg", mplot)
InspectDR.write_eps("export_multiplot.eps", mplot)
InspectDR.write_pdf("export_multiplot.pdf", mplot)

:DONE


text:
using CSV
using DataFrames
using GZip
using Distributions
using NPZ
using FixedPointNumbers

function load_individuals(path::AbstractString)::DataFrame
  df = GZip.open(path,"r") do io
    df = CSV.read(io, copycols=true)  # read CSV as DataFrame
    # create new DataFrame with
    df = DataFrame(
      age=Int8.(df.age),
      gender = df.gender .== 1,
      household_index = Int32.(df.household_index),
      social_competence = Normed{UInt16,16}.(df.social_competence),
      ishealthcare = df.ishealthcare .== 1
    )
    sort!(df, :household_index)
  end
end

sample2dist(samples) = countuniquesorted(samples) |> x -> DiscreteNonParametric(x[1], x[2]./sum(x[2]) )
load_dist_from_samples(path::AbstractString) = npzread(path) |> sample2dist


text:
import Base: ==, setindex!, copy

readonly_vector_wrapper_functions = (:getindex, :firstindex, :lastindex, :iterate, :length, :size, :keys)

"""
A clause is a set of (non-conflicting) literals.
Literal identifiers are number; positive literals correspond to positive numbers, negative literals to negative numbers.

A set of literals is represented by a sorted vector.
It follows that negative literals precede positive literals.
"""
struct Clause
  literals::Vector{Int}

  Clause(literals::Vector{Int}) = new(sort(literals))
end

Clause(literals::Vararg{Int}) = Clause([literals...])

for f in readonly_vector_wrapper_functions
  @eval Base.$f(cl::Clause, args...) = Base.$f(cl.literals, args...)
end

function ==(c1::Clause, c2::Clause)
  c1.literals == c2.literals
end

"""
A formula is a set of clauses.

A set of clauses is represented by a vector.
Unlike `Clause`, no sorting is being done.

See also: [`Clause`](@ref).
"""
struct Formula
  variable_count::UInt
  clauses::Vector{Clause}
end

function determine_variable_count(clauses::Vector{Clause})::UInt
  maximum(clauses) do cl
    maximum(abs, cl)
  end
end

Formula(clauses::Vector{Clause}) = Formula(determine_variable_count(clauses),
                                           clauses)
Formula(clauses::Vararg{Clause}) = Formula([clauses...])
Formula(var_count) = Formula(var_count, [])

for f in readonly_vector_wrapper_functions
  @eval Base.$f(fla::Formula, args...) = Base.$f(fla.clauses, args...)
end

# TODO: Implement in subquadratic time.
function ==(f1::Formula, f2::Formula)
  for c1 in f1.clauses
    matched = false

    for c2 in f2.clauses
      if c1 == c2
        matched = true
        break
      end
    end

    if !matched
      return false
    end
  end

  return true
end

struct PartialAssignment
  assignment::Vector{Int8}
end

PartialAssignment(n::Number) = PartialAssignment(zeros(Int8, n))

for f in readonly_vector_wrapper_functions
  @eval Base.$f(pa::PartialAssignment, args...) = Base.$f(pa.assignment, args...)
end

function setindex!(pa::PartialAssignment, value::T, key) where {T <: Integer}
  @assert -1 <= value <= 1
  pa.assignment[key] = value
end

function ==(pa1::PartialAssignment, pa2::PartialAssignment)
  pa1.assignment == pa2.assignment
end

function copy(pa::PartialAssignment)
  PartialAssignment(copy(pa.assignment))
end

function is_clause_satisfied(cl::Clause, pa::PartialAssignment)
  any(cl) do l
    sign(l) == pa[abs(l)]
  end
end

function is_literal_undecided(l::Int, pa::PartialAssignment)::Bool
  0 == pa[abs(l)]
end

"""
    clause_size(cl::Clause, pa::PartialAssignment)

Counts the number of unsatisfied literals in `cl` with respect to partial
assignment `pa`. If `cl` is satisfied, returns 0.

See also: [`Clause`](@ref), [`PartialAssignment`](@ref).
"""
function clause_size(cl::Clause, pa::PartialAssignment)
  if is_clause_satisfied(cl, pa)
    return 0
  end

  count((l) -> is_literal_undecided(l, pa), cl)
end

function is_clause_unit(cl::Clause, pa::PartialAssignment)
  1 == clause_size(cl, pa)
end

"""
    is_clause_falsified(cl::Clause, pa::PartialAssignment)::Bool

Returns true if all literals of `cl` are decided and `cl` is not satisfied.
"""
function is_clause_falsified(cl::Clause, pa::PartialAssignment)::Bool
  !is_clause_satisfied(cl, pa) && 0 == clause_size(cl, pa)
end

"""
    find_unit_literal(cl::Clause, pa::PartialAssignment)

If `cl` is a unit clause, returns the “unit literal”, otherwise returns 0.

See also: [`PartialAssignment`](@ref)
"""
function find_unit_literal(cl::Clause, pa::PartialAssignment)::Int
  if !is_clause_unit(cl, pa)
    return 0
  end

  cl[findfirst((l) -> is_literal_undecided(l, pa), cl)]
end

function find_unit_literal(fla::Formula, pa::PartialAssignment)
  for cl in fla
    tmp = find_unit_literal(cl, pa)
    if 0 != tmp
      return tmp
    end
  end

  0
end

abstract type SatSolverResult end

struct SatisfiableFormula <: SatSolverResult
  satisfying_assignment::PartialAssignment
end

struct UnsatisfiableFormula <: SatSolverResult
end

function is_formula_satisfiable(res::SatisfiableFormula)
  true
end

function is_formula_satisfiable(res::UnsatisfiableFormula)
  false
end


text:
"""
# Population genetics analyses in Julia
Repository:    https://www.github.com/biojulia/PopGen.jl/
Documentation: https://biojulia.net/PopGen.jl/

\nA few things things you can do to get started:

## Import Data
- `PopGen.read(filename; kwargs...)`
- `genepop(infile; kwargs...)`  or similar file-specific importer
- use available `@gulfsharks` or `@nancycats` datasets

## Explore PopData
- `populations(PopData)` to view population information
- `loci(PopData)` to view locus names
- `samplenames(PopData)` to view sample names
- `missingdata(PopData, by = ...)` to view missing information

## Manipulate PopData
- `populations!(PopData, ...)` to rename populations
- `locations!(PopData, ...)` to add geographical coordinates
- `exclude!(PopData; kwargs...)` to selectively remove data

## Analyses
- `richness(PopData)` to calculate allelic richness
- `relatedness(PopData, method = ...)` to get pairwise relatedness of individuals
- `summary(PopData)` to calculate F-statistics, heterozygosity, etc.
- `hwetest(PopData)` to test for Hardy-Weinberg Equilibrium
- `pairwisefst(PopData)` to calculate FST between pairs of populations
"""
module PopGen

#   o O       o O       o O       o O       o O
# o | | O   o | | O   o | | O   o | | O   o | | O
# | | | | O | | | | Dependencies| | | | O | | | | O
# O | | o   O | | o   O | | o   O | | o   O | | o
#   O o       O o       O o       O o       O o

using Reexport
using PopGenCore
@reexport module PopGenCore
    export PopData, PopDataInfo, GenoArray, Genotype, SNP, Msat, genodata, sampleinfo, locusinfo
    export isbiallelic, ishom, ishet
    export delimited, csv, genepop, vcf, bcf, @nancycats, @gulfsharks
    export sampleinfo, locusinfo, samplenames, loci
    export copy, size, sort, dropmonomorphic, dropmonomorphic!
    export dropmultiallelic, dropmultiallelic!
    export loci, samples, locations, locations!
    export locidataframe, locimatrix
    export genotypes, genotypes, genotype
    export populations, populations!
    export exclude, remove, omit, exclude!, remove!, omit!, keep, keep!, filter, filter!
end
@reexport import PopGenCore: read, write

using Distributions, DataFrames, PooledArrays
using Random: shuffle
using ProgressMeter
using MultipleTesting, StatsBase

#   o O       o O       o O       o O       o O
# o | | O   o | | O   o | | O   o | | O   o | | O
# | | | | O | | | |Include Files| | | | O | | | | O
# O | | o   O | | o   O | | o   O | | o   O | | o
#   O o       O o       O o       O o       O o

include("Utils.jl")
# heterozygosity functions
include("Heterozygosity.jl")
# manipulation and exploration
include("DataExploration.jl")
# summary information
include("SummaryInfo.jl")
#Analyses
include("HardyWeinberg.jl")
include("FStatistics/FstGlobal.jl")
include("FStatistics/FstByLocus.jl")
include("FStatistics/PairwiseFST.jl")
include("FStatistics/FstPermutations.jl")
include("Relatedness/PairwiseRelatedness.jl")
include("Relatedness/RelatednessMoments.jl")
include("Relatedness/RelatednessPostHocs.jl")

end # module PopGen


text:
# power analysis functions for two sample t-test

"""
    powerT2nTest(d::Real = 0, n1::Real = 0, n2::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of an unbalanced two samples of `n1` and `n2` observations to test the effect size `d`
at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂, `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest, `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function powerT2nTest(;
    d::Real = 0.0,
    n1::Real = 0,
    n2::Real = 0,
    alpha::Float64 = 0.05,
    alternative::String = "two")

    check_args(d=d,n1=n1,n2=n2,alpha=alpha)

    if alternative == "less"
        tside = ttside = 1
    elseif alternative in ("two","twosided","two.sided","two-sided","two sided")
        tside = ttside = 2
        d = abs(d)
    elseif alternative == "greater"
        ttside = 3
        tside = 1
    end

    ν = ceil(Int64,n1 + n2 - 2)
    λ = d*(1/sqrt(1/n1 + 1/n2))
    if ttside == 1
        return cdf(NoncentralT(ν,λ),quantile(TDist(ν),alpha))
    elseif ttside == 2
        qu = cquantile(TDist(ν),alpha/tside)
        return ccdf(NoncentralT(ν,λ),qu) + cdf(NoncentralT(ν,λ),-qu)
    elseif ttside == 3
        return ccdf(NoncentralT(ν,λ),cquantile(TDist(ν),alpha))
    end
end

"""
    samplesizeT2nTest(d::Real = 0, n1::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size `n2` for a unbalanced two samples of which one has `n1` observations
to test the effect size `d` with power > `power` (default: 0.8) at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function samplesizeT2nTest(;
    n1::Real = 0,
    d::Real = 0.0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    alternative::String = "two")

    check_args(d=d,n1=n1,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerT2nTest(n1 = n1, n2 = x, d = d, alpha = alpha, alternative = alternative) - power, 2+1e-09, 1e+10))
end

"""
    effectsizeT2nTest(n1::Real = 0, n2::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that unbalanced two samples of `n1` and `n2` observations can detect
with power > `power` (default: 0.8) at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function effectsizeT2nTest(;
    n1::Real = 0,
    n2::Real = 0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(n1=n1,alpha=alpha,power=power)

    return fzero(x -> powerT2nTest(n1 = n1, n2 = n2, d = x, alpha = alpha, alternative = alternative) - power,.001,100)
end

"""
    alphaT2nTest(n1::Real = 0, n2::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that unbalanced two samples of `n1` and `n2` observations have
in detecting the effect size `d` with power > `power` (default: 0.8) at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function alphaT2nTest(;
    n1::Real = 0,
    n2::Real = 0,
    d::Real = 0.0,
    power::Float64 = 0.0,
    alternative::String = "two"
    )

    check_args(d=d,n1=n1,n2=n2,power=power)

    return fzero(x->powerT2nTest(n1 = n1, n2 = n2, d = d, alpha = x, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    pwr.T2nTest(d::Real = 0, n1::Real = 0, n2::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n2`),
effect size (`d`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function T2nTest(;
    n1::Real = 0,
    n2::Real = 0,
    d::Real = 0.0,
    alpha::Float64 = 0.0,
    power::Float64 = 0.0,
    alternative::String = "two")

    if sum([x == 0 for x in (n1,n2,d,alpha,power)]) != 1
        error("exactly one of `n2`, `d`, `power`, and `alpha` must be zero")
    end

    if power == 0.0
        power = powerT2nTest(n1 = n1, n2 = n2, d = d, alpha = alpha, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaT2nTest(n1 = n1, n2 = n2, d = d, power = power, alternative = alternative)
    elseif d == 0.0
        d = effectsizeT2nTest(n1 = n1, n2 = n2, alpha = alpha, power = power, alternative = alternative)
    elseif n2 == 0
        n2 = samplesizeT2nTest(n1 = n1, d = d, alpha = alpha, power = power, alternative = alternative)
    end

    stype = Dict("onesample" => "One-sample", "one.sample" => "One-sample", "one sample" => "One-sample",
        "twosample" => "Two-sample", "two.sample" => "Two-sample", "two sample" => "Two-sample",
        "paired" => "Paired")
    alt = Dict("two" => "two-sided", "two.sided" => "two-sided", "less" => "less", "greater" => "greater")
    note = ""

    return htest(
        string("T-test power calculation"),
        OrderedDict(
            "n1" => n1,
            "n2" => n2,
            "d" => d,
            "alpha" => alpha,
            "power" => power,
            "alternative" => alt[alternative]
            )
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr