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 = "Cairo"
version = v"1.16.0"

sources = [
    ArchiveSource("https://www.cairographics.org/releases/cairo-$(version).tar.xz",
                  "5e7b29b3f113ef870d1e3ecf8adf21f923396401604bda16d44be45e66052331"),
    DirectorySource("./bundled"),
]

version = v"1.16.1" # <-- This version number is a lie to build for the experimental platforms

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cairo-*/

if [[ "${target}" == *-mingw* ]]; then
    # Link against libssp to fix errors like
    #     /opt/x86_64-w64-mingw32/bin/../lib/gcc/x86_64-w64-mingw32/8.1.0/../../../../x86_64-w64-mingw32/bin/ld: .libs/cairo-output-stream.o: in function `memcpy':
    #     /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/string.h:202: undefined reference to `__memcpy_chk'
    #     /opt/x86_64-w64-mingw32/bin/../lib/gcc/x86_64-w64-mingw32/8.1.0/../../../../x86_64-w64-mingw32/bin/ld: .libs/cairo-win32-font.o: in function `memcpy':
    #     /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/string.h:202: undefined reference to `__memcpy_chk'
    #     /opt/x86_64-w64-mingw32/bin/../lib/gcc/x86_64-w64-mingw32/8.1.0/../../../../x86_64-w64-mingw32/bin/ld: .libs/cairo-pdf-interchange.o: in function `strcat':
    #     /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/string.h:234: undefined reference to `__strcat_chk'
    atomic_patch -p1 ../patches/mingw-libssp.patch
    autoreconf -fiv
fi

# Because `zlib` doesn't have a proper `.pc` file, configure fails to find.
export CPPFLAGS="-I${includedir}"

if [[ "${target}" == *-apple-* ]]; then
    BACKEND_OPTIONS="--enable-quartz --enable-quartz-image --disable-xcb --disable-xlib"
elif [[ "${target}" == *-mingw* ]]; then
    BACKEND_OPTIONS="--enable-win32 --disable-xcb --disable-xlib"
elif [[ "${target}" == *-linux-* ]] || [[ "${target}" == *freebsd* ]]; then
    BACKEND_OPTIONS="--enable-xlib --enable-xcb --enable-xlib-xcb"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --disable-static \
    --enable-ft \
    --enable-tee \
    --enable-svg \
    --enable-ps \
    --enable-pdf \
    --enable-gobject \
    --disable-dependency-tracking \
    ${BACKEND_OPTIONS}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libcairo-gobject", :libcairo_gobject),
    LibraryProduct("libcairo-script-interpreter", :libcairo_script_interpreter),
    LibraryProduct("libcairo", :libcairo),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"),
    Dependency("Glib_jll"),
    Dependency("Pixman_jll"),
    Dependency("libpng_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Bzip2_jll", v"1.0.8"; compat="1.0.8"),
    Dependency("Xorg_libXext_jll"),
    Dependency("Xorg_libXrender_jll"),
    Dependency("LZO_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8", julia_compat="1.6")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using RecipesBase
import ..Resampling: ResampledDatasetStatistic

@recipe function f(::Type{ResampledDatasetStatistic{T,B}},
        x::ResampledDatasetStatistic{T,B}) where {T, B}
    x.estimates
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 71a70a62-28fb-11eb-2092-7d3640b6ddf5
using Pkg, DrWatson

# ╔═╡ 71d2cdaa-28fb-11eb-3826-cf35b99cde1f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using SpecialFunctions
	using StatisticalRethinking
end

# ╔═╡ d31428e6-85d2-11eb-396b-c9d6da91e9ae
md" ### Fig-10.2s.jl"

# ╔═╡ 1a849306-88c0-11eb-1b19-85694099f38f
function alpha(sigma, beta)
	sigma * sqrt(gamma(1/beta) / gamma(3/beta))
end

# ╔═╡ 61e6d20a-88bf-11eb-32c3-216979f1c4fb
begin
	β = [1.0, 1.5, 2.0, 2.5]
	α = [alpha(1.0, β[i]) for i in 1:length(β)]
	σ = [α[i]^2 * gamma(3/β[i]) / gamma(1/β[i]) for i in 1:length(β)]

	x = -4.0:0.01:4.0
	fig1 = plot(xlims=(-4,0, 4.0), xlab="x", ylab="Density")
	for i in 1:length(β)
		y = pdf(PGeneralizedGaussian(0.0, α[i], β[i]), x)
		plot!(x, y, lab = "α=$(round(α[i], digits=1)), β[$i]=$(β[i])")
	end
	plot!(x, pdf(Normal(), x), line=(:dash, 2), color=:darkblue,
		lab = "Normal(0, 1)")
	fig1
end

# ╔═╡ 864c11f0-8900-11eb-1a83-87d3d8a6297d
begin
	x2 = -4.0:0.01:4.0
	β2 = 1.0:0.01:4.0
	y2 = []
	fig2 = plot(xlims=(1.0, 4.0), leg=false)
	for i in 1:length(β2)
		α2 = alpha(1.0, β2[i])
		σ2 = α2^2 * gamma(3/β2[i]) / gamma(1/β2[i])
		local p = pdf(PGeneralizedGaussian(0.0, α2, β2[i]), x)
		h2 = -sum([p[i] * log(p[i]) for i in 1:length(p)])
		append!(y2, [h2])
	end
	plot!(β2, y2, xlab = "Shape (β2)", ylab="Entropy")
	vline!([2.0], line=:dash)
	fig2
end

# ╔═╡ deb0ea70-8901-11eb-37c6-5dc0a9b2d0dd
plot(fig1, fig2, layout=(1,2))

# ╔═╡ 22d43054-8661-11eb-0dd7-fb846edb81c4
md" ### End Fig-10.2s.jl"

# ╔═╡ Cell order:
# ╟─d31428e6-85d2-11eb-396b-c9d6da91e9ae
# ╠═71a70a62-28fb-11eb-2092-7d3640b6ddf5
# ╠═71d2cdaa-28fb-11eb-3826-cf35b99cde1f
# ╠═1a849306-88c0-11eb-1b19-85694099f38f
# ╠═61e6d20a-88bf-11eb-32c3-216979f1c4fb
# ╠═864c11f0-8900-11eb-1a83-87d3d8a6297d
# ╠═deb0ea70-8901-11eb-37c6-5dc0a9b2d0dd
# ╟─22d43054-8661-11eb-0dd7-fb846edb81c4


text:
@testset "SBP elements" begin

    tol = 200*eps()

    N = 3
    f(N,r,s) = r^N + s^N # function for testing differentiation

    @testset "Hicken's 2N triangular SBP operators" begin
        rd = RefElemData(Tri(), SBP{Hicken}(), N)
        @unpack rq,sq,wq = rd
        rq2,sq2,wq2 = quad_nodes(Tri(),2*N)
        @test sum(wq.*f.(2*N,rq,sq)) ≈ sum(wq2.*f.(2*N,rq2,sq2))
        @test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
        @test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
        @test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
    end

    @testset "Kubatko's Legendre face node triangular SBP operators" begin
        rd = RefElemData(Tri(), SBP{Kubatko{LegendreFaceNodes}}(), N)
        @unpack rq,sq,wq = rd
        rq2,sq2,wq2 = quad_nodes(Tri(),2*N-1)
        @test sum(wq.*f.(2*N-1,rq,sq)) ≈ sum(wq2.*f.(2*N-1,rq2,sq2))
        @test rd.Dr*rd.r.^N ≈ N*rd.r.^(N-1)
        @test rd.Ds*rd.s.^N ≈ N*rd.s.^(N-1)
        @test norm(rd.Dr*rd.s + rd.Ds*rd.r) < tol
    end

    @testset "Warning for N=6 Kubatko Lobatto SBP nodes" begin
        @test_logs (:warn,"N=6 SBP operators with quadrature strength 2N-1 and Lobatto face nodes may require very small timesteps.") RefElemData(Tri(),SBP{Kubatko{LobattoFaceNodes}}(),6)
    end

    @testset "TensorProductLobatto quad" begin
        rd = RefElemData(Quad(),SBP(),N)
        @test rd.r == rd.rst[1]
        @test rd.Np == length(rd.r)
        @test rd.Nq == length(rd.rq)
        @test abs(sum(rd.wq)) ≈ 4
        @test abs(sum(rd.wf)) ≈ 8
        @test abs(sum(rd.wf .* rd.nrJ)) + abs(sum(rd.wf .* rd.nsJ)) < tol
        @test rd.Pq*rd.Vq ≈ I
        @test all(vec.(rd.rstf) .≈ (x->getindex(x,rd.Fmask)).(rd.rst))
        @test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
    end

    @testset "TensorProductLobatto Hex" begin
        rd = RefElemData(Hex(),SBP(),N)
        @test propertynames(rd)[1] == :elementType
        @test rd.t == rd.rst[3]
        @test rd.tf == rd.rstf[3]
        @test rd.tq == rd.rstq[3]
        @test rd.rp == rd.rstp[1]
        @test rd.sp == rd.rstp[2]
        @test rd.tp == rd.rstp[3]
        @test rd.Np == length(rd.r)
        @test rd.Nq == length(rd.rq)
        @test abs(sum(rd.wq)) ≈ 8
        @test abs(sum(rd.wf)) ≈ 6*4
        @test abs(sum(rd.wf .* rd.nrJ)) < tol
        @test abs(sum(rd.wf .* rd.nsJ)) < tol
        @test abs(sum(rd.wf .* rd.ntJ)) < tol
        @test rd.Pq*rd.Vq ≈ I
        @test invoke(inverse_trace_constant,Tuple{RefElemData},rd) ≈ inverse_trace_constant(rd)
    end

end


text:
@testset "IntrusiveList" begin

mutable struct MultiTagListNode
    value::Int64
    next_a::MultiTagListNode
    prev_a::MultiTagListNode
    next_b::MultiTagListNode
    prev_b::MultiTagListNode

    function MultiTagListNode(value)
        x = new(value)
        x.next_a = x
        x.prev_a = x
        x.next_b = x
        x.prev_b = x
        x
    end
end

IntrusiveCollections.setnext!(x::MultiTagListNode, next::MultiTagListNode, ::Val{:a}) = x.next_a = next
IntrusiveCollections.setnext!(x::MultiTagListNode, next::MultiTagListNode, ::Val{:b}) = x.next_b = next
IntrusiveCollections.setprev!(x::MultiTagListNode, prev::MultiTagListNode, ::Val{:a}) = x.prev_a = prev
IntrusiveCollections.setprev!(x::MultiTagListNode, prev::MultiTagListNode, ::Val{:b}) = x.prev_b = prev
IntrusiveCollections.getnext(x::MultiTagListNode, ::Val{:a}) = x.next_a
IntrusiveCollections.getnext(x::MultiTagListNode, ::Val{:b}) = x.next_b
IntrusiveCollections.getprev(x::MultiTagListNode, ::Val{:a}) = x.prev_a
IntrusiveCollections.getprev(x::MultiTagListNode, ::Val{:b}) = x.prev_b

list_types = [
    IntrusiveList{IntrusiveListNode{Int}},
    IntrusiveSList{IntrusiveSListNode{Int}}
]

@testset "push and pop $ListType" for ListType in list_types
    # TODO: empty!, append!, prepend!
    NodeType = eltype(ListType)
    list = ListType()
    n1 = NodeType(1)
    pushfirst!(list, n1)
    @test !isempty(list)
    @test all([x for x in list] .=== [n1])
    n2 = NodeType(2)
    pushfirst!(list, n2)
    @test !isempty(list)
    @test all([x for x in list] .=== [n2, n1])
    @test popfirst!(list) === n2
    @test !isempty(list)
    @test all([x for x in list] .=== [n1])
    @test popfirst!(list) === n1
    @test isempty(list)
    push!(list, n1)
    @test !isempty(list)
    @test all([x for x in list] .=== [n1])
    push!(list, n2)
    @test !isempty(list)
    @test all([x for x in list] .=== [n1, n2])
    @test pop!(list) === n2
    @test !isempty(list)
    @test all([x for x in list] .=== [n1])
    @test pop!(list) === n1
    @test isempty(list)
    @test_throws BoundsError pop!(list)
end

@testset "empty $ListType" for ListType in list_types
    list = ListType()
    @test isempty(list)
    push!(list, eltype(ListType)(1))
    @test !isempty(list)
    empty!(list)
    @test isempty(list)
end

@testset "length $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list = ListType()
    @test length(list) == 0
    push!(list, NodeType(1))
    @test length(list) == 1
    push!(list, NodeType(2))
    @test length(list) == 2
    pop!(list)
    @test length(list) == 1
end

@testset "delete $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list = ListType()
    n1 = NodeType(1)
    n2 = NodeType(2)
    n3 = NodeType(3)
    push!(list, n1)
    delete!(list, n1)
    @test isempty(list)
    push!(list, n1, n2)
    delete!(list, n1)
    @test all([x for x in list] .=== [n2])
    pushfirst!(list, n1)
    delete!(list, n2)
    @test all([x for x in list] .=== [n1])
    push!(list, n2, n3)
    delete!(list, n2)
    @test all([x for x in list] .=== [n1, n3])
end

@testset "deleteafter $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list = ListType()
    n1 = NodeType(1)
    n2 = NodeType(2)
    n3 = NodeType(3)
    push!(list, n1)
    deleteafter!(list, n1)
    @test isempty(list)
    push!(list, n1, n2)
    deleteafter!(list, n2)
    @test all([x for x in list] .=== [n2])
    pushfirst!(list, n1)
    deleteafter!(list, n1)
    @test all([x for x in list] .=== [n1])
    push!(list, n2, n3)
    deleteafter!(list, n1)
    @test all([x for x in list] .=== [n1, n3])
end

@testset "insertafter $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list = ListType()
    n1 = NodeType(1)
    n2 = NodeType(2)
    n3 = NodeType(3)
    push!(list, n1)
    insertafter!(list, n1, n2)
    @test all([x for x in list] .=== [n1, n2])
    insertafter!(list, n1, n3)
    @test all([x for x in list] .=== [n1, n3, n2])
end

@testset "append $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list1 = ListType()
    list2 = ListType()
    n1 = NodeType(1)
    n2 = NodeType(2)
    n3 = NodeType(3)
    n4 = NodeType(4)
    # empty -> empty
    append!(list1, list2)
    @test isempty(list1)
    @test isempty(list2)
    # empty -> [1]
    push!(list2, n1)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1])
    # empty -> [1, 2]
    list1 = ListType()
    list2 = ListType()
    push!(list2, n1, n2)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1, n2])
    # [1] -> empty
    list1 = ListType()
    list2 = ListType()
    push!(list1, n1)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1])
    # [1] -> [2]
    list1 = ListType()
    list2 = ListType()
    push!(list1, n1)
    push!(list2, n2)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1, n2])
    # [1, 2] -> [3]
    list1 = ListType()
    list2 = ListType()
    push!(list1, n1, n2)
    push!(list2, n3)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1, n2, n3])
    # [1, 2] -> [3, 4]
    list1 = ListType()
    list2 = ListType()
    push!(list1, n1, n2)
    push!(list2, n3, n4)
    append!(list1, list2)
    @test isempty(list2)
    @test all([x for x in list1] .=== [n1, n2, n3, n4])
end

@testset "accessors $ListType" for ListType in list_types
    # TODO: isempty, length, in, indexin, first, last
    NodeType = eltype(ListType)
    list = ListType()
    n1 = NodeType(1)
    n2 = NodeType(2)
    @test_throws BoundsError first(list)
    @test_throws BoundsError last(list)
    push!(list, n1)
    @test first(list) === n1
    @test last(list) === n1
    push!(list, n2)
    @test first(list) === n1
    @test last(list) === n2
end

@testset "iteration $ListType" for ListType in list_types
    NodeType = eltype(ListType)
    list = ListType()
    @test isempty(collect(list))
    @test typeof(collect(list)) == Vector{NodeType}
    pushfirst!(list, NodeType(1))
    pushfirst!(list, NodeType(2))
    pushfirst!(list, NodeType(3))
    @test [x.value for x in list] == [3, 2, 1]
end

@testset "tags $ListType" for ListType in (TaggedIntrusiveList, TaggedIntrusiveSList)
    a = ListType{MultiTagListNode, :a}()
    b = ListType{MultiTagListNode, :b}()
    n1 = MultiTagListNode(1)
    n2 = MultiTagListNode(2)
    n3 = MultiTagListNode(3)
    push!(a, n1)
    @test !isempty(a)
    @test isempty(b)
    @test all([x for x in a] .=== [n1])
    push!(b, n2)
    @test !isempty(a)
    @test !isempty(b)
    @test all([x for x in a] .=== [n1])
    @test all([x for x in b] .=== [n2])
    push!(a, n3)
    @test !isempty(a)
    @test !isempty(b)
    @test all([x for x in a] .=== [n1, n3])
    @test all([x for x in b] .=== [n2])
    @test pop!(b) === n2
    @test !isempty(a)
    @test isempty(b)
    @test all([x for x in a] .=== [n1, n3])
    @test pop!(a) === n3
    @test !isempty(a)
    @test isempty(b)
    @test all([x for x in a] .=== [n1])
    @test pop!(a) === n1
    @test isempty(a)
    @test isempty(b)
end

@testset "output" begin
    list = IntrusiveList{IntrusiveListNode{Int}}()
    @test "$list" == "[]"
    push!(list, IntrusiveListNode(1))
    @test "$list" == "[1]"
    push!(list, IntrusiveListNode(2))
    @test "$list" == "[1↔2]"
    push!(list, IntrusiveListNode(3))
    @test "$list" == "[1↔2↔3]"
    list = IntrusiveSList{IntrusiveSListNode{Int}}()
    @test "$list" == "[]"
    push!(list, IntrusiveSListNode(1))
    @test "$list" == "[1]"
    push!(list, IntrusiveSListNode(2))
    @test "$list" == "[1→2]"
    push!(list, IntrusiveSListNode(3))
    @test "$list" == "[1→2→3]"
end

end


text:
using ChainRulesCore: @thunk, ChainRulesCore, NO_FIELDS
using ArgCheck: @argcheck, @check
import OffsetArrays

# WCN order

struct TargetY end
struct TargetX end
struct TargetKer end

function argcheck_wcn_conv1d!(y, x, ker)
    @argcheck ndims(ker) == 3
    @argcheck ndims(y) == ndims(x) == 3
    y_axw, y_axc, y_axb = axes(y)
    x_axw, x_axc, x_axb = axes(x)
    ker_axw, ker_axc_x, ker_axc_y = axes(ker)
    @argcheck x_axb == y_axb
    @argcheck x_axc == ker_axc_x
    @argcheck y_axc == ker_axc_y
    @argcheck y_axw == first(x_axw)+last(ker_axw):last(x_axw)+first(ker_axw)
end

@noinline function wcn_conv1d!(y, x, ker, ::TargetY)
    y .= zero(eltype(y))
    for ib in axes(y, 3)
        for ic_x in axes(x, 2)
            for ic_y in axes(y, 2)
                # @inbounds for iw_ker in axes(ker, 1)
                for iw_ker in axes(ker, 1)
                    c = ker[iw_ker, ic_x, ic_y]
                    @simd for iw_y in axes(y, 1)
                        iw_x = iw_y - iw_ker
                        y[iw_y, ic_y, ib] += c * x[iw_x, ic_x, ib]
                    end
                end
            end
        end
    end
    return y
end

@noinline function wcn_conv1d!(y, x, ker, ::TargetX)
    x .= zero(eltype(x))
    for ib in axes(y, 3)
        for ic_x in axes(x, 2)
            for ic_y in axes(y, 2)
                # @inbounds for iw_ker in axes(ker, 1)
                for iw_ker in axes(ker, 1)
                    c = ker[iw_ker, ic_x, ic_y]
                    @simd for iw_y in axes(y, 1)
                        iw_x = iw_y - iw_ker
                        x[iw_x, ic_x, ib] += c * y[iw_y, ic_y, ib]
                    end
                end
            end
        end
    end
    return x
end

@noinline function wcn_conv1d!(y, x, ker, ::TargetKer)
    ker .= zero(eltype(ker))
    for ib in axes(y, 3)
        for ic_x in axes(x, 2)
            for ic_y in axes(y, 2)
                # for iw_ker in axes(ker, 1)
                @inbounds for iw_ker in axes(ker, 1)
                    @simd for iw_y in axes(y, 1)
                        iw_x = iw_y - iw_ker
                        ker[iw_ker, ic_x, ic_y] +=
                            x[iw_x, ic_x, ib] * y[iw_y, ic_y, ib]
                    end
                end
            end
        end
    end
    return ker
end

function calc_axes(y::Nothing, x, ker)
    axw_x = axes(x)[1]
    axw_ker = axes(ker)[1]
    axw = last(axw_ker)+first(axw_x):first(axw_ker)+last(axw_x)
    axc = axes(ker, 3)
    axb = axes(x, 3)
    return (axw, axc, axb)
end

function calc_axes(y, x::Nothing, ker)
    ay = axes(y)[1]
    ak = axes(ker)[1]
    ay_begin = first(ay)
    ay_end = last(ay)
    ak_begin = first(ak)
    ak_end = last(ak)
    # y[begin] = ker[end]+x[begin]
    # y[end] = ker[begin]+x[end]
    axw = (ay_begin - ak_end):(ay_end - ak_begin)
    axc = axes(ker, 2)
    axb = axes(y, 3)
    return (axw, axc, axb)
end

function calc_axes(y, x, ker::Nothing)
    ax = axes(x, 1)
    ay = axes(y, 1)
    # y[begin] = ker[end]+x[begin]
    # y[end] = ker[begin]+x[end]
    ay_begin = first(ay)
    ay_end = last(ay)
    ax_begin = first(ax)
    ax_end = last(ax)
    ak = (ay_end - ax_end):(ay_begin - ax_begin)
    return (ak, axes(x, 2), axes(y, 2))
end

target(y::Nothing,x,ker) = TargetY()
target(y,x::Nothing,ker) = TargetX()
target(y,x,ker::Nothing) = TargetKer()

function allocate(y, x, ker)
    template = something(y, x)
    alloc = similar(template, calc_axes(y, x, ker))
    return (something(y,alloc),
            something(x,alloc),
            something(ker, alloc))
end

function wcn_conv1d(y, x, ker)
    t = target(y,x,ker)
    y_, x_, ker_ = allocate(y, x, ker)
    argcheck_wcn_conv1d!(y_, x_, ker_)
    wcn_conv1d!(y_, x_, ker_, t)
end

conv(x, ker) = conv_hole(nothing, x, ker)
conv_hole(x, y, ker) = parent(wcn_conv1d(x,y,ker))

function loopsaxes_wcn_conv1d(y, x, ker)
    loopaxes = (axes(y, 3),
        axes(x, 2),
        axes(y, 2),
        axes(ker, 1),
        axes(y, 1)
    )
end

function nobs_wcn_conv1d(y, x, ker)
    2*mapreduce(length, *, loopsaxes_wcn_conv1d(y, x, ker))
end

################################################################################
##### AD
################################################################################
function ChainRulesCore.rrule(::typeof(conv), x, ker)
    y = conv(x, ker)
    function pullback_conv(ȳ)
        (NO_FIELDS,
            @thunk(conv_hole(ȳ, nothing, ker)),
            @thunk(conv_hole(ȳ, x, nothing)),
        )
    end
    (y, pullback_conv)
end


text:
struct Discretization{LType, UType, PType}
  low_mat::LType                        # Lower triangular matrix of LU composition
  upp_mat::UType                        # Upper triangular matrix of LU composition
  perm_array::PType                     # Permutation array of LU composition
  num_steps_year::Int64                 # Number of time steps per astronomical year
  mesh::Mesh
  model::Model                          # Physical model
  annual_temperature::Array{Float64, 2} # Annual temperature at every time step in an astronomical year
  rhs::Vector{Float64}
  last_rhs::Vector{Float64}
end

function Discretization(mesh, model, num_steps_year; run_garbage_collector = true)
    low_mat, upp_mat, perm_array = compute_lu_matrices(mesh, model, num_steps_year)

    if run_garbage_collector
      GC.gc()
    end

    annual_temperature = fill(5.0, mesh.dof, num_steps_year) # Magic initialization
    rhs     = zeros(mesh.dof)  # TODO: The EBM Fortran code initializes the RHS to zero... Maybe we want to initialize it differently
    last_rhs = zeros(mesh.dof)

    Discretization(low_mat, upp_mat, perm_array, num_steps_year, mesh, model, annual_temperature, rhs, last_rhs)
end

Base.size(discretization::Discretization) = size(discretization.mesh)

function Base.show(io::IO, discretization::Discretization)
  nx, ny = size(discretization)
  print(io, "Discretization() with ", nx, "×", ny, " degrees of freedom")
end

function compute_lu_matrices(mesh, model, num_steps_year)
  mat = compute_matrix(mesh,num_steps_year,model)
  lu_dec = lu(mat)
  low_mat = sparse(lu_dec.L)
  upp_mat = sparse(lu_dec.U)

  return low_mat, upp_mat, lu_dec.p
end


text:
import ..Datasets: dataset, datafile, reader, Mode, Test

snli_init() = register(DataDep(
    "GLUE-SNLI",
    """
    The Stanford Natural Language Inference (SNLI) task (GLUE version)
    """,
    "https://firebasestorage.googleapis.com/v0/b/mtl-sentence-representations.appspot.com/o/data%2FSNLI.zip?alt=media&token=4afcfbb2-ff0c-4b2d-a09a-dbf07926f4df",
    "48c972c3d3590cb79227cd91fda7319ac14068ce804e703364524e171b53dc16";
    post_fetch_method = fn -> begin
      mv(fn, "SNLI.zip")
      DataDeps.unpack("SNLI.zip")
      innerdir = "SNLI"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct SNLI <: Dataset end

function dataset(::Type{M}, d::SNLI) where M <: Mode
    ds = reader(datafile(M, d))
    header = split(take!(ds), '\t')
    field_num = length(header)
    needed_field = (8,9,11)
    rds = get_channels(String, length(needed_field))

    task = @async begin
        for s ∈ ds
            s = split(s, '\t')
            length(s) != field_num && continue
            for (i, j) ∈ enumerate(needed_field)
                put!(rds[i], s[j])
            end
        end
    end
    for rd ∈ rds
        bind(rd, task)
    end

    rds
end

function trainfile(::SNLI)
    datadep"GLUE-SNLI/train.tsv"
end

function devfile(::SNLI)
    datadep"GLUE-SNLI/dev.tsv"
end

function testfile(::SNLI)
    datadep"GLUE-SNLI/test.tsv"
end

get_labels(::SNLI) = ("entailment", "neutral", "contradiction")


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra

np = pyimport("numpy")
nx, ny =  8,10
nnodes, neles = (nx + 1)*(ny + 1), nx*ny
x = np.linspace(0.0, 0.5, nx + 1)
y = np.linspace(0.0, 0.5, ny + 1)
X, Y = np.meshgrid(x, y)
nodes = zeros(nnodes,2)
nodes[:,1], nodes[:,2] = X'[:], Y'[:]
ndofs = 2

EBC, g = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)

EBC[collect(1:nx+1), :] .= -1

function ggt(t)
    return zeros(sum(EBC.==-2)), zeros(sum(EBC.==-2))
end
gt = ggt

#pull in the y direction
FBC, fext = zeros(Int64, nnodes, ndofs), zeros(nnodes, ndofs)

# todo PARAMETER
FORCE_TYPE = "nonconstant"

if FORCE_TYPE == "constant"
    #pull in the y direction
    FBC[collect((nx+1)*ny + 1:(nx+1)*ny + nx+1), 2] .= -1
    fext[collect((nx+1)*ny + 1:(nx+1)*ny + nx+1), 2] = collect(range(2.0, stop=5.0, length=nx+1))*1e8*(0.1*tid+0.5)
    fext[(nx+1)*ny + 1, 2] /= 2.0
    fext[(nx+1)*ny + nx+1, 2] /= 2.0
else
    FBC[collect((nx+1)*ny + 1:(nx+1)*ny + nx+1), 2] .= -2
end

#force load function
function fft(t)
    local f
    if (@isdefined data_for) && data_for=="elasticity"
        @info "Generating data using small force"
        f = 1.0e6 *(1.0*tid) * sin(pi*t/T) * ones(nx + 1) # used for learning elasticity H
    else
        f = 1.0e8 *(1.0*tid) * sin(pi*t/T) * ones(nx + 1)
    end
    f[1] /= 2.0
    f[end] /= 2.0
    return f
end
ft = fft

elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx+1)*(j-1) + i
        elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        coords = nodes[elnodes,:]
        push!(elements,SmallStrainContinuum(coords,elnodes, prop,2))
    end
end

T = 0.0005
NT = 100
Δt = T/NT
stress_scale = 1.0e10
strain_scale = 1


text:
using Distributions
using HypothesisTests
using Random
using Plots

function create_recomb_sites(l::Float64, t::Int64)
  recomb_sites::Vector{Float64} = [0.0];
  if !(t == 0)
    while (site = recomb_sites[end] + rand(Exponential(1/t))) < l
      push!(recomb_sites, site);
    end
  end
  push!(recomb_sites, l);
  return recomb_sites
end

function create_ancestor(l::Float64, t::Int64, pₑ::Float64, pₐ::Float64, pₙ::Float64)

  recomb_sites = create_recomb_sites(l, t);
  pops::Vector{String} = wsample(["european", "african", "native"], [pₑ, pₐ, pₙ], length(recomb_sites)-1);
  return recomb_sites, pops
end

function admix!(recomb_offspring::Vector{Float64}, pops_offspring::Vector{String}, recomb_ancestors::Array{Array{Float64,1},1}, pops_ancestors::Array{Array{String,1},1}, new_recomb_sites::Vector{Float64}, start₁::Int64, start₂::Int64)
  if length(new_recomb_sites) == 1
    append!(recomb_offspring, recomb_ancestors[1][start₁:end]);
    append!(pops_offspring, pops_ancestors[1][max((start₁-1),1):end]);
  else
    let ind₁ = searchsortedlast(recomb_ancestors[1], new_recomb_sites[1]), ind₂ = searchsortedfirst(recomb_ancestors[2], new_recomb_sites[1])
      append!(recomb_offspring, recomb_ancestors[1][start₁:ind₁]);
      push!(recomb_offspring, new_recomb_sites[1]);
      append!(pops_offspring, pops_ancestors[1][max((start₁-1),1):ind₁]);
      admix!(recomb_offspring, pops_offspring, recomb_ancestors[[2,1]], pops_ancestors[[2,1]], new_recomb_sites[2:end], ind₂, ind₁ + 1);
    end
  end
end

function create_admix(recomb_ancestors::Array{Array{Float64,1},1}, pops_ancestors::Array{Array{String,1},1})
  recomb_offspring::Vector{Float64} = [];
  pops_offspring::Vector{String} = [];
  let start₁ = 1, start₂ = 1, l = recomb_ancestors[1][end], flip_coin = shuffle(1:2)
    new_recomb_sites = create_recomb_sites(l, 1);
    admix!(recomb_offspring, pops_offspring, recomb_ancestors[flip_coin], pops_ancestors[flip_coin], new_recomb_sites[2:end], 1, 1);
  end
  return recomb_offspring, pops_offspring
end

function get_hap_lengths(recomb_sites::Vector{Float64}, pops::Vector{String})
  hap_lengths = Vector{Float64}(undef, 0);
  hap_pops = Vector{String}(undef, 0);
  let ind = 1, l = length(recomb_sites)
    while ind < l
      j = 0;
      while ((ind + j + 1) <= (l-1))  && (pops[ind+j+1] == pops[ind])
        j +=1;
      end
      push!(hap_lengths, recomb_sites[ind + j + 1] - recomb_sites[ind]);
      push!(hap_pops, pops[ind]);
      ind = ind + j + 1;
    end
  end
  return hap_lengths, hap_pops
end

function family_tree_admix!(recomb_tree::Vector{Vector{Float64}}, tree_pop::Vector{Vector{String}})
  len = length(recomb_tree);
  if !(len == 1)
    for ind = 1:div(len,2)
      let son = create_admix(collect((pop!(recomb_tree), pop!(recomb_tree))), collect((pop!(tree_pop), pop!(tree_pop))))
        pushfirst!(recomb_tree, son[1]);
        pushfirst!(tree_pop, son[2]);
      end
    end
    family_tree_admix!(recomb_tree, tree_pop);
  end
end

function family_tree(l::Float64, t::Int64, t_ancestors::Vector{Int64}, pₑ::Vector{Float64}, pₐ::Vector{Float64}, pₙ::Vector{Float64})
  recomb_tree = Array{Array{Float64, 1}, 1}(undef, 2^t);
  tree_pop = Array{Array{String, 1}, 1}(undef, 2^t);
  for ind = 1:2^t
    let ancestor = create_ancestor(l, t_ancestors[ind], pₑ[ind], pₐ[ind], pₙ[ind])
      recomb_tree[ind] = ancestor[1];
      tree_pop[ind] = ancestor[2];
    end
  end
  family_tree_admix!(recomb_tree, tree_pop);
  return pop!(recomb_tree), pop!(tree_pop)
end

function max_hap(recombs::Vector{Float64}, pops::Vector{String}, target_pop::String = "native")
  hap_lengths, hap_pops = get_hap_lengths(recombs, pops);
  return sum(hap_pops .== target_pop) == 0 ? 0.0 : maximum(hap_lengths[hap_pops .== target_pop])
end

function simulate_cdf_max_hap(max_iter::Int64, l::Float64, t::Int64, t_ancestors::Vector{Int64}, pₑ::Vector{Float64}, pₐ::Vector{Float64}, pₙ::Vector{Float64})
  sim = Vector{Float64}(undef,max_iter);
  for ind = 1:max_iter
    sim[ind] = max_hap(family_tree(l, t, t_ancestors, pₑ, pₐ, pₙ)...)
  end
  return sim
end

function sum_hap(recombs::Vector{Float64}, pops::Vector{String}, target_pop::String = "native")
  hap_lengths, hap_pops = get_hap_lengths(recombs, pops);
  return sum(hap_pops .== target_pop) == 0 ? 0.0 : sum(hap_lengths[hap_pops .== target_pop])
end

function simulate_cdf_sum_hap(max_iter::Int64, l::Float64, t::Int64, t_ancestors::Vector{Int64}, pₑ::Vector{Float64}, pₐ::Vector{Float64}, pₙ::Vector{Float64})
  sim = Vector{Float64}(undef, max_iter);
  for ind = 1:max_iter
    sim[ind] = sum_hap(family_tree(l, t, t_ancestors, pₑ, pₐ, pₙ)...);
  end
  return sim
end

####

function get_atoms(sim::Vector{Float64}, l::Float64)::Vector{Float64}
  N = length(sim)
  p0 = sum(sim .== 0.)/N
  pL = sum(sim .== l)/N
  return [p0, pL]
end

function get_p(sim::Vector{Float64}, mi::Float64)::Float64
  N = length(sim)
  return sum(sim .≤ mi)/N
end

function get_pᵢ(p::Float64, atoms::Vector{Float64})
  return (p-atoms[1])/(1.0 - atoms[1])
end

###
# Functions that simulate statistics under the null hypothesis,
# for both maximum and sum statistics
# p_e, p_a and p_n are vectors that represent average percentage of genetic
# information of each ancestral population

function simulate_max_all(max_iter::Int64, L::Vector{Float64}, t::Int64, t_ancestors::Vector{Int64},
        pₑ::Vector{Float64}, pₐ::Vector{Float64}, pₙ::Vector{Float64})::Matrix{Float64}
    sim_all = Matrix{Float64}(undef, 22, max_iter)
    for i = 1:22
      sim_all[i, :] = simulate_cdf_max_hap(max_iter, L[i], t, t_ancestors, pₑ, pₐ, pₙ)
    end
    return sim_all
end

function simulate_sum_all(max_iter::Int64, L::Vector{Float64}, t::Int64, t_ancestors::Vector{Int64},
        pₑ::Vector{Float64}, pₐ::Vector{Float64}, pₙ::Vector{Float64})::Matrix{Float64}
    sim_all = Matrix{Float64}(undef, 22, max_iter)
    for i = 1:22
      sim_all[i, :] = simulate_cdf_sum_hap(max_iter, L[i], t, t_ancestors, pₑ, pₐ, pₙ)
    end
    return sim_all
end

###

function atoms_matrix(sim_all::Matrix{Float64}, L::Vector{Float64})
   atoms_all = Matrix{Float64}(undef, 22, 2)
   for i = 1:22
     atoms_all[i, :] = get_atoms(sim_all[i, :], L[i])
   end
   return atoms_all
 end

function vector_pᵢ(m_obs::Vector{Float64}, sim_all::Matrix{Float64},
  atoms_all::Matrix{Float64})::Vector{Float64}
  p_all = Vector{Float64}(undef, 22)
  for i = 1:22
    p_all[i] = get_p(sim_all[i, :], m_obs[i])
  end
  pi_all = (p_all - atoms_all[:, 1])./(1.0 .- atoms_all[:, 1])
  return pi_all
end

###
# Functions to obtain the test p-value from all chromosome scores

function get_pmax(pi_all::Vector{Float64}, atoms_all::Matrix{Float64})::Float64
  pm = maximum(pi_all)
  pvalue = 1.
  for i = 1:22
    pvalue *= atoms_all[i,1] + pm * (1.0 - atoms_all[i,1] - atoms_all[i,2])
  end
  return pvalue
end

function get_psum(pi_all::Vector{Float64}, atoms_all::Matrix{Float64})::Float64
  function dist_pᵢ(u::Float64, atoms::Vector{Float64})::Float64
    if u < 1. - atoms[1] - atoms[2]
      return(u / (1. - atoms[1] - atoms[2]))
    elseif u < 1. - atoms[2]
      return 0.
    else
      return 1.
    end
  end

  N = 10000
  pm = sum(pi_all)
  sim = rand(22, N)
  for i = 1:22
    for j = 1:N
      sim[i,j] = dist_pᵢ(sim[i,j], atoms_all[i, :])
    end
  end
  sim_pi = vec(sum(sim, dims = 1))
  return sum(sim_pi .< pm)/N
end

###

function test_maxmax(m_obs::Vector{Float64}, max_iter::Int64, t::Int64,
   L::Vector{Float64})::Float64
  t_ancestors = fill(0, 2^t);
  pₑ = push!(fill(1.0, 2^t - 1), 0.0);
  pₐ = fill(0.0, 2^t);
  pₙ = push!(fill(0.0, 2^t - 1), 1.0);
  sim_all = simulate_max_all(max_iter, L, t, t_ancestors, pₑ, pₐ, pₙ)
  atoms_all = atoms_matrix(sim_all, L)
  pi_all = vector_pᵢ(m_obs, sim_all, atoms_all)
  pvalue = get_pmax(pi_all, atoms_all)
  return pvalue
end

function test_maxsum(m_obs::Vector{Float64}, max_iter::Int64, t::Int64,
   L::Vector{Float64})::Float64
  t_ancestors = fill(0, 2^t);
  pₑ = push!(fill(1.0, 2^t - 1), 0.0);
  pₐ = fill(0.0, 2^t);
  pₙ = push!(fill(0.0, 2^t - 1), 1.0);
  sim_all = simulate_max_all(max_iter, L, t, t_ancestors, pₑ, pₐ, pₙ)
  atoms_all = atoms_matrix(sim_all, L)
  pi_all = vector_pᵢ(m_obs, sim_all, atoms_all)
  pvalue = get_psum(pi_all, atoms_all)
  return pvalue
end

function test_summax(m_obs::Vector{Float64}, max_iter::Int64, t::Int64,
   L::Vector{Float64})::Float64
  t_ancestors = fill(0, 2^t);
  pₑ = push!(fill(1.0, 2^t - 1), 0.0);
  pₐ = fill(0.0, 2^t);
  pₙ = push!(fill(0.0, 2^t - 1), 1.0);
  sim_all = simulate_sum_all(max_iter, L, t, t_ancestors, pₑ, pₐ, pₙ)
  atoms_all = atoms_matrix(sim_all, L)
  pi_all = vector_pᵢ(m_obs, sim_all, atoms_all)
  pvalue = get_pmax(pi_all, atoms_all)
  return pvalue
end

function test_sumsum(m_obs::Vector{Float64}, max_iter::Int64, t::Int64,
   L::Vector{Float64})::Float64
  t_ancestors = fill(0, 2^t);
  pₑ = push!(fill(1.0, 2^t - 1), 0.0);
  pₐ = fill(0.0, 2^t);
  pₙ = push!(fill(0.0, 2^t - 1), 1.0);
  sim_all = simulate_sum_all(max_iter, L, t, t_ancestors, pₑ, pₐ, pₙ)
  atoms_all = atoms_matrix(sim_all, L)
  pi_all = vector_pᵢ(m_obs, sim_all, atoms_all)
  pvalue = get_psum(pi_all, atoms_all)
  return pvalue
end

##### test_method_all is the main function to run any test
# This function depends on all previous functions.
# method can take values in: test_maxmax, test_maxsum, test_summax, test_sumsum,
# depending on the choices of chromosome statistic and scores grouping
# Important note: T is a vector with every generation that we want to test;
# the first mating for generation T is trivial and thus not simulated,
# so if you want to test for two generations ago, then T = 1.
function test_method_all(M_obs::Matrix{Float64}, max_iter::Int64, T::Vector{Int64},
  L::Vector{Float64}, method::Function)::Matrix{Float64}
  N_donors = size(M_obs)[2]
  N_times = length(T)
  Pvalues = Matrix{Float64}(undef, N_donors, N_times)
  for t = 1:N_times
    for d = 1:N_donors
      Pvalues[d, t] = method(M_obs[:, d], max_iter, T[t], L)
    end
    print("t = ", T[t] + 1, " done \n")
  end
  return Pvalues
end

#####
## Importing the data
#####

using  DelimitedFiles
cd("/My/Work/Directory")

# est_max has maximum of all lengths for each chromosome
est_max = readdlm("est_max.csv", ' ', Float64);

# Now we take the maximum among pairs of chromosomes,
# we obtain a statistic for each chromosome pair
est_max_fin = Matrix{Float64}(undef, 22, 20);
for i = 1:20
  est_max_fin[:, i] = max.(est_max[2*i-1, :], est_max[2*i, :])
end

# Now, the same thing for the sum of all lengths for each chromosome
est_sum = readdlm("est_sum.csv", ' ', Float64);
est_sum_fin = Matrix{Float64}(undef, 22, 20);
for i = 1:20
  est_sum_fin[:, i] = max.(est_sum[2*i-1, :], est_sum[2*i, :])
end

# L will have the length of every chromosome pair (in centimorgans)
# We scale to work in morgans
L = vec(readdlm("largos_chr.csv", ' ', Float64)./100);

#####
# Plots for test p-values for every individual
res = test_method_all(est_max_fin, 10000, [1, 2,3, 4], L, test_maxmax)
scatter(res, ylim = [0,1], label = ["t = 2" "t = 3" "t = 4" "t = 5"], legend = :topleft,
  xlabel = "Individuals", ylabel = "p-value", mode = "markers", marker = [:diamond :star :circle :square])
plot!([1,20], [0.05, 0.05], width = 2, label = "p-value = 0.05")

#####
# Power of the test: first scenario

# t is number of generations minus one
t = 4;
# We set the structure of anccestors.
t_ancestors = fill(4, 2^t);

# we set the chromosome coloring probabilities
pₑ = fill(1.0 - 2.0/2^t, 2^t);
pₐ = fill(0.0, 2^t);
pₙ = fill(2.0/2^t, 2^t);

# We simulate observations under this scenario
n_obs = 100
sims_sum = Matrix{Float64}(undef, 22, n_obs)
for i = 1:22
  sims_sum[i,:] = simulate_cdf_sum_hap(n_obs, L[i], t, t_ancestors, pₑ, pₐ, pₙ);
end

# We compute p-values for the chosen methodology
max_iter = 1000;
res = test_method_all(sims_sum, max_iter, t, L, test_summax)
sum(res .< 0.05, dims = 1)/max_iter

#####
# Power of the test: second scenario

# sample of chromosomes that will not have native american genetic information
chr_bool = [fill(1, 11); fill(0, 11)][randperm(22)]

max_iter = 1000
res_maxmax = Matrix{Float64}(undef, max_iter, 4);
for t = 1:4
  t_ancestors = fill(0, 2^t);
  sims_max = Matrix{Float64}(undef, 22, max_iter)
  for i = 1:22
    if chr_bool[i] < 0.5
      pₑ = push!(fill(1.0 , 2^t - 1), 0.);
      pₐ = fill(0.0, 2^t);
      pₙ = push!(fill(0., 2^t - 1), 1.);
    else
      pₑ = fill(1.0 , 2^t);
      pₐ = fill(0.0, 2^t);
      pₙ = fill(0., 2^t);
    end
    sims_max[i,:] = simulate_cdf_max_hap(max_iter, L[i], t, t_ancestors, pₑ, pₐ, pₙ);
  end
  res_maxmax[:, t] = test_method_all(sims_max, 1000, [t], L, test_maxmax)
end

sum(res_maxmax .< 0.05, dims = 1)/1000

#####
# Power of the test: third scenario

max_iter = 1000
res_summax = Matrix{Float64}(undef, max_iter, 4);
for t = 2:5
  t_ancestors = fill(0, 2^t);
  pₑ = push!(fill(1.0 , 2^t - 1), 0.);
  pₐ = fill(0.0, 2^t);
  pₙ = push!(fill(0., 2^t - 1), 1.);

  sims_sum = Matrix{Float64}(undef, 22, max_iter)
  for i = 1:22
    sims_sum[i,:] = simulate_cdf_sum_hap(max_iter, L[i], t, t_ancestors, pₑ, pₐ, pₙ);
  end
  res_summax[:, t-1] = test_method_all(sims_sum, 1000, [t - 1], L, test_summax)
end

sum(res_summax .< 0.05, dims = 1)/1000


text:
using Diffusion
Diffusion.main()


text:
################################################################################
# Toy Example
################################################################################
using Plots
T = Float64

# Define the dynamics of the system
p = 3 # Number of players
model = DoubleIntegratorGame(p=p) # game with 3 players with unicycle dynamics
# model = UnicycleGame(p=p) # game with 3 players with unicycle dynamics
# model = BicycleGame(p=p) # game with 3 players with unicycle dynamics
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(1*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([2, 0.0,0,0]),
      SVector{model.ni[2],T}([2, 0.0,0,0]),
      SVector{model.ni[3],T}([2, 0.0,0,0]),
      # SVector{model.ni[4],T}([3,+0.8,0,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]

# # Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
# radius = 0.5*ones(p)
# μ = 4.0*ones(p)
# μ = 4.0*[1.0, 0.1]
# μ = 4.0*[0.1, 1.0]
# add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.19
add_collision_avoidance!(game_con, radius)

# # Add control bounds
# u_max =  5*ones(SVector{m,T})
# u_min = -5*ones(SVector{m,T})
# add_control_bound!(game_con, u_max, u_min)

# # Add wall constraint
# walls = [Wall([-1.0,-1.75], [1.0,0.25], [1.,-1.]/sqrt(2))]
# add_wall_constraint!(game_con, walls)

# # Add circle constraint
# xc = [1., 2., 3.]
# yc = [1., 2., 3.]
# radius = [0.1, 0.2, 0.3]
# add_circle_constraint!(game_con, xc, yc, radius)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.0, 0.0, 0.0, #0.0,
   -0.4, 0.0, 0.4, #0.6,
    0.0, 0.0, 0.0, #0.0,
    0.0, 0.0, 0.0, #0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-5
opts.α_dual = 1.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-5
opts.ϵ_sta = 1e-5
opts.ϵ_con = 1e-5
opts.ϵ_opt = 1e-5
opts.regularize = true
opts.αx_dual = ones(p)

# s = 0.5
# opts.αx_dual = [sqrt(s), 1/sqrt(s)]
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
@time newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

function kkt_residual(prob::GameProblem)
    residual!(prob, prob.pdtraj)
    col_res = collision_residual(prob)
    kkt_res = vcat(prob.core.res, col_res)
    return kkt_res
end

function collision_residual(prob::GameProblem)
    probsize = prob.probsize
    N = probsize.N
    p = probsize.p
    px = probsize.px
    n_col = Int(p*(p-1)/2)
    col_res = zeros((N-1)*n_col)

    off = 1
    for k = 2:N #fix for correct indexing
        for i = 1:p
            for j = i+1:p
                ind = findfirst(x -> all(x.con.x2 .== px[j]), prob.game_con.state_conval[i])
                conval = prob.game_con.state_conval[i][ind]
                col_res[off] = conval.vals[k-1][1]
                off += 1
            end
        end
    end
    return col_res
end

function kkt_residual_jacobian(prob::GameProblem)
    probsize = prob.probsize
    N = probsize.N
    p = probsize.p
    px = probsize.px
    n_col = Int(p*(p-1)/2)
    n_λ = p*(p-1)
    n_res = length(prob.core.res)

    residual_jacobian!(prob, prob.pdtraj)
    kkt_jac = spzeros(n_res + (N-1)*n_col, n_res + (N-1)*n_λ)
    kkt_jac[1:n_res, 1:n_res] .= prob.core.jac
    off = 1
    for k = 2:N
        for i = 1:p
            stamp = stampify(:opt, i, :x, 1, k)
            vind = prob.core.verti_inds[stamp]
            for j ∈ setdiff(1:p,i)
                ind = findfirst(x -> all(x.con.x2 .== px[j]), prob.game_con.state_conval[i])
                conval = prob.game_con.state_conval[i][ind]
                hind = n_res + off
                kkt_jac[vind, hind:hind] .= conval.jac[k-1]'
                off += 1
            end
        end
    end

    off = 1
    for k = 2:N
        for i = 1:p
            hind = prob.core.horiz_inds[:x][1][k]
            for j = i+1:p
                ind = findfirst(x -> all(x.con.x2 .== px[j]), prob.game_con.state_conval[i])
                conval = prob.game_con.state_conval[i][ind]
                vind = n_res + off
                kkt_jac[vind:vind, hind] .= conval.jac[k-1]
                off += 1
            end
        end
    end
    return kkt_jac
end

function update_λcol!(prob::GameProblem, Δλcol)
    probsize = prob.probsize
    N = probsize.N
    p = probsize.p
    px = probsize.px
    off = 1
    for k = 2:N
        for i = 1:p
            for j ∈ setdiff(1:p,i)
                ind = findfirst(x -> all(x.con.x2 .== px[j]), prob.game_con.state_conval[i])
                conval = prob.game_con.state_conval[i][ind]
                # @show off
                # @show "before" conval.λ[k-1][1]
                if conval.λ[k-1][1] != 0
                    conval.λ[k-1][1] += Δλcol[off]
                    conval.λ[k-1][1] = max(0, conval.λ[k-1][1])
                end
                # @show "after " conval.λ[k-1][1]
                off += 1
            end
        end
    end
    return nothing
end

function display_nullspace(prob_; M::Int=5, amplitude::T=1e-3, atol::T=1e-10) where {T}
    kkt_jac = kkt_residual_jacobian(prob_)
    kkt_jac = Matrix(kkt_jac)
    n_res = length(prob_.core.res)
    ns = nullspace(kkt_jac, atol=atol)
    n_vec = size(ns)[2]
    nv = ns[:,rand(1:n_vec)]
    nv = nv./mean(abs.(nv))*amplitude
    plt = plot(legend=false)
    anim = @animate for k = 1:M
        prob = deepcopy(prob_)
        Δtraj = nv[1:n_res]
        Δλcol = nv[n_res+1:end]
        set_traj!(prob.core, prob.Δpdtraj, k*Δtraj)
        # update_traj!(prob.pdtraj, prob.pdtraj, 1.0, prob.Δpdtraj)
        update_λcol!(prob, k*Δλcol)
        # @show Δλcol

        prob.opts.shift = 0
        prob.opts.α_dual = 0.0
        prob.opts.αx_dual = zeros(p)
        prob.opts.dual_reset = false
        prob.opts.ρ_0 = prob.pen.ρ[1]
        prob.opts.ρ_max = prob.pen.ρ[1]
        prob.opts.inner_print = false
        prob.opts.outer_print = false
        prob.opts.inner_iter = 20
        prob.opts.outer_iter = 3
        prob.opts.reg_0 = 1e-8
        newton_solve!(prob)
        residual!(prob)
        @show mean(abs.(prob.core.res))
        # @show mean(abs.(Δtraj))
        # @show mean(abs.(Δλcol))
        opt = mean(abs.(prob.stats.opt_vio[end].vio))
        sta = mean(abs.(prob.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
        println("opt = "*scn(opt)*"  ", opt <= prob.opts.ϵ_opt )
        println("sta = "*scn(sta)*"  ", opt <= prob.opts.ϵ_opt )
        println("dyn = "*scn(dyn)*"  ", opt <= prob.opts.ϵ_opt )
        # plot!(prob.stats)
        if opt < 1e-5
            plot!(prob.model, prob.pdtraj.pr, plt=plt)
        end
        # plot!(prob.stats)
    end
    path = joinpath("~/Documents/nullspace_p3_din.gif")
    gif(anim, path, fps = 10)
    return nothing
end

function follow_nullspace(prob; M::Int=5, amplitude::T=1e-3, atol::T=1e-10) where {T}

    plt = plot(legend=false)
    anim = @animate for k = 1:M
        Random.seed!(k)
        kkt_jac = kkt_residual_jacobian(prob)
        kkt_jac = Matrix(kkt_jac)
        n_res = length(prob.core.res)
        ns = nullspace(kkt_jac, atol=atol)
        n_vec = size(ns)[2]
        nv = ns[:,10]
        # nv = ns[:,rand(1:n_vec)]
        nv = nv./mean(abs.(nv))*amplitude

        Δtraj = nv[1:n_res]
        Δλcol = nv[n_res+1:end]
        set_traj!(prob.core, prob.Δpdtraj, Δtraj)
        update_traj!(prob.pdtraj, prob.pdtraj, 1.0, prob.Δpdtraj)
        update_λcol!(prob, Δλcol)

        prob.opts.shift = 0
        prob.opts.α_dual = 0.0
        prob.opts.αx_dual = zeros(p)
        prob.opts.dual_reset = false
        prob.opts.ρ_0 = prob.pen.ρ[1]
        prob.opts.ρ_max = prob.pen.ρ[1]
        prob.opts.inner_print = false
        prob.opts.outer_print = false
        prob.opts.inner_iter = 20
        prob.opts.outer_iter = 10
        prob.opts.reg_0 = 1e-8
        newton_solve!(prob)
        residual!(prob)
        println("***********************")
        opt = mean(abs.(prob.stats.opt_vio[end].vio))
        sta = mean(abs.(prob.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
        println("opt = "*scn(opt)*"  ", opt <= prob.opts.ϵ_opt )
        println("sta = "*scn(sta)*"  ", opt <= prob.opts.ϵ_opt )
        println("dyn = "*scn(dyn)*"  ", opt <= prob.opts.ϵ_opt )
        if opt < 1e-5
            plot!(prob.model, prob.pdtraj.pr, plt=plt)
        end
        # plot!(prob.stats)
    end
    path = joinpath("~/Documents/follow_nullspace_p3_uni.gif")
    gif(anim, path, fps = 20)
    return nothing
end

# evaluate!(prob.game_con, prob.pdtraj.pr)
# collision_residual(prob)
# kkt_residual(prob)
# kkt_jac = kkt_residual_jacobian(prob)
# ns = nullspace(Matrix(kkt_jac), atol=1e-10)
# nv = ns[:,1]

prob_copy = deepcopy(prob)
# display_nullspace(prob_copy, M=10, amplitude=1e-2)
follow_nullspace(prob_copy, M=10, amplitude=1e-4)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LoadBalancingRulePropertiesFormat <: SwaggerModel
    frontendIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendIPConfiguration
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    probe::Any # spec type: Union{ Nothing, SubResource } # spec name: probe
    protocol::Any # spec type: Union{ Nothing, TransportProtocol } # spec name: protocol
    loadDistribution::Any # spec type: Union{ Nothing, String } # spec name: loadDistribution
    frontendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: frontendPort
    backendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: backendPort
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    enableFloatingIP::Any # spec type: Union{ Nothing, Bool } # spec name: enableFloatingIP
    enableTcpReset::Any # spec type: Union{ Nothing, Bool } # spec name: enableTcpReset
    disableOutboundSnat::Any # spec type: Union{ Nothing, Bool } # spec name: disableOutboundSnat
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function LoadBalancingRulePropertiesFormat(;frontendIPConfiguration=nothing, backendAddressPool=nothing, probe=nothing, protocol=nothing, loadDistribution=nothing, frontendPort=nothing, backendPort=nothing, idleTimeoutInMinutes=nothing, enableFloatingIP=nothing, enableTcpReset=nothing, disableOutboundSnat=nothing, provisioningState=nothing)
        o = new()
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        setfield!(o, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("probe"), probe)
        setfield!(o, Symbol("probe"), probe)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("loadDistribution"), loadDistribution)
        setfield!(o, Symbol("loadDistribution"), loadDistribution)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("frontendPort"), frontendPort)
        setfield!(o, Symbol("frontendPort"), frontendPort)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("backendPort"), backendPort)
        setfield!(o, Symbol("backendPort"), backendPort)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("enableFloatingIP"), enableFloatingIP)
        setfield!(o, Symbol("enableFloatingIP"), enableFloatingIP)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("enableTcpReset"), enableTcpReset)
        setfield!(o, Symbol("enableTcpReset"), enableTcpReset)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("disableOutboundSnat"), disableOutboundSnat)
        setfield!(o, Symbol("disableOutboundSnat"), disableOutboundSnat)
        validate_property(LoadBalancingRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type LoadBalancingRulePropertiesFormat

const _property_map_LoadBalancingRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("frontendIPConfiguration")=>Symbol("frontendIPConfiguration"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("probe")=>Symbol("probe"), Symbol("protocol")=>Symbol("protocol"), Symbol("loadDistribution")=>Symbol("loadDistribution"), Symbol("frontendPort")=>Symbol("frontendPort"), Symbol("backendPort")=>Symbol("backendPort"), Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("enableFloatingIP")=>Symbol("enableFloatingIP"), Symbol("enableTcpReset")=>Symbol("enableTcpReset"), Symbol("disableOutboundSnat")=>Symbol("disableOutboundSnat"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_LoadBalancingRulePropertiesFormat = Dict{Symbol,String}(Symbol("frontendIPConfiguration")=>"SubResource", Symbol("backendAddressPool")=>"SubResource", Symbol("probe")=>"SubResource", Symbol("protocol")=>"TransportProtocol", Symbol("loadDistribution")=>"String", Symbol("frontendPort")=>"Int32", Symbol("backendPort")=>"Int32", Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("enableFloatingIP")=>"Bool", Symbol("enableTcpReset")=>"Bool", Symbol("disableOutboundSnat")=>"Bool", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ LoadBalancingRulePropertiesFormat }) = collect(keys(_property_map_LoadBalancingRulePropertiesFormat))
Swagger.property_type(::Type{ LoadBalancingRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancingRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ LoadBalancingRulePropertiesFormat }, property_name::Symbol) =  _property_map_LoadBalancingRulePropertiesFormat[property_name]

const _allowed_LoadBalancingRulePropertiesFormat_loadDistribution = ["Default", "SourceIP", "SourceIPProtocol"]

function check_required(o::LoadBalancingRulePropertiesFormat)
    (getproperty(o, Symbol("protocol")) === nothing) && (return false)
    (getproperty(o, Symbol("frontendPort")) === nothing) && (return false)
    true
end

function validate_property(::Type{ LoadBalancingRulePropertiesFormat }, name::Symbol, val)
    if name === Symbol("loadDistribution")
        Swagger.validate_param(name, "LoadBalancingRulePropertiesFormat", :enum, val, _allowed_LoadBalancingRulePropertiesFormat_loadDistribution)
    end
end


text:
using PairedComparisons
using Test
using DataFrames

@testset "PairedComparisons.jl" begin
    # Elo tests
    elo = Elo(kfac = 25.0) # A standard Elo with k of 25 and default rating of 1500

    @test predict(elo, 1, 2) == 0.5

    games1 = DataFrame(Player1 = [1], Player2 = [2], Player1Wins = [1], Player2Wins = [0])
    games2 = copy(games1)
    games2.Day = [1]

    fit!(elo, games1)

    @test length(elo.ratings) == 2
    @test abs(predict(elo, 1, 2) - 0.5359159269451023) < 0.001

    #Glicko tests
    glicko = Glicko(c = 85.0, default_rating = (1500.0, 350.0, -1))

    @test predict(glicko, 1, 2, 1) == 0.5

    fit!(glicko, games2)
    @test length(glicko.ratings) == 2
    @test abs(predict(glicko, 1, 2, 1) - 0.72720) < 0.001

    #BradleyTerry tests
    bt = BT()

    @test predict(bt, 1, 2) == 0.5
    add_games!(bt, games2)
    @test length(bt.playergames) == 2
    @test abs(loglikelihood(bt) + 2.0794) < 0.001
    iterate!(bt, 50)
    @test abs(loglikelihood(bt) + 1.6847756) < 0.001
    @test abs(predict(bt, 1, 2)  - 0.7419441266134) < 0.001

    #WHR tests

    whr = WHR()
    add_games!(whr, games2)
    @test length(whr.playerdaygames) == 2
    @test abs(predict(whr, 1, 2) - 0.5) < 0.0001
    @test abs(loglikelihood(whr) + 4.63934111119) < 0.001
    iterate!(whr, 5)
    @test abs(predict(whr, 1, 2) - 0.684114747) < 0.001

end


text:
"""
$(README)
"""
module Syncer

using Comonicon
using DocStringExtensions

"""
Simple `rsync` wrapper.

# Arguments

- `source`: source folder
- `dest`: destination folder

# Flags

- `--common`: only sync common folders
- `--notdry`: perform syncing
"""
@main function main(
        source::String,
        dest::String;
        common::Bool = false,
        notdry::Bool = false,
    )
    @assert !endswith(source, "/")
    @assert !endswith(dest, "/")

    cmd = `rsync --recursive --compress --archive --delete --verbose --progress --stats`

    dry = !notdry
    if dry
        cmd = `$cmd --dry-run`
    end

    if common
        # Source folders to sync over
        sourcefolders = basename.(filter(isdir, readdir(source; join = true)))
        destfolders = basename.(filter(isdir, readdir(dest; join = true)))
        commonfolders = intersect(sourcefolders, destfolders)

        println("Only in source:")
        for fld in setdiff(sourcefolders, commonfolders)
            println(joinpath(source, fld))
        end

        println("Only in dest:")
        for fld in setdiff(destfolders, commonfolders)
            println(joinpath(dest, fld))
        end

        srclist = tempname()
        open(srclist; write = true) do io
            foreach(commonfolders) do fld
                println(io, fld)
            end
        end
        cmd = `$cmd --files-from=$srclist`
    end

    cmd = `$cmd $source/ $dest/`

    if dry
        run(pipeline(cmd, `less`))
    else
        run(cmd)
    end

    return nothing
end

end # module


text:
function testAlgorithmA(m::Int64, C::Float64, p::Float64, ntrials::Int64)
  r = C*p
  v = r*(1-r^(m-1))/(1-r^m)
  testFunction(BernoulliFactory._algorithmA, ntrials, v, makef(p), m, C,
    GLOBAL_RNG)
end

function testHighPowerLogistic(m::Int64, β::Float64, C::Float64, p::Float64,
  ntrials::Int64)
  r = C*p
  v = (β*r)^m./sum((β*r).^(0:m))
  testFunction(BernoulliFactory._highPowerLogistic, ntrials, v, makef(p), m, β, C,
    GLOBAL_RNG)
end

function testLinearSmall(f::F, m::Int64, p::Float64, C::Float64) where
  F<:Function
  ϵ = (1.0 - C*p)
  @assert ϵ > 0.5
  ϵ = (0.5 + ϵ)/2.0
  testFunction(f, m, C*p, makef(p), C, ϵ, GLOBAL_RNG)
end

Random.seed!(192837465)

testAlgorithmA(2, 2.0, 0.2, ntrials)
testAlgorithmA(3, 2.0, 0.2, ntrials)
testAlgorithmA(4, 2.0, 0.2, ntrials)

testHighPowerLogistic(2, 2.5, 2.0, 0.2, ntrials)
testHighPowerLogistic(3, 4.5, 3.0, 0.4, ntrials)

testLinear(BernoulliFactory._huber2017, ntrials, 0.1, 2.0)
testLinear(BernoulliFactory._huber2017, ntrials, 0.4, 2.0)
testLinear(BernoulliFactory._huber2017, ntrials, 0.1, 3.0)
testLinear(BernoulliFactory._huber2017, ntrials, 0.6, 1.5)

testLinearSmall(BernoulliFactory._huber2017Small, ntrials, 0.1, 2.0)
testLinearSmall(BernoulliFactory._huber2017Small, ntrials, 0.2, 2.0)
testLinearSmall(BernoulliFactory._huber2017Small, ntrials, 0.1, 4.0)


text:
@testset "@q macro" begin
    @test 4 == @q 2 + 2
    @test [0, 1] == @q til(2)
    @test begin
        @q begin
          a = 2
          b = 42
        end
        @q(a) == 2 && @q(b) == 42
    end
    @test begin
        f = Dict(:x=>42)
        e = K(:(x+y))
        Q.resolve!(e, f)
        e == K(:(42+y))
    end
    @test (x = 42; 43 == @q let x;x+1 end)
    @test q(".z.K") == @q _z.K
end


text:
using MPI
using Test
using StaticArrays
using Logging, Printf

using ClimateMachine
using ClimateMachine.LinearSolvers
using ClimateMachine.GeneralizedConjugateResidualSolver
using ClimateMachine.GeneralizedMinimalResidualSolver
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGmethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.DGmethods
import ClimateMachine.DGmethods:
    BalanceLaw,
    vars_state_auxiliary,
    vars_state_conservative,
    vars_state_gradient,
    vars_state_gradient_flux,
    flux_first_order!,
    flux_second_order!,
    source!,
    boundary_state!,
    numerical_boundary_flux_second_order!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    init_state_auxiliary!,
    init_state_conservative!,
    LocalGeometry

import ClimateMachine.DGmethods.NumericalFluxes:
    NumericalFluxSecondOrder, numerical_flux_second_order!

if !@isdefined integration_testing
    const integration_testing = parse(
        Bool,
        lowercase(get(ENV, "JULIA_CLIMA_INTEGRATION_TESTING", "false")),
    )
end

struct PoissonModel{dim} <: BalanceLaw end

vars_state_auxiliary(::PoissonModel, T) = @vars(rhs_ϕ::T)
vars_state_conservative(::PoissonModel, T) = @vars(ϕ::T)
vars_state_gradient(::PoissonModel, T) = @vars(ϕ::T)
vars_state_gradient_flux(::PoissonModel, T) = @vars(∇ϕ::SVector{3, T})

boundary_state!(nf, bl::PoissonModel, _...) = nothing

function flux_first_order!(
    ::PoissonModel,
    flux::Grad,
    state::Vars,
    state_auxiliary::Vars,
    t::Real,
)
    nothing
end

function flux_second_order!(
    ::PoissonModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    state_auxiliary::Vars,
    t::Real,
)
    flux.ϕ = diffusive.∇ϕ
end

struct PenaltyNumFluxDiffusive <: NumericalFluxSecondOrder end

# There is no boundary since we are periodic
numerical_boundary_flux_second_order!(nf::PenaltyNumFluxDiffusive, _...) =
    nothing

function numerical_flux_second_order!(
    ::PenaltyNumFluxDiffusive,
    bl::PoissonModel,
    fluxᵀn::Vars{S},
    n::SVector,
    state⁻::Vars{S},
    diff⁻::Vars{D},
    hyperdiff⁻::Vars{HD},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    diff⁺::Vars{D},
    hyperdiff⁺::Vars{HD},
    aux⁺::Vars{A},
    t,
) where {S, HD, D, A}

    numerical_flux_second_order!(
        CentralNumericalFluxSecondOrder(),
        bl,
        fluxᵀn,
        n,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        t,
    )

    Fᵀn = parent(fluxᵀn)
    FT = eltype(Fᵀn)
    tau = FT(1)
    Fᵀn .-= tau * (parent(state⁻) - parent(state⁺))
end

function compute_gradient_argument!(
    ::PoissonModel,
    transformstate::Vars,
    state::Vars,
    state_auxiliary::Vars,
    t::Real,
)
    transformstate.ϕ = state.ϕ
end

function compute_gradient_flux!(
    ::PoissonModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    state_auxiliary::Vars,
    t::Real,
)
    diffusive.∇ϕ = ∇transform.ϕ
end

source!(::PoissonModel, _...) = nothing

# note, that the code assumes solutions with zero mean
sol1d(x) = sin(2pi * x)^4 - 3 / 8
dxx_sol1d(x) =
    -16 * pi^2 * sin(2pi * x)^2 * (sin(2pi * x)^2 - 3 * cos(2pi * x)^2)

function init_state_auxiliary!(
    ::PoissonModel{dim},
    aux::Vars,
    g::LocalGeometry,
) where {dim}
    aux.rhs_ϕ = 0
    @inbounds for d in 1:dim
        x1 = g.coord[d]
        x2 = g.coord[1 + mod(d, dim)]
        x3 = g.coord[1 + mod(d + 1, dim)]
        x23 = SVector(x2, x3)
        aux.rhs_ϕ -= dxx_sol1d(x1) * prod(sol1d, view(x23, 1:(dim - 1)))
    end
end

function init_state_conservative!(
    ::PoissonModel{dim},
    state::Vars,
    aux::Vars,
    coords,
    t,
) where {dim}
    state.ϕ = prod(sol1d, view(coords, 1:dim))
end

function run(
    mpicomm,
    ArrayType,
    FT,
    dim,
    polynomialorder,
    brickrange,
    periodicity,
    linmethod,
)

    topology = BrickTopology(mpicomm, brickrange, periodicity = periodicity)
    grid = DiscontinuousSpectralElementGrid(
        topology,
        polynomialorder = polynomialorder,
        DeviceArray = ArrayType,
        FloatType = FT,
    )

    dg = DGModel(
        PoissonModel{dim}(),
        grid,
        CentralNumericalFluxFirstOrder(),
        PenaltyNumFluxDiffusive(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    Qrhs = dg.state_auxiliary
    Qexact = init_ode_state(dg, FT(0))

    linearoperator!(y, x) = dg(y, x, nothing, 0; increment = false)

    linearsolver = linmethod(Q)

    iters = linearsolve!(linearoperator!, linearsolver, Q, Qrhs)

    error = euclidean_distance(Q, Qexact)

    @info @sprintf """Finished
    error = %.16e
    iters = %d
    """ error iters
    error, iters
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = 4
    base_num_elem = 4
    tol = 1e-9

    linmethods = (
        b -> GeneralizedConjugateResidual(3, b, rtol = tol),
        b -> GeneralizedMinimalResidual(b, M = 7, rtol = tol),
    )

    expected_result = Array{Float64}(undef, 2, 2, 3) # method, dim-1, lvl

    # GeneralizedConjugateResidual
    expected_result[1, 1, 1] = 5.0540243616448058e-02
    expected_result[1, 1, 2] = 1.4802275366044011e-03
    expected_result[1, 1, 3] = 3.3852821775121401e-05
    expected_result[1, 2, 1] = 1.4957957657736219e-02
    expected_result[1, 2, 2] = 4.7282369781541172e-04
    expected_result[1, 2, 3] = 1.4697449643351771e-05

    # GeneralizedMinimalResidual
    expected_result[2, 1, 1] = 5.0540243587512981e-02
    expected_result[2, 1, 2] = 1.4802275409186211e-03
    expected_result[2, 1, 3] = 3.3852820667079927e-05
    expected_result[2, 2, 1] = 1.4957957659220951e-02
    expected_result[2, 2, 2] = 4.7282369895963614e-04
    expected_result[2, 2, 3] = 1.4697449516628483e-05

    lvls = integration_testing ? size(expected_result)[end] : 1

    for (m, linmethod) in enumerate(linmethods), FT in (Float64,)
        result = Array{Tuple{FT, Int}}(undef, lvls)
        for dim in 2:3
            for l in 1:lvls
                Ne = ntuple(d -> 2^(l - 1) * base_num_elem, dim)
                brickrange =
                    ntuple(d -> range(FT(0), length = Ne[d], stop = 1), dim)
                periodicity = ntuple(d -> true, dim)

                @info (ArrayType, FT, m, dim)
                result[l] = run(
                    mpicomm,
                    ArrayType,
                    FT,
                    dim,
                    polynomialorder,
                    brickrange,
                    periodicity,
                    linmethod,
                )

                @test isapprox(
                    result[l][1],
                    FT(expected_result[m, dim - 1, l]),
                    rtol = sqrt(tol),
                )
            end

            if integration_testing
                @info begin
                    msg = ""
                    for l in 1:(lvls - 1)
                        rate = log2(result[l][1]) - log2(result[l + 1][1])
                        msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
                    end
                    msg
                end
            end
        end
    end
end

nothing


text:
using Interpolations: LinearInterpolation
import ..ContinuumAbsorption: total_continuum_absorption

"""
    synthesize(atm, linelist, λ_start, λ_stop, [λ_step=0.01]; metallicity=0, abundances=Dict(), vmic=0, ... )

Solve the transfer equation in the model atmosphere `atm` with the transitions in `linelist` at the
from `λ_start` to `λ_stop` in steps of `λ_step` (all Å) to get the resultant astrophysical flux at
each wavelength.  An `AbstractRange` can also be probvided directly in place of `λ_start`, `λ_stop`,
and `λ_step`.

Returns a named tuple with keys:
- `flux`: the output spectrum
- `alpha`: the linear absorption coefficient at each wavelenth and atmospheric layer a Matrix of
   size (layers x wavelengths)
- `number_densities`: A dictionary mapping `Species` to vectors of number densities at each
   atmospheric layer
- `wavelengths`: The vacuum wavelenths (in Å) over which the synthesis was performed.  If
  `air_wavelengths=true` this will not be the same as the input wavelenths.

Optional arguments:
- `metallicity`, i.e. [metals/H] is log_10 solar relative
- `abundances` is a `Dict` mapping atomic symbols to ``A(X)`` format abundances, i.e.
   ``A(x) = \\log_{10}(n_X/n_\\mathrm{H}) + 12``, where ``n_X`` is the number density of ``X``.
   These override `metallicity`.
- `vmic` (default: 0) is the microturbulent velocity, ``\\xi``, in km/s.
- `air_wavelengths` (default: `false`): Whether or not the input wavelengths are air wavelenths to
   be converted to vacuum wavelengths by Korg.  The conversion will not be exact, so that the
   wavelenth range can internally be represented by an evenly-spaced range.  If the approximation
   error is greater than `wavelength_conversion_warn_threshold`, an error will be thrown.
- `wavelength_conversion_warn_threshold` (default: 1e-4): see `air_wavelengths`.
- `line_buffer` (default: 10): the farthest (in Å) any line can be from the provided wavelenth range
   before it is discarded.  If the edge of your window is near a strong line, you may have to turn
   this up.
- `cntm_step` (default 1): the distance (in Å) between point at which the continuum opacity is
  calculated.
- `hydrogen_lines` (default: `true`): whether or not to include H lines in the synthesis.
- `mu_grid`: the range of (surface) μ values at which to calculate the surface flux when doing
   transfer in spherical geometry (when `atm` is a `ShellAtmosphere`).
- `line_cutoff_threshold` (default: `1e-3`): the fraction of the continuum absorption coefficient
   at which line profiles are truncated.  This has major performance impacts, since line absorption
   calculations dominate more syntheses.  Turn it down for more precision at the expense of runtime.
   The default value should effect final spectra below the 10^-3 level.
- `ionization_energies`, a `Dict` mapping `Species` to their first three ionization energies,
   defaults to `Korg.ionization_energies`.
- `partition_funcs`, a `Dict` mapping `Species` to partition functions. Defaults to data from
   Barklem & Collet 2016, `Korg.partition_funcs`.
- `equilibrium_constants`, a `Dict` mapping `Species` representing diatomic molecules to their
   molecular equilbrium constants in partial pressure form.  Defaults to data from
   Barklem and Collet 2016, `Korg.equilibrium_constants`.
"""
function synthesize(atm::ModelAtmosphere, linelist, λ_start, λ_stop, λ_step=0.01
                    ; air_wavelengths=false, wavelength_conversion_warn_threshold=1e-4, kwargs...)
    wls = if air_wavelengths
        len = Int(round((λ_stop - λ_start)/λ_step))+1
        vac_start, vac_stop = air_to_vacuum.((λ_start, λ_stop))
        vac_step = (vac_stop - vac_start) / (len-1)
        wls = StepRangeLen(vac_start, vac_step, len)
        max_diff = maximum(abs.(wls .- air_to_vacuum.(λ_start:λ_step:λ_stop)))
        if max_diff > wavelength_conversion_warn_threshold
            throw(ArgumentError("A linear air wavelength range can't be approximated exactly with a"
                                *"linear vacuum wavelength range. This solution differs by up to " *
                                "$max_diff Å.  Adjust wavelength_conversion_warn_threshold if you" *
                                "want to suppress this error."))
        end
        wls
    else
        StepRangeLen(λ_start, λ_step, Int(round((λ_stop - λ_start)/λ_step))+1)
    end
    synthesize(atm, linelist, wls; kwargs...)
end
function synthesize(atm::ModelAtmosphere, linelist, λs::AbstractRange; metallicity::Real=0.0,
                    vmic::Real=1.0, abundances=Dict(), line_buffer::Real=10.0, cntm_step::Real=1.0,
                    hydrogen_lines=true, mu_grid=0:0.05:1, line_cutoff_threshold=1e-3,
                    ionization_energies=ionization_energies,
                    partition_funcs=partition_funcs, equilibrium_constants=equilibrium_constants)
    #work in cm
    λs = λs * 1e-8
    cntm_step *= 1e-8
    line_buffer *= 1e-8
    cntmλs = (λs[1] - line_buffer - cntm_step) : cntm_step : (λs[end] + line_buffer + cntm_step)

    #sort the lines if necessary and check that λs is sorted
    issorted(linelist; by=l->l.wl) || sort!(linelist, by=l->l.wl)
    if step(λs) < 0
        throw(ArgumentError("λs must be in increasing order."))
    end

    #discard lines far from the wavelength range being synthesized
    linelist = filter(l-> λs[1] - line_buffer*1e-8 <= l.wl <= λs[end] + line_buffer*1e-8, linelist)

    abundances = get_absolute_abundances(metallicity, abundances)
    MEQs = molecular_equilibrium_equations(abundances, ionization_energies, partition_funcs,
                                           equilibrium_constants)

    sorted_cntmνs = c_cgs ./ reverse(cntmλs) #frequencies at which to calculate the continuum

    #float-like type general to handle dual numbers
    α_type = typeof(promote(atm.layers[1].temp, length(linelist) > 0 ? linelist[1].wl : 1.0, λs[1],
                            metallicity, vmic, abundances[1])[1])
    #the absorption coefficient, α, for each wavelength and atmospheric layer
    α = Matrix{α_type}(undef, length(atm.layers), length(λs))
    α_cntm = Vector(undef, length(atm.layers))     #vector of continuum-absorption interpolators
    α5 = Vector{α_type}(undef, length(atm.layers)) #each layer's absorption at reference λ (5000 Å)
    n_dicts = Vector(undef, length(atm.layers))    #vector of (species -> number density) Dicts
    for (i, layer) in enumerate(atm.layers)
        n_dicts[i] = molecular_equilibrium(MEQs, layer.temp, layer.number_density,
                                           layer.electron_number_density)

        α_cntm_vals = reverse(total_continuum_absorption(sorted_cntmνs, layer.temp,
                                                         layer.electron_number_density,
                                                         n_dicts[i], partition_funcs))
        α_cntm[i] = LinearInterpolation(cntmλs, α_cntm_vals)
        α[i, :] .= α_cntm[i].(λs)

        α5[i] = total_continuum_absorption([c_cgs/5e-5], layer.temp, layer.electron_number_density,
                                           n_dicts[i], partition_funcs)[1]

        if hydrogen_lines
            α[i, :] .+= hydrogen_line_absorption(λs, layer.temp, layer.electron_number_density,
                                                 n_dicts[i][species"H_I"],
                                                 partition_funcs[species"H_I"],
                                                 hline_stark_profiles, vmic*1e5)
        end
    end

    #put number densities in a dict of vectors, rather than a vector of dicts.
    number_densities = Dict([spec=>[n[spec] for n in n_dicts] for spec in keys(n_dicts[1])])

    #add contribution of line absorption to α
    line_absorption!(α, linelist, λs, [layer.temp for layer in atm.layers],
                     [layer.electron_number_density for layer in atm.layers], number_densities,
                     partition_funcs, vmic*1e5; α_cntm=α_cntm,
                     cutoff_threshold=line_cutoff_threshold)

    source_fn = blackbody.((l->l.temp).(atm.layers), λs')
    flux = radiative_transfer(atm, α, source_fn, α5, mu_grid)

    (flux=flux, alpha=α, number_densities=number_densities, wavelengths=λs.*1e8)
end

"""
    get_absolute_abundances(metallicity, A_X)

Calculate ``n_X/n_\\mathrm{total}`` for each element X given some specified abundances, ``A(X)``.  Use the
metallicity [``X``/H] to calculate those remaining from the solar values (except He).
"""
function get_absolute_abundances(metallicity, A_X::Dict) :: Vector{Number}
    if "H" in keys(A_X)
        throw(ArgumentError("A(H) set, but A(H) = 12 by definition. Adjust \"metallicity\" and "
                           * "\"abundances\" to implicitly set the amount of H"))
    end

    #populate dictionary of absolute abundaces
    abundances = map(0x01:Natoms) do elem
        if elem == 0x01 #hydrogen
            1.0
        elseif atomic_symbols[elem] in keys(A_X)
            10^(A_X[atomic_symbols[elem]] - 12.0)
        else
            #I'm accessing the module global solar_abundances here, but it doesn't make sense to
            #make this an optional argument because this behavior can be completely overridden by
            #specifying all abundances explicitely.
            Δ = (elem == 0x02 #= helium =#) ? 0.0 : metallicity
            10^(solar_abundances[elem] + Δ - 12.0)
        end
    end
    #now normalize so that sum(N_x/N_total) = 1
    abundances ./= sum(abundances)
    abundances
end

"""
    blackbody(T, λ)

The value of the Planck blackbody function for temperature `T` at wavelength `λ`.
"""
function blackbody(T, λ)
    h = hplanck_cgs
    c = c_cgs
    k = kboltz_cgs

    2*h*c^2/λ^5 * 1/(exp(h*c/λ/k/T) - 1)
end


text:
export s_eval_strain_on_gauss_points, s_compute_stiffness_matrix,
s_compute_internal_force_term,
f_eval_strain_on_gauss_points, f_compute_internal_force_term,
s_compute_stiffness_matrix1

@doc raw"""
    s_compute_stiffness_matrix(k::Union{Array{Float64,3}, PyObject})

Computes the small strain stiffness matrix. $k$ is a $n\times 3\times 3$ matrix, where $n$ is the total number of Gauss points.
Returns a SparseTensor.
"""
function s_compute_stiffness_matrix(k::Union{Array{Float64,3}, PyObject}, domain::Domain)
    small_continuum_stiffness_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","small_continuum_stiffness", multiple=true)
    k = convert_to_tensor([k], [Float64]); k = k[1]
    ii, jj, vv = small_continuum_stiffness_(k)
    SparseTensor(ii+1, jj+1, vv, domain.neqs, domain.neqs)
end

@doc raw"""
    s_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject})

Computes the strain on Gauss points in the small strain case. `state` is the full displacement vector.
"""
function s_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject}, domain::Domain)
    small_continuum_strain_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","small_continuum_strain")
    state = convert_to_tensor([state], [Float64]); state = state[1]
    ep = small_continuum_strain_(state)
    set_shape(ep, (getNGauss(domain), 3))
end

@doc raw"""
    f_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject})

Computes the strain on Gauss points in the finite strain case. `state` is the full displacement vector.
"""
function f_eval_strain_on_gauss_points(state::Union{Array{Float64,1}, PyObject}, domain::Domain)
    finit_continuum_strain_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","finit_continuum_strain")
    state = convert_to_tensor([state], [Float64]); state = state[1]
    ep = finit_continuum_strain_(state)
    set_shape(ep, (getNGauss(domain), 3))
end

@doc raw"""
    s_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject})

Computes the internal force
$$\int_\Omega \sigma : \delta \epsilon dx$$
Only active DOFs are considered.
"""
function s_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject}, domain::Domain)
    small_continuum_fint_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","small_continuum_fint")
    stress = convert_to_tensor([stress], [Float64]); stress = stress[1]
    out = small_continuum_fint_(stress)
    set_shape(out, (domain.neqs,))
end

@doc raw"""
    f_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject},
        state::Union{Array{Float64,1}, PyObject},
        domain::Domain)

Computes the internal force for finite strain continuum

$$\int_\Omega \sigma : \delta \epsilon dx$$

Only active DOFs are considered.
"""
function f_compute_internal_force_term(stress::Union{Array{Float64,2}, PyObject},
    state::Union{Array{Float64,1}, PyObject},
    domain::Domain)
    finite_continuum_fint_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","finite_continuum_fint")
    stress,state = convert_to_tensor([stress,state], [Float64,Float64])
    out = finite_continuum_fint_(stress,state)
    set_shape(out, (domain.neqs,))
end

@doc raw"""
    s_compute_stiffness_matrix1(k::Union{PyObject, Array{Float64,3}}, domain::Domain)

Computes the stiffness matrix
```
\int_\Omega (K\nabla u) \cdot \nabla \delta u dx
```
"""
function s_compute_stiffness_matrix1(k::Union{PyObject, Array{Float64,3}}, domain::Domain)
    small_continuum_stiffness1_ = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib","small_continuum_stiffness1", multiple=true)
    k = convert_to_tensor(Any[k], [Float64]); k = k[1]
    ii, jj, vv = small_continuum_stiffness1_(k)
    A = SparseTensor(ii+1, jj+1, vv, domain.neqs, domain.neqs)
end


text:
module TestNamespace

using Omega
using Omega.Inference: swapsinglesite, normalkernel
using UnicodePlots
using Statistics: mean

function test_ssmc_1()
  μ = uniform(0.0, 1.0)
  x = normal(μ, 5.0)
  samples = rand(μ, x ==ₛ 7.0, 10000; alg= SSMH)
  println(histogram([samples...]))
  println(mean(samples))
end

test_ssmc_1()

function test_ssmc_2()
  x = logistic(1.0, 2.0, (2, 3, 4))
  y = sum(x)
  rand(x, y ==ₛ 5.0; alg = SSMH)
end

test_ssmc_2()

function test_sshm_drift()
  x = normal(0.0, 1.0)
  y = normal(1.0, 1.0)
  z = rand(x - y, x ==ₛ y, 10000; alg=SSMH, proposal = (rng, ω) -> swapsinglesite(rng, ω) do x
    normalkernel(rng, x, .1)
  end)
  z = convert(Array{Float64}, z)
  println(histogram([z[300:end]...]))
  println(mean(z[300:end]))
end

test_sshm_drift()

end


text:
function temp_pkg_dir(fn::Function)
  # Used in tests below to setup and teardown a sandboxed package directory
  const tmpdir = ENV["JULIA_PKGDIR"] = abspath(string("tmp.",randstring()))
  @test !isdir(Pkg.dir())
  try
    Pkg.init()
    @test isdir(Pkg.dir())
    Pkg.resolve()

    fn()
  finally
    run(`rm -rf $tmpdir`)
  end
end

# Test adding a removing a package
temp_pkg_dir() do
  @test isempty(Pkg.installed())
  Pkg.add("Example")
  @test [keys(Pkg.installed())...] == ["Example"]
  Pkg.rm("Example")
  @test isempty(Pkg.installed())
end

# testing a package with test dependencies causes them to be installed for the duration of the test
temp_pkg_dir() do
  Pkg.generate("PackageWithTestDependencies", "MIT")
  @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

  isdir(Pkg.dir("PackageWithTestDependencies","test")) || mkdir(Pkg.dir("PackageWithTestDependencies","test"))
  open(Pkg.dir("PackageWithTestDependencies","test","REQUIRE"),"w") do f
    println(f,"Example")
  end

  open(Pkg.dir("PackageWithTestDependencies","test","runtests.jl"),"w") do f
    println(f,"using Base.Test")
    println(f,"@test haskey(Pkg.installed(), \"Example\")")
  end

  Pkg.resolve()
  @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]

  Pkg.test("PackageWithTestDependencies")

  @test [keys(Pkg.installed())...] == ["PackageWithTestDependencies"]
end

# testing a package with no runtests.jl errors
temp_pkg_dir() do
  Pkg.generate("PackageWithNoTests", "MIT")

  if isfile(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
    rm(Pkg.dir("PackageWithNoTests", "test", "runtests.jl"))
  end

  try
    Pkg.test("PackageWithNoTests")
  catch err
    @test err.msg == "PackageWithNoTests did not provide a test/runtests.jl file"
  end
end

# testing a package with failing tests errors
temp_pkg_dir() do
  Pkg.generate("PackageWithFailingTests", "MIT")

  isdir(Pkg.dir("PackageWithFailingTests","test")) || mkdir(Pkg.dir("PackageWithFailingTests","test"))
  open(Pkg.dir("PackageWithFailingTests", "test", "runtests.jl"),"w") do f
    println(f,"using Base.Test")
    println(f,"@test false")
  end

  try
    Pkg.test("PackageWithFailingTests")
  catch err
    @test err.msg == "PackageWithFailingTests had test errors"
  end
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
using RvSpectML

""" prepare_line_list_pass1( linelist_fn, spectra, pipeline; Δv_to_avoid_tellurics, v_center_to_avoid_tellurics )
"""
function prepare_line_list_pass1( linelist_fn::String, all_spectra::AbstractVector{SpecT}, pipeline::PipelinePlan; recalc::Bool = false,
         Δv_to_avoid_tellurics::Real =RvSpectMLBase.max_bc , v_center_to_avoid_tellurics::Real = 0.0 ) where { SpecT <: AbstractSpectra }
   @assert length(linelist_fn) >= 1
   @assert length(all_spectra) >= 1
   if need_to(pipeline,:read_line_list) || recalc
      lambda_range_with_good_data = get_λ_range(all_spectra)
      if verbose println("# Reading line list for CCF: ", linelist_for_ccf_filename, ".")  end
      espresso_filename = joinpath(pkgdir(RvSpectML),"data","masks",linelist_for_ccf_filename)
      espresso_df = RvSpectML.read_linelist_espresso(espresso_filename)
      #inst_module = RvSpectML.get_inst_module(typeof(first(all_spectra).inst))
      #line_list_df = EXPRES.filter_line_list(espresso_df,first(all_spectra).inst)
      line_list_df = filter_line_list(espresso_df,first(all_spectra).inst)
      if eltype(all_spectra) <: AnyEXPRES
         RvSpectML.discard_pixel_mask(all_spectra)
         RvSpectML.discard_excalibur_mask(all_spectra)
      end
      set_cache!(pipeline,:read_line_list,line_list_df)
      dont_need_to!(pipeline,:read_line_list);
    end

   if need_to(pipeline,:clean_line_list_tellurics) || recalc
      if verbose println("# Removing lines with telluric contamination.")  end
      @assert !need_to(pipeline,:read_line_list)
      @assert !need_to(pipeline,:read_spectra)
      if typeof(first(all_spectra).inst) <: AnyEXPRES
         line_list_no_tellurics_df = make_clean_line_list_from_tellurics_expres(line_list_df, all_spectra, Δv_to_avoid_tellurics = Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center_to_avoid_tellurics)
               # RvSpectML.discard_tellurics(all_spectra)  # Keep, since individual line fits use the tellurics info data later
               set_cache!(pipeline,:clean_line_list_tellurics,line_list_no_tellurics_df)
         else
         @warn("Removing lines with telluric contamination currently only works with EXPRES data.")
      end
      dont_need_to!(pipeline,:clean_line_list_tellurics);
    end

    if has_cache(pipeline,:clean_line_list_tellurics) return read_cache(pipeline,:clean_line_list_tellurics)
    elseif has_cache(pipeline,:read_line_list)        return read_cache(pipeline,:read_line_list)
    else   @error("Invalid pipeline state.")          end
end


text:
"""
    dfr = session_frame(ai; odortiming=1, vaccuum=2, soundlick=3, camera=4, sniff=5, lickometer=6, odordirection=7, lickbinary=8, fs=$default_fs)

Construct a trial session DataFrame from an input matrix `ai`. The keyword arguments provide the column index for each signal.
The return value is a `DataFrameRate` which also encodes the sampling frequency (as supplied by `fs`).
"""
function session_frame(ai; odortiming=1, vaccuum=2, soundlick=3, camera=4, sniff=5, lickometer=6, odordirection=7, lickbinary=8, fs=default_fs)
    df = DataFrame()
    if soundlick ∈ axes(ai, 2)
        df.soundlick = view(ai, :, soundlick)
    end
    if lickometer ∈ axes(ai, 2)
        df.lickometer = view(ai, :, lickometer)
    end
    if odortiming ∈ axes(ai, 2)
        df.odortiming = view(ai, :, odortiming)
    end
    if odordirection ∈ axes(ai, 2)
        df.odordirection = view(ai, :, odordirection)
    end
    if sniff ∈ axes(ai, 2)
        df.sniff = view(ai, :, sniff)
    end
    return DataFrameRate(df, fs)
end

"""
    trialdatas(dfr::DataFrameRate; fcarrierlick=60Hz, soundduration=1s, trialoffset=-1s, trialduration=10s, minlickduration=0.005s, maxlickamp=0.1, lohi_sound=(1, 3), lohi_odor=(0.5, 3.0), tprecision=0.1)
    trialdatas(df::DataFrame; fs=$default_fs, fcarrierlick=60Hz, soundduration=1s, trialoffset=-1s, trialduration=10s, minlickduration=0.005s, maxlickamp=0.1, lohi_sound=(1, 3), lohi_odor=(0.5, 3.0), tprecision=0.1)

Collect data on all trials in `df`, which should be a DataFrame with colums "soundlick", "lickometer", "odortiming", and "odordirection".
Each should be an analog signal with timing data.
Supplying `dfr` provides the DataFrame and the sampling rate.

The keyword arguments allow you to set:
- `fs`: the sampling frequency (affects conversion to physical units)
- `fcarrierlick`: the frequency of the carrier signal for lick detection (grounding this signal indicates a lick)
- `soundduration`: the expected duration of the sound pulse. This is used to find trials.
- `tprecision`: fractional "slop" in timing permitted for identifying a sound pulse.
- `trialoffset`: start of the trial relative to the sound pulse onset
- `trialduration`: length of the trial
- `minlickduration`: the minimum length of grounding the lick carrier signal to count as a lick
- `maxlickamp`: maximum allowed deviation of carrier voltage during a lick (as a fraction of lick carrier stddev)
- `lohi_sound`: voltage thresholds for sound lo and hi (a sound pulse starts when this goes hi)
- `lohi_odor`: voltage threshold for odor timing (odor pulse starts when this goes hi)
"""
function trialdatas(df::DataFrame; fs=default_fs, fcarrierlick=60Hz, soundduration=1s, trialoffset=-1s, trialduration=10s, minlickduration=0.005s, maxlickamp=0.1, lohi_sound=(1, 3), lohi_odor=(0.5, 3.0), tprecision=0.1, kwargs...)
    soundlick, lick, odortiming, dir = df.soundlick, df.lickometer, df.odortiming, df.odordirection
    checkcarrier(lick, fs, fcarrierlick; kwargs...)
    lickthresh = std(lick) * maxlickamp
    alltrs = TrialData[]
    i, n = 1, length(soundlick)
    i = advance_to_lo(soundlick, i, n, lohi_sound)
    while i < n
        i = advance_to_hi(soundlick, i, n, lohi_sound)
        i >= n && break
        # We're at a rising edge. When does it fall again?
        j = advance_to_lo(soundlick, i, n, lohi_sound)
        thi = (j-i)/fs   # duration, in time units, of the hi pulse
        if (1-tprecision)*soundduration <= thi <= (1+tprecision)*soundduration
            # We're in a sound marker. Initiate a trial.
            tr = clamp(i + round(Int, float(trialoffset*fs)), 1, n) : clamp(i + round(Int, (trialduration + trialoffset) * float(fs)), 1, n)
            si = (i-first(tr))/fs .. (j-first(tr))/fs
            icue = advance_to_hi(soundlick, j, n, lohi_sound)
            iodor = advance_to_hi(odortiming, tr, lohi_odor)
            jodor = advance_to_lo(odortiming, iodor, last(tr), lohi_odor)
            oi = (iodor-first(tr))/fs .. (jodor-first(tr))/fs
            trs = TrialData(tr, si, (icue - first(tr))/fs, oi, dir[i])
            addlicks!(trs.lickonsets, view(lick, tr); fs, minlickduration, lickthresh)

            push!(alltrs, trs)
            i = last(tr)
        else
            error("soundlick pulse of duration $(uconvert(s, thi)) starting at $i not recognized")
        end
    end
    return alltrs
end
trialdatas(dfr::DataFrameRate; kwargs...) = trialdatas(dfr.df; fs=dfr.fs, kwargs...)

# Detect all licks in the trial. `lick` should be just the voltages recorded during a single trial.
function addlicks!(out, lick; fs, minlickduration, lickthresh)
    k = k1 = firstindex(lick)
    while k < lastindex(lick)
        k += 1
        abs(lick[k]) < lickthresh || continue        # if we're not close to ground, don't look for a lick
        kend = k
        while kend < lastindex(lick)
            kend += 1
            abs(lick[kend]) < lickthresh || break    # keep looking forward until an above-threshold deviation from ground
        end
        if kend - k > fs*minlickduration
            push!(out, (k - k1)/fs)                  # we stayed grounded for long enough, count it as a lick
        end
        k = kend
    end
    return out
end

"""
    is_lick(td::TrialData; deadline=4s)

Return `true` if the animal licked. Only consider licks occurring after the cue time and within `deadline` thereafter.
"""
function is_lick(td::TrialData; deadline=4s)
    tstart = td.cuetime # minimum(td.odorinterval)
    firstidx = searchsortedfirst(td.lickonsets, tstart)
    firstidx > lastindex(td.lickonsets) && return false
    return td.lickonsets[firstidx] - tstart < deadline
end

"""
    collect_by_trialtype(ispos::Function, data, tds::AbstractVector{TrialData})

Categorize `data` by trial type (true positive, true negative, false positive, false negative).
`ispos(td)` should return `true` if this trial should be a positive (lick) response.
"""
function collect_by_trialtype(ispos, data, tds::AbstractVector{TrialData})
    T = typeof(first(data))
    out = Dict("TP" => T[], "TN" => T[], "FP" => T[], "FN" => T[])
    for (val, td) in zip(data, tds)
        isp, isl = ispos(td), is_lick(td)
        key =  isp &  isl ? "TP" :
               isp & !isl ? "FN" :
              !isp &  isl ? "FP" : "TN"
        push!(out[key], val)
    end
    return out
end

"""
    rng = trialrange(ispos::AbstractVector{Bool}, fractionpositive=0.5, pval=0.01)

Compute the largest range over which the probability of `true` values within `ispos` is consistent with `fractionpositive`.
`pval` sets the p-value for determining inconsistency.
"""
function trialrange(ispos, fractionpositive=0.5, pval=0.05)
    # It's not obvious that one can easily compute the distribution of consistent stretches,
    # so do it by simulation. We cache simulation results to save computation time.
    key = (ceil(Int, log2(length(ispos))), Float64(fractionpositive))
    mindist, maxdist = get!(_trialdistribution, key) do
        simulate_trialdist(key...)
    end
    minbounds = getindex_fraction.(mindist, pval/2)
    maxbounds = getindex_fraction.(maxdist, 1-pval/2)
    # Block out the bad stretches
    cispos = cumsum(ispos)
    isok = fill(true, eachindex(ispos))
    n = length(ispos)
    for k = 1:n-1
        mnb, mxb = minbounds[k], maxbounds[k]
        for start = 1:n-k
            all(isok[start:start+k]) || continue
            np = cispos[start+k] - cispos[start]
            if !(mnb <= np <= mxb)
                isok[start:start+k] .= false
            end
        end
    end
    # Find the longest contiguous stretch
    startbest, stopbest = 1, 0
    istart = firstindex(isok)
    while istart < lastindex(isok)
        if !isok[istart]
            istart += 1
            continue
        end
        iend = istart
        while iend < lastindex(isok) && isok[iend+1]
            iend += 1
        end
        if iend - istart > stopbest - startbest
            startbest, stopbest = istart, iend
        end
        istart = iend+1
    end
    return startbest:stopbest
end

const VVI = Vector{Vector{Int}}
const _trialdistribution = Dict{Tuple{Int,Float64},Tuple{VVI,VVI}}()  # cached simulation results

function simulate_trialdist(log2n::Int, fractionpositive::Float64, nsim::Int=10^4)
    n = 2^log2n
    npmindists = [fill(n+1, nsim) for _ = 1:n-1]
    npmaxdists = [fill( -1, nsim) for _ = 1:n-1]
    for isim = 1:nsim
        b = rand(n) .<= fractionpositive
        cb = cumsum(b)
        for k = 1:n-1
            npmin, npmax = typemax(Int), typemin(Int)
            @inbounds @simd for start = 1:n-k
                np = cb[start+k] - cb[start]
                npmin = min(npmin, np)
                npmax = max(npmax, np)
            end
            npmindists[k][isim] = npmin
            npmaxdists[k][isim] = npmax
        end
    end
    foreach(sort!, npmindists)
    foreach(sort!, npmaxdists)
    return npmindists, npmaxdists
end


text:
module Params

using Util

export t_slot, batch_run, verbose, n_runs, t_total, n_channel,
       r_init, B, pkt_size, b_size, pkt_min, base_freq, chan_bw, chan_trans_prob, noise, P_levels

# Simulation parameters:
const batch_run = true
const verbose = true && !batch_run
const debug = false
# number of runs
const n_runs = 10
# time slot
const t_slot = 10e-3 # s
# total simulation time (as time slots)
const t_total = int(floor(600 / t_slot)) # convert seconds to time slots
# number of channels
const n_channel = int8(8)
# radius of initial map
const r_init = 1000

# Buffer parameters

# number of buffer slots
const B = int16(1024) # packets
# Size of a buffer slot (also packet size)
const pkt_size = 1024 # bits
# buffer size in bits
const b_size = B * pkt_size
# package rate for buffer traffic
# defined as a discrete uniform distribution with parameters:
const pkt_min = 0 # pkg / slot, inclusive
const pkt_redundancy = 0.1 # redundancy (error resilience ratio) of packets

# durations
const t_sense = 0.1 * t_slot
const t_sw = 0.05 * t_slot
const t_backoff = Params.t_sense * 3

# powers
const P_tx = 200e-3 # W
const P_levels = [0.5 1 2 4] * P_tx
const P_sense = 0.5*P_tx
const P_sw = 0.5*P_tx
const P_idle = 0.2*P_tx
const P_rec = 40 # W

# channel parameters
const base_freq = 1.8e9 # 1.8 GHz
const chan_bw = 1e6 # 1 MHz
# allowed bitrate
const bitrate = 3.75e6 # 1 Mbps
# state transition probabilities, first state is good state
const chan_trans_prob = [0.7 0.3; 0.7 0.3]'
# channel noises for each channel type (type 1, type 2 etc)
const P_sig = todBm(P_tx/chan_bw) # signal power density in dBm
const noise = P_sig - (32.4 + 20 * (log10(base_freq / 1e9) + log10(r_init))) + [-8 -5; -11 -6] # [-5 -3; -11 -4]

#traffic parameters
const traffic_trans_prob =  [0.7 0.3; 0.9 0.1]'
const traffic_probs = [0.3; 0.7]

# Q learning parameters
module QParams

export beta_overflow, beta_md, beta_loss, epsilon, discount

const beta_overflow = 100
const beta_md = 1 # misdetection punishment coefficient
const beta_loss = 2 # punishment for data loss in channel
const epsilon = 0.05 # exploration probability
const discount = 0.4 # discount factor, gamma

# weights of capability levels
const cap_weights = [1; 1.1]

end

module Indoor

const d_indoor = 2 # m
const walls = 2

end

# saturation time of initial values
const t_saturation = 1000
# capacity of control channel, assuming SNR of channel is 7 dB
const controlcapacity = 1e6 * log2(1 + 10 .^ 0.7)
const twotiers = true
# detection and false alarm probabilities for diferent agent types
const pd = twotiers ? [0.90, 0.90] : [0.88, 0.97]
const pf = twotiers ? [0.10, 0.10] : [0.12, 0.03]
# Location accuracy for different agent types
const eps_accuracy = twotiers ? [10 10] : [50 1]

const energysaving = false

export ParamT, genparams

type ParamT
    goodratio :: (Int, Int)
    beta_idle :: Float64
    sharingperiod :: Int
    buf_levels :: Int
    pkt_max :: Int
    n_stationary_agent :: Int
    iteration :: Int
    prefix :: ASCIIString
    n_agent :: Int8
    n_good_channel :: Int8
    caplevels :: Vector{Int}
    δ :: Float64 # distance factor
    d_svd :: Int # # of dimensions used by SVD compression (SVD compression level)
    trustQ :: Float64 # trust to others' experiences
end

function gencaplevels(n_agent, k, n)
    @assert 0 <= k <= n
    caplevels = ones(Int, n_agent)
    indexes = [1:n_agent] .% n
    indexes[indexes .== 0] = n
    caplevels[indexes .<= k] += 1
    caplevels
end

function genparams()
    params = Array(ParamT,0)
    i = 1
    n_good_channel = int8(3)
    buf_levels = 3
    sharingperiod = 1000

    for d_svd in [11, 12], trustQ in [0.1]
        for δ in [0.6, 1]
            for density in [(6, 8), (10, 10)]
                n_agent = int8(density[1])
                pkt_max = int8(density[2])
                beta_idles = n_agent < 7 ? [1, 2] : [4, 8]
                for beta_idle in beta_idles
                    for goodratio in [(1, 2), (3, 4), (1, 1)]
                        caplevels = gencaplevels(n_agent, goodratio[1], goodratio[2])
                        prefix = @sprintf("%d-%d-%d-%d-%.2f-%d-%d-%.2f-%d-%.2f", n_runs, t_total, n_agent, pkt_max, beta_idle, goodratio[1], goodratio[2], δ, d_svd, trustQ)
                        push!(params, ParamT(goodratio, beta_idle, sharingperiod, buf_levels, pkt_max, div(n_agent, 2), i, prefix, n_agent, n_good_channel, caplevels, δ, d_svd, trustQ))
                        i += 1
                    end
            end
            end
        end
    end
    params
end

end