text:
import SimpleValueGraphs: tuple_of_types

@testset "edge iterator $G" for
    G in (EdgeValGraph, EdgeValOutDiGraph, EdgeValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G((s,d) -> rand_sample(E_VALS), gs, tuple_of_types(E_VALS))

        iter = edges(g)

        E_should_be = (is_directed(g) ? ValDiEdge : ValEdge){V, E_VALS}

        @testset "eltype" begin
            @test eltype(iter) <: E_should_be
        end

        @testset "length" begin
            @test length(iter) == ne(g)
        end

        @testset "correct edges" begin
            @test all(iter) do e
                s, d, vs = src(e), dst(e), vals(e)
                e isa E_should_be &&
                has_edge(g, s, d) &&
                get_edgevals(g, s, d) == vs
            end
        end

        @testset "lexically strictly increasing" begin
            s_prev, d_prev = 0, 0
            @test all(iter) do (e)
                s, d = src(e), dst(e)
                if (s < s_prev ||  (s == s_prev && d <= d_prev))
                    return false
                end
                s_prev, d_prev = s, d
                return true
            end
        end

    end
end

@testset "inneighbors $G" for
    G in (EdgeValGraph, EdgeValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G((s,d) -> rand_sample(E_VALS), gs, tuple_of_types(E_VALS))

        @testset "inneighbors are unique" begin
            @test all(vertices(g)) do u
                allunique(inneighbors(g, u))
            end
        end

        @testset "inneighbors same as simple graph" begin
            @test all(vertices(g)) do u
                issetequal(inneighbors(g, u), inneighbors(gs, u))
            end
        end

    end
end

@testset "outneighbors $G" for
    G in (EdgeValGraph, EdgeValOutDiGraph, EdgeValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G((s,d) -> rand_sample(E_VALS), gs, tuple_of_types(E_VALS))

        @testset "outneighbors are unique" begin
            @test all(vertices(g)) do u
                allunique(outneighbors(g, u))
            end
        end

        @testset "outneighbors same as simple graph" begin
            @test all(vertices(g)) do u
                issetequal(outneighbors(g, u), outneighbors(gs, u))
            end
        end

    end
end


text:
# Structured dataset coordinates can be specified using either a 4D array
# (3, Ni, Nj, Nk), or a tuple (x, y, z).
const Array4 = AbstractArray{T, 4} where T
const Array3Tuple3 = NTuple{3, A} where A <: AbstractArray{T, 3} where T
const StructuredCoords = Union{Array4, Array3Tuple3}

structured_dims(xyz::Array4) = size(xyz)[2:4]
structured_dims(xyz::Array3Tuple3) = size(first(xyz))

function structured_grid(filename::AbstractString,
                         xyz::StructuredCoords;
                         compress=true, append::Bool=true, extent=nothing)
    Ni, Nj, Nk = structured_dims(xyz)
    Npts = Ni * Nj * Nk
    Ncomp = num_components(xyz, Npts)
    Ncls = num_cells_structured(Ni, Nj, Nk)
    ext = extent_attribute(Ni, Nj, Nk, extent)

    if Ncomp != 3  # three components (x, y, z)
        msg = "Coordinate array `xyz` has incorrect dimensions.\n" *
              "Expected dimensions: (3, Ni, Nj, Nk).\n" *
              "Actual dimensions: $(size(xyz))"
        throw(ArgumentError(msg))
    end

    xvtk = XMLDocument()
    vtk = DatasetFile(xvtk, add_extension(filename, ".vts"), "StructuredGrid",
                      Npts, Ncls, compress, append)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # StructuredGrid node
    xGrid = new_child(xroot, vtk.grid_type)
    set_attribute(xGrid, "WholeExtent", ext)

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "Extent", ext)

    # Points node
    xPoints = new_child(xPiece, "Points")

    # DataArray node
    data_to_xml(vtk, xPoints, xyz, "Points", 3)

    return vtk::DatasetFile
end

# 3D variant of vtk_grid with 4D array xyz.
vtk_grid(filename::AbstractString, xyz::AbstractArray{T,4};
         kwargs...) where T =
    structured_grid(filename, xyz; kwargs...)

# 3D variant of vtk_grid with 3D arrays x, y, z.
function vtk_grid(filename::AbstractString, x::AbstractArray{T,3},
                  y::AbstractArray{T,3}, z::AbstractArray{T,3};
                  kwargs...) where T
    if !(size(x) == size(y) == size(z))
        throw(ArgumentError("Size of x, y and z arrays must be the same."))
    end
    structured_grid(filename, (x, y, z); kwargs...)
end

# 2D variant of vtk_grid with 3D array xy
function vtk_grid(filename::AbstractString, xy::AbstractArray{T,3};
                  kwargs...) where T
    Ncomp, Ni, Nj = size(xy)
    if Ncomp != 2
        msg = "Coordinate array `xy` has incorrect dimensions.\n" *
              "Expected dimensions: (2, Ni, Nj).\n" *
              "Actual dimensions: $(size(xy))"
        throw(ArgumentError(msg))
    end
    Nk = 1
    xyz = zeros(T, 3, Ni, Nj, Nk)
    for j = 1:Nj, i = 1:Ni, n = 1:2
        xyz[n, i, j, 1] = xy[n, i, j]
    end
    structured_grid(filename, xyz; kwargs...)
end

# 2D variant of vtk_grid with 2D arrays x, y.
function vtk_grid(filename::AbstractString, x::AbstractArray{T,2},
                  y::AbstractArray{T,2}; kwargs...) where T
    if size(x) != size(y)
        throw(ArgumentError("Size of x and y arrays must be the same."))
    end
    Ni, Nj = size(x)
    Nk = 1
    xyz = zeros(T, 3, Ni, Nj, Nk)
    for j = 1:Nj, i = 1:Ni
        xyz[1, i, j, 1] = x[i, j]
        xyz[2, i, j, 1] = y[i, j]
    end
    structured_grid(filename, xyz; kwargs...)
end

ext: jl


text:
module TestConstructors
    using Base.Test
    using DataTables, DataTables.Index

    #
    # DataTable
    #

    dt = DataTable()
    @test isequal(dt.columns, Any[])
    @test isequal(dt.colindex, Index())

    dt = DataTable(Any[NullableCategoricalVector(zeros(3)),
                       NullableCategoricalVector(ones(3))],
                   Index([:x1, :x2]))
    @test size(dt, 1) == 3
    @test size(dt, 2) == 2

    @test isequal(dt, DataTable(Any[NullableCategoricalVector(zeros(3)),
                                    NullableCategoricalVector(ones(3))]))
    @test isequal(dt, DataTable(x1 = NullableArray([0.0, 0.0, 0.0]),
                                x2 = NullableArray([1.0, 1.0, 1.0])))

    dt2 = convert(DataTable, NullableArray([0.0 1.0;
                                            0.0 1.0;
                                            0.0 1.0]))
    names!(dt2, [:x1, :x2])
    @test isequal(dt[:x1], NullableArray(dt2[:x1]))
    @test isequal(dt[:x2], NullableArray(dt2[:x2]))

    @test isequal(dt, DataTable(x1 = NullableArray([0.0, 0.0, 0.0]),
                                x2 = NullableArray([1.0, 1.0, 1.0])))
    @test isequal(dt, DataTable(x1 = NullableArray([0.0, 0.0, 0.0]),
                                x2 = NullableArray([1.0, 1.0, 1.0]),
                                x3 = NullableArray([2.0, 2.0, 2.0]))[[:x1, :x2]])

    dt = DataTable(Nullable{Int}, 2, 2)
    @test size(dt) == (2, 2)
    @test eltypes(dt) == [Nullable{Int}, Nullable{Int}]

    dt = DataTable([Nullable{Int}, Nullable{Float64}], [:x1, :x2], 2)
    @test size(dt) == (2, 2)
    @test eltypes(dt) == [Nullable{Int}, Nullable{Float64}]

    @test isequal(dt, DataTable([Nullable{Int}, Nullable{Float64}], 2))

    @test_throws BoundsError SubDataTable(DataTable(A=1), 0)
    @test_throws BoundsError SubDataTable(DataTable(A=1), 0)
    @test isequal(SubDataTable(DataTable(A=1), 1), DataTable(A=1))
    @test isequal(SubDataTable(DataTable(A=1:10), 1:4), DataTable(A=1:4))
    @test isequal(view(SubDataTable(DataTable(A=1:10), 1:4), 2), DataTable(A=2))
    @test isequal(view(SubDataTable(DataTable(A=1:10), 1:4), [true, true, false, false]), DataTable(A=1:2))

    @test DataTable(a=1, b=1:2) == DataTable(a=[1,1], b=[1,2])

    @testset "associative" begin
        dt = DataTable(Dict(:A => 1:3, :B => 4:6))
        @test dt == DataTable(A = 1:3, B = 4:6)
        @test eltypes(dt) == [Int, Int]
    end

    @testset "recyclers" begin
        @test DataTable(a = 1:5, b = 1) == DataTable(a = collect(1:5), b = fill(1, 5))
        @test DataTable(a = 1, b = 1:5) == DataTable(a = fill(1, 5), b = collect(1:5))
    end

    @testset "constructor errors" begin
        @test_throws DimensionMismatch DataTable(a=1, b=[])
        @test_throws DimensionMismatch DataTable(Any[collect(1:10)], DataTables.Index([:A, :B]))
        @test_throws DimensionMismatch DataTable(A = rand(2,2))
        @test_throws DimensionMismatch DataTable(A = rand(2,1))
    end

    @testset "column types" begin
        dt = DataTable(A = 1:3, B = 2:4, C = 3:5)
        answer = [Array{Int,1}, Array{Int,1}, Array{Int,1}]
        @test map(typeof, dt.columns) == answer
        dt[:D] = NullableArray([4, 5, Nullable()])
        push!(answer, NullableArray{Int,1})
        @test map(typeof, dt.columns) == answer
        dt[:E] = 'c'
        push!(answer, Array{Char,1})
        @test map(typeof, dt.columns) == answer
    end
end


text:
export findelement

# 3D version
function findelement(xC, yC, zC, location, p, lin)
    ex, ey, ez = size(lin)
    # i
    currentmin = ones(1)
    minind = ones(Int64, 1)
    currentmin[1] = abs.(xC[p[lin[1,1,1]]] .- location[1])
    for i in 2:ex
        current = abs.(xC[p[lin[i,1,1]]] .- location[1])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    i = minind[1]
    # j
    currentmin[1] = abs.(yC[p[lin[1,1,1]]] .- location[2])
    minind[1] = 1
    for i in 2:ey
        current = abs.(yC[p[lin[1,i,1]]] .- location[2])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    j = minind[1]
    # k
    currentmin[1] = abs.(zC[p[lin[1,1,1]]] .- location[3])
    minind[1] = 1
    for i in 2:ez
        current = abs.(zC[p[lin[1,1, i]]] .- location[3])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    k = minind[1]
    return p[lin[i,j,k]]
end

# 2D version
function findelement(xC, yC, location, p, lin)
    ex, ey = size(lin)
    # i
    currentmin = ones(1)
    minind = ones(Int64, 1)
    currentmin[1] = abs.(xC[p[lin[1,1]]] .- location[1])
    for i in 2:ex
        current = abs.(xC[p[lin[i,1]]] .- location[1])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    i = minind[1]
    # j
    currentmin[1] = abs.(yC[p[lin[1,1]]] .- location[2])
    minind[1] = 1
    for i in 2:ey
        current = abs.(yC[p[lin[1,i]]] .- location[2])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    j = minind[1]
    return p[lin[i,j]]
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
# Step 2: Redirect to provider.
function auth(r::HTTP.Request)
    pkey = getquery(r, "provider")
    provider = get(PROVIDERS, pkey, nothing)
    provider === nothing && return html(400, "Requested uknown OAuth provider")

    # If the user has already authenticated, skip.
    state = getcookie(r, "state")
    if !isempty(state) && haskey(USERS, state)
        F = typeof(USERS[state].forge)
        if pkey == "github" && F === GitHubAPI || pkey == "gitlab" && F === GitLabAPI
            # TODO: This does not support custom providers.
            return HTTP.Response(307, ["Location" => ROUTES[:SELECT]])
        end
    end

    state = String(rand('a':'z', 32))
    return HTTP.Response(307, [
        "Set-Cookie" => String(HTTP.Cookie("state", state; path="/"), false),
        "Location" => provider.auth_url * "?" * HTTP.escapeuri(Dict(
            :response_type => "code",
            :client_id => provider.client_id,
            :redirect_uri => callback_url(pkey),
            :scope => provider.scope,
            :state => state,
        )),
    ])
end


text:
Base.@kwdef struct LMDDiffusivity{T} <: AbstractParameters
     CKSL :: T = 0.1   # Surface layer fraction
       Cτ :: T = 0.4   # Von Karman constant

    Cstab :: T = 2.0   # Stable buoyancy flux parameter for wind-driven turbulence
    Cunst :: T = 6.4   # Unstable buoyancy flux parameter for wind-driven turbulence

       Cn :: T = 1.0   # Exponent for effect of stable buoyancy forcing on wind mixing
    Cmτ_U :: T = 0.25  # Exponent for effect of unstable buoyancy forcing on wind mixing of U
    Cmτ_T :: T = 0.5   # Exponent for effect of unstable buoyancy forcing on wind mixing of T
    Cmb_U :: T = 1/3   # Exponent for the effect of wind on convective mixing of U
    Cmb_T :: T = 1/3   # Exponent for effect of wind on convective mixing of T

     Cd_U :: T = 0.5   # Wind mixing regime threshold for momentum
     Cd_T :: T = 2.5   # Wind mixing regime threshold for tracers

     Cb_U :: T = 0.599 # Buoyancy flux parameter for convective turbulence
     Cb_T :: T = 1.36  # Buoyancy flux parameter for convective turbulence
    Cτb_U :: T = (Cτ / Cb_U)^(1/Cmb_U) * (1 + Cunst*Cd_U)^(Cmτ_U/Cmb_U) - Cd_U  # Wind stress parameter for convective turbulence
    Cτb_T :: T = (Cτ / Cb_T)^(1/Cmb_T) * (1 + Cunst*Cd_T)^(Cmτ_T/Cmb_T) - Cd_T  # Wind stress parameter for convective turbulence

      KU₀ :: T = 1e-6 # Interior viscosity for velocity
      KT₀ :: T = 1e-7 # Interior diffusivity for temperature
      KS₀ :: T = 1e-9 # Interior diffusivity for salinity
end

## ** The K-Profile-Parameterization **
K_KPP(h, 𝒲, d::T, p) where T = 0<d<1 ? max(zero(T), h * 𝒲 * shape(d, p)) : -zero(T)

𝒲_Holtslag(Cτ, Cτb, ωτ, ωb, d) = Cτ * (ωτ^3 + Cτb * d * ωb^3)^(1/3)
𝒲_Holtslag(m, i) = 𝒲_Holtslag(m.diffusivity.Cτ, m.diffusivity.Cτb, KPP.ωτ(m), KPP.ωb(m), KPP.d(m, i))

𝒲_LMD_unstable_U(m, i) = KPP.𝒲_unstable(
    m.diffusivity.CKSL, m.diffusivity.Cd_U,
    m.diffusivity.Cτ, m.diffusivity.Cunst,
    m.diffusivity.Cb_U, m.diffusivity.Cτb_U,
    m.diffusivity.Cmτ_U, m.diffusivity.Cmb_U,
    ωτ(m), ωb(m), d(m, i)
    )

𝒲_LMD_unstable_T(m, i) = KPP.𝒲_unstable(
    m.diffusivity.CKSL, m.diffusivity.Cd_T,
    m.diffusivity.Cτ, m.diffusivity.Cunst,
    m.diffusivity.Cb_T, m.diffusivity.Cτb_T,
    m.diffusivity.Cmτ_T, m.diffusivity.Cmb_T,
    ωτ(m), ωb(m), d(m, i)
    )

𝒲_LMD_stable(m, i) = KPP.𝒲_stable(
    m.diffusivity.Cτ, m.diffusivity.Cstab, m.diffusivity.Cn,
    ωτ(m), ωb(m), d(m, i)
    )

"Return the vertical velocity scale for momentum at face point i"
function 𝒲_LMD_U(m, i)
    if !isforced(m)
        return 0
    elseif isunstable(m)
        return 𝒲_LMD_unstable_U(m, i)
    else
        return 𝒲_LMD_stable(m, i)
    end
end

"Return the vertical velocity scale for tracers at face point i."
function 𝒲_LMD_T(m, i)
    if !isforced(m)
        return 0
    elseif isunstable(m)
        return 𝒲_LMD_unstable_T(m, i)
    else
        return 𝒲_LMD_stable(m, i)
    end
end

const 𝒲_LMD_V = 𝒲_LMD_U
const 𝒲_LMD_S = 𝒲_LMD_T

Base.@kwdef struct HoltslagDiffusivity{T} <: AbstractParameters
     Cτ :: T = 0.4
    Cτb :: T = 15.6
    KU₀ :: T = 1e-6 # Interior viscosity for velocity
    KT₀ :: T = 1e-7 # Interior diffusivity for temperature
    KS₀ :: T = 1e-9 # Interior diffusivity for salinity
end


text:
using DynamicCorporateModels
using Test

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


text:
module NaturalES
using Random
using LinearAlgebra
abstract type OptMethod end

include("utils.jl")
include("separable_nes.jl")
include("exponential_nes.jl")

"""
optimize(f, μ, σ, [method=sNES;options...])

minimizes the function `f` according to:

    `f` : function to optimize
        μ::Vector -> cost::Real
    `μ` : initial condition
        μ::Vector
    `σ` : initial uncertainty on μ
        σ::{Real | Vector | Matrix}
    `method` : xNES or sNES
        xNES = exponential evolution strategies, expensive but powerful on non separable objective
        sNES = separable evolution strategies, lightweight very powerful for separable or very high dimensional objectives
    `options` :
             ημ = learning rate for μ,
             ησ = learning rate for uncertainties,
           atol = tolerance on uncertainties (default 1e-8),
        samples = number of samples used to build Natural Gradient approximation,
        iterations = upper limit on the number of iterations, default 10^4)

"""
function optimize(f,μ,σ,method::Type{<:OptMethod}=sNES;options...)
    pickmethod(method)(f,μ,σ;options...)
end

export sNES,xNES,optimize
end # module


text:
# Exact representation of algebraic numbers
# (Numbers that are roots of polynomials with integer coefficients)
# And also arithmetic on algebraic numbers,
# including +, -, *, /, and radicals.

using Nemo
import PolynomialRoots
import PolynomialRoots:roots

import Base.zero,Base.one
import Base.+,Base.-,Base.*,Base./,Base.inv
import Base.abs,Base.conj
import Base.real,Base.imag
import Base.==

# see: http://nemocas.org/nemo-0.4.pdf

# An algebraic number,
# consisting of the minimal polynomial of the number,
# an arbitrary-precision approximation of the number,
# and prec which specifies the minimal distance between roots of p
# TODO: apprx has to be complex.
struct AlgebraicNumber{T<:Integer,F<:AbstractFloat} <: Number
	coeff::Vector{T}
	apprx::Complex{F}
	prec::F
end

# algebraic number from just poly and approximation.
# computes precision and simplifies as well.
function AlgebraicNumber(coeff::Vector{T}, apprx::Complex{F}) where {T<:Integer,F<:AbstractFloat}
	an = AlgebraicNumber{T,F}(coeff, apprx, calc_precision(coeff, apprx))
	return simplify(an)
end

# AlgebraicNumber from any integer type
AlgebraicNumber(x::T) where {T<:Integer} =
    AlgebraicNumber(BigInt[-x,one(T)], Complex{BigFloat}(x))

# AlgebraicNumber from rationails
AlgebraicNumber(x::Rational) =
    AlgebraicNumber(BigInt[-numerator(x), denominator(x)], Complex{BigFloat}(x))

function poly_from_coeff(a)
	R,x=PolynomialRing(Nemo.FlintZZ,"x")
	sum([a[i]*x^(i-1) for i=1:length(a)])
end

import Base.show
# TODO: only show up to precision
function show(io::IO, an::AlgebraicNumber)
	print(io,"≈")
	#ndigits = max(10, round(Int,ceil(convert(Float64,log(an.prec)/log(10)))))
	show(io,convert(Complex{Float64},an.apprx))
	#print(io,"...")
end

#get_coeffs(p::Nemo.fmpz_poly) = pointer_to_array(convert(Ptr{Int64}, p.coeffs), (p.length,))
get_coeffs(p::Nemo.fmpz_poly) = [BigInt(Nemo.coeff(p,i)) for i=0:Nemo.degree(p)]
prec_roots(a::Vector{T}) where {T<:Integer} = PolynomialRoots.roots(convert(Array{BigFloat},a))
# TODO: make sure roots returns distinct roots

# Given an algebraic number, find minimum precision required
# to specify it among roots of an.p
# TODO: handle case of repeated roots precisely
function calc_precision(coeff::Vector{T}, apprx::Complex{F}) where {T<:Integer,F<:AbstractFloat}
	# compute smallest distance between all pairs of elements in x
	function min_pairwise_dist(x)
		biginf = convert(BigFloat,Inf)
		n = length(x)
		if n<=1
			return Inf
		else
			pdists = [i < j ? abs(x[i]-x[j]) : biginf for i=1:n,j=1:n]
			return minimum(pdists)
		end
	end

	# first, find all roots of p
	rts = prec_roots(coeff)

	# first, trivial case
	if length(rts)==1
		return convert(BigFloat, Inf)
	end

	# find minimum pairwise distance between roots;
	# multiply by 0.5 safety factor
	return 0.5*min_pairwise_dist(rts)
end

# simplify an algebraic number by reducing p to the minimal polynomial.
# This assumes that calc_precision! has already been called.
function simplify(an::AlgebraicNumber)
	# for all factors of an.p, find the one that matches our roots

	# If linear polynomial, then already irreducible.
	if length(an.coeff)<=2
		return an
	end

	# Otherwise, factor out.
	R, x = PolynomialRing(Nemo.FlintZZ, "x")
	p = R(map(Nemo.FlintZZ, an.coeff))
	fctr_dict = Nemo.factor(p)
	#fctrs = keys(fctr_dict)
	fctrs = [p for (p,e) in fctr_dict]

	# first, trivial case
	if length(fctrs)==1
		# irreducible case
		if first(values(fctr_dict))==1
			return AlgebraicNumber(get_coeffs(first(fctrs)),an.apprx,an.prec)
		end
		# reducible case
		coeffs1 = get_coeffs(first(fctrs))
		apprx1  = an.apprx
		return AlgebraicNumber(coeffs1,apprx1,calc_precision(coeffs1,apprx1))
	end
	# case where more than one factor exists
	mindists = [minimum(abs.(an.apprx .- prec_roots(get_coeffs(fctr)))) for fctr in fctrs]
	(newprec, i) = findmin(mindists)
	fctr = collect(fctrs)[i]
	an = AlgebraicNumber(get_coeffs(fctr),an.apprx,newprec)
	return an
end

function ==(an1::AlgebraicNumber,an2::AlgebraicNumber)
	cf1 = an1.coeff
	cf2 = an2.coeff
	(cf1/cf1[end])==(cf2/cf2[end]) || return false
	prec1 = calc_precision(an1.coeff, an1.apprx)
	prec2 = calc_precision(an2.coeff, an2.apprx)
	return abs(an1.apprx-an2.apprx)<min(prec1,prec2)
end

inv(an::AlgebraicNumber) = AlgebraicNumber(reverse(an.coeff), inv(an.apprx))

# interleave each elemnet of a with n zeros
interleave(a,n) =  vec(vcat(a',zeros(Int64,n,length(a))))
function root(an::AlgebraicNumber,n::Int64)
	if n==0
		throw(ArgumentError("n must be nonzero"))
	end
	if n==1
		return an
	end
	if n < 0
		an = inv(an)
		n = -n
	end
	# TODO: quickly calculate precision
	return AlgebraicNumber(interleave(an.coeff, n-1), an.apprx^(1/n))
end

import Base.sqrt
import Base.cbrt
sqrt(an::AlgebraicNumber) = root(an,2)
cbrt(an::AlgebraicNumber) = root(an,3)

# TODO: special, more efficient cases for ^2 and ^3
function pow2(an::AlgebraicNumber)
	cfs = an.coeff
	# first check if it is already in the form of a square root.
	if all(cfs[2:2:end] .== 0)
		pp_cfs = cfs
	else
		cfs2 = [iseven(i) ? -cfs[i] : cfs[i] for i=1:length(cfs)]
		pp = poly_from_coeff(cfs)*poly_from_coeff(cfs2)
		pp_cfs = get_coeffs(pp)
	end
	p2 = pp_cfs[1:2:end]
	return AlgebraicNumber(p2, an.apprx*an.apprx)
end

# partially simplify a polynomial b
# eliminating repeated factors
reduce_repeated_factors(p::Nemo.fmpz_poly) = prod(keys(Nemo.factor(p)))

# multiplication
function *(an1::AlgebraicNumber,an2::AlgebraicNumber)
	if an1==0 || an2==0
		# TODO: don't handle this explicitly
		return zero(AlgebraicNumber)
	end
	# check if p==q, if then use a more optimized and correct routine
	#if an1.coeff == an2.coeff
	#	return
	#end
	p = composed_product(an1.coeff, an2.coeff)
	return AlgebraicNumber(p, an1.apprx * an2.apprx)
end

function +(an1::AlgebraicNumber,an2::AlgebraicNumber)
	p = composed_sum(an1.coeff, an2.coeff)
	return AlgebraicNumber(p, an1.apprx + an2.apprx)
end

function -(an1::AlgebraicNumber)
	cfs = copy(an1.coeff)
	for i=1:2:length(cfs)
		cfs[i]=-cfs[i]
	end
	return AlgebraicNumber(cfs, -an1.apprx, an1.prec)
end

-(an1::AlgebraicNumber,an2::AlgebraicNumber) = an1+(-an2)
/(an1::AlgebraicNumber,an2::AlgebraicNumber) = an1*(inv(an2))

# the complex conjugate of an algebraic number has the same minimal polynomial
conj(an::AlgebraicNumber) = AlgebraicNumber(an.coeff,conj(an.apprx),an.prec)
abs(an::AlgebraicNumber) = sqrt(an*conj(an))

zero(::Type{AlgebraicNumber}) = AlgebraicNumber(BigInt[0, 1],Complex{BigFloat}(0.0),BigFloat(1.0))
one(::Type{AlgebraicNumber})  = AlgebraicNumber(BigInt[-1,1],Complex{BigFloat}(1.0),BigFloat(1.0))

real(an::AlgebraicNumber) = (an+conj(an))*AlgebraicNumber(BigInt[1,-2], BigFloat(0.5)+0im,BigFloat(0.5))
imag(an::AlgebraicNumber) = (an-conj(an))*AlgebraicNumber(BigInt[1,0,4],BigFloat(-0.5)*im,BigFloat(0.5))

# take roots of a polynomial,
# and return them as algebraic numbers
function alg_roots(coeff::Vector{Integer})
	#TODO
end

confirm_algnumber(b) = sum(b.coeff .* [b.apprx^(i-1) for i=1:length(b.coeff)])

# compute exp(pi*i*a),
# which is algebraic if a is rational.
function exp_alg(a::Rational)
	# first, obtain polynomial
	p = interleave(BigInt[-1,1], 2*denominator(a)-1)
	# now, select root.
	apprx = exp(im*BigFloat(pi)*a)
	# Finally, return minimal polynomial w.r.t. that root
	return AlgebraicNumber(p,apprx)
end

cos_alg(a::Rational) = real(exp_alg(a))
sin_alg(a::Rational) = imag(exp_alg(a))


text:
using GpABC, OrdinaryDiffEq, Distances, Distributions, Dates, Logging

# seed!(4)

#
# ABC settings
#

n_particles = 1000
progress_every = 1000

#
# Emulation settings
#
n_design_points = 100

#
# True parameters
#
true_params =  [2.0, 2.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.),
            Uniform(75., 125.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)]
param_indices = [1,2,3]
n_var_params = length(param_indices)

#
# ODE solver settings
#
Tspan = (0.0, 10.0)
x0 = [0.0, 0.0, 0.0]
solver = RK4()
# saveat = 0.1
# threshold_schedule = [3.0, 2.0, 1.0, 0.5, 0.2]
saveat = 0.001
threshold_schedule = [3.0, 2.0, 1.0, 0.5, 0.2] .* 5

function ODE_3GeneReg(dx, x, pars, t)
  dx[1] = pars[1]/(1+pars[7]*x[3]) - pars[4]*x[1]
  dx[2] = pars[2]*pars[8]*x[1]/(1+pars[8]*x[1]) - pars[5]*x[2]
  dx[3] = pars[3]*pars[9]*x[1]*pars[10]*x[2]/((1+pars[9]*x[1])*(1+pars[10]*x[2])) - pars[6]*x[3]
end

#
# Returns the solution to the toy model as solved by OrdinaryDiffEq
#
function GeneReg(params::AbstractArray{Float64,1},
    Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1},
    solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

  prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params)
  Obs = solve(prob, solver, saveat=saveat)

  return Array{Float64, 2}(Obs)
end

reference_data = GeneReg(true_params, Tspan, x0, solver, saveat)
# reference_data += randn(size(reference_data)) * 0.1

function simulator_function(var_param_idx::Int, var_param::Real)
    params = copy(true_params)
    params[var_param_idx] = var_param
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulator_function(var_param_indices::AbstractArray{Int, 1}, var_params::AbstractArray{T, 1}) where {T<:Real}
    params = copy(true_params)
    params[var_param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulator_function(var_params::AbstractArray{T, 1}) where {T<:Real}
    simulator_function(param_indices, var_params)
    params = copy(true_params)
    params[param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulate_distance(var_params)
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_params)))
end

function simulate_distance(var_param_idx::Int, var_param::Real)
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_param_idx, var_param)))
end

function simulate_distance(var_param_indices::AbstractArray{Int, 1}, var_params::AbstractArray{T, 1}) where {T<:Real}
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_param_indices, var_params)))
end

function run_emulation()
    @info("EMULATION")
    start_ts = now()
    emu_out = EmulatedABCSMC(reference_data,
        simulator_function,
        priors[param_indices],
        threshold_schedule,
        n_particles,
        n_design_points;
        batch_size=1000,
        emulator_retraining = PreviousPopulationThresholdRetraining(n_design_points, min(100, n_design_points/2), 10),
        emulated_particle_selection = MeanVarEmulatedParticleSelection(),
        # write_progress=true
        )
    elapsed_msg = canonicalize(Dates.CompoundPeriod(now() - start_ts))
    @info "Emulation took $(elapsed_msg)"
    emu_out
end

function run_simulation()
    println("SIMULATION")
    start_ts = now()
    sim_out = SimulatedABCSMC(reference_data,
        simulator_function,
        priors[param_indices],
        threshold_schedule,
        n_particles,
        # write_progress=true
        )
    elapsed_msg = canonicalize(Dates.CompoundPeriod(now() - start_ts))
    @info "Simulation took $(elapsed_msg)"
    sim_out
end

emu_out = run_emulation()
sim_out = run_simulation()
tparams = true_params[param_indices]

#=
import Pkg; Pkg.activate("plots", shared=true)
using Plots
pyplot()
plot(emu_out, sim_out)
plot(emu_out)
plot(sim_out)
=#


text:
module OceanColorData

include("Conversions.jl")
include("Classifiers.jl")

export RemotelySensedReflectance, RrsToChla
export FuzzyClassification, Moore2009, Jackson2017

end # module


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

include("selection/tag_expressions_test.jl")
include("gherkin/experimental/runtests.jl")

include("gherkin/feature_test.jl")
include("gherkin/scenario_test.jl")
include("gherkin/scenario_outline_test.jl")
include("gherkin/scenario_descriptions_test.jl")
include("gherkin/tag_test.jl")
include("gherkin/data_tables_test.jl")

include("selection/tag_selection_test.jl")

include("executor_test.jl")
include("executor_presentation_test.jl")
include("executor_feature_test.jl")
include("executor_datatables_test.jl")
include("step_def_test.jl")
include("result_accumulator_test.jl")
include("asserts_test.jl")
include("outlines_test.jl")
include("exec_env_test.jl")
include("variables_test.jl")
include("suggestion_test.jl")


text:
module MolecularConformation

export  preprocessing,NMRInfo,NMRType,conformation,ConformationSetup,ConformationOutput,
AtomType,MoleculeType,classicBP, classicBP_closure,
quaternionBP,≈,generate_virtual_path, bondangle,torsionmatrix,badtorsionangle,
pruningtest,LDE,build_distance_matrix,outputfilter,writefile,convert_to_dataframe

# loading basic packages
using LinearAlgebra,DelimitedFiles,SparseArrays,DataFrames
#using Quaternions
import Base.show

include("utils.jl")
include("solvers.jl")
include("show.jl")

# main function to run the solver and related problem
"""
```
conformation(NMRdata,cs::ConformationSetup)

```
This is the main function in order to get the conformation of a molecule. To run this function, we need to setup a NMR file and then define some options using the ConformationSetup type.

## Example

```julia-repl
julia> options = ConformationSetup(0.001,classicBP,true)

julia> conformation(NMRdata,options)
```
as return a ConformationOutput type is provided.

There are others parameters to setup, for example, need to complete.
"""
function conformation(NMRdata::NMRType,
		     cs::ConformationSetup)

	solutions = cs.solver(NMRdata,cs.precision,cs.virtual_precision,cs.allsolutions)
	s = ConformationOutput(cs.solver,solutions[1],solutions[2],solutions[3])
	if cs.evalLDE == true
		map(i->MolecularConformation.LDE(s.molecules[i],NMRdata),[1:1:s.number;])
	end
	return s
end

end


text:
export Interval

## Standard interval
# T Must be a Vector space
immutable Interval{T} <: IntervalDomain{T}
	a::T
	b::T
	Interval()=new(-one(T),one(T))
	Interval(a,b)=new(a,b)
end

Interval()=Interval{Float64}()
Interval(a::Int,b::Int) = Interval(Float64(a),Float64(b)) #convenience method
Interval(a,b) = Interval{promote_type(typeof(a),typeof(b))}(a,b)
Interval(a::Tuple,b::Tuple)=Interval(Vec(a...),Vec(b...))

function Interval{T<:Number}(d::AbstractVector{T})
    @assert length(d)==2
    @assert isfinite(d[1]) && isfinite(d[2])
    Interval(d...)
end

Base.convert{T<:Number}(::Type{Interval{T}}, d::Interval) = Interval{T}(d.a,d.b)

AnyInterval{T}(::Type{T}) = Interval{T}(NaN,NaN)
AnyInterval() = AnyInterval(Float64)
isambiguous(d::Interval) = all(isnan(d.a)) && all(isnan(d.b))
Base.convert{T<:Number}(::Type{Interval{T}},::AnyDomain) = AnyInterval(T)
Base.convert(::Type{Interval},::AnyDomain) = AnyInterval()

## Information

Base.first(d::Interval)=d.a
Base.last(d::Interval)=d.b
Base.isempty(d::Interval)=isapprox(d.a,d.b;atol=200eps(eltype(d)))

Base.issubset(a::Interval,b::Interval)=first(a)∈b && last(a)∈b

## Map interval
# The first definition  is the more general

mobius(S::Space,x...) = mobius(domain(S),x...)

tocanonical{T}(d::Interval{T},x::T) = 2norm(x-d.a)/arclength(d)-1
tocanonical{T<:Complex}(d::Interval{T},x::Number) = 2norm(x-d.a)/arclength(d)-1
tocanonical{T}(d::Interval{T},x::AbstractVector{T}) = map(x->tocanonical(d,x),x)
mobius(d::Interval,x) = (d.a + d.b - 2x)/(d.a - d.b)
tocanonical{T<:Real}(d::Interval{T},x::AbstractVector{T}) = mobius(d,x)
tocanonical{T<:Real}(d::Interval{T},x) = mobius(d,x)
tocanonicalD{T<:Real}(d::Interval{T},x) = 2/(d.b- d.a)
fromcanonical{T<:Number}(d::Interval{T},v::AbstractArray)=promote_type(T,eltype(v))[fromcanonical(d,vk) for vk in v]
fromcanonical{T<:Number}(d::Interval{T},x) = (d.a + d.b)/2 + (d.b - d.a)x/2
fromcanonical{T<:Vec}(d::Interval{T},x::Number) = (d.a + d.b)/2 + (d.b - d.a)x/2
fromcanonicalD(d::Interval,x) = (d.b- d.a) / 2

arclength(d::Interval) = norm(d.b - d.a)
Base.angle(d::Interval)=angle(d.b-d.a)
complexlength(d::Interval)=d.b-d.a

==(d::Interval,m::Interval) = d.a == m.a && d.b == m.b
function Base.isapprox(d::Interval,m::Interval)
    tol=10E-12
    norm(d.a-m.a)<tol&&norm(d.b-m.b)<tol
end

## algebra

for op in (:*,:+,:-,:.*,:.+,:.-,:.^)
    @eval begin
        $op(c::Number,d::Interval)=Interval($op(c,d.a),$op(c,d.b))
        $op(d::Interval,c::Number)=Interval($op(d.a,c),$op(d.b,c))
    end
end

for op in (:/,:./)
    @eval $op(d::Interval,c::Number)=Interval($op(d.a,c),$op(d.b,c))
end

Base.sqrt(d::Interval)=Interval(sqrt(d.a),sqrt(d.b))

+(d1::Interval,d2::Interval)=Interval(d1.a+d2.a,d1.b+d2.b)

## intersect/union

Base.reverse(d::Interval)=Interval(d.b,d.a)

function Base.intersect{T<:Real,V<:Real}(a::Interval{T},b::Interval{V})
    if first(a) > last(a)
        intersect(reverse(a),b)
    elseif first(b) > last(b)
        intersect(a,reverse(b))
    elseif first(a) > first(b)
        intersect(b,a)
    elseif last(a) <= first(b)
        EmptyDomain()
    elseif last(a)>=last(b)
        b
    elseif isapprox(first(b),last(a);atol=100eps(promote_type(T,V))/max(arclength(a),arclength(b)))
        EmptyDomain()
    else
        Interval(first(b),last(a))
    end
end

function Base.setdiff{T<:Real,V<:Real}(a::Interval{T},b::Interval{V})
    # ensure a/b are well-ordered
    if first(a) > last(a)
        intersect(reverse(a),b)
    elseif first(b) > last(b)
        intersect(a,reverse(b))
    elseif first(a)< first(b)
        if last(a) <= first(b)
            a
        else # first(a) ≤ first(b) ≤last(a)
            #TODO: setdiff in the middle
            if last(a) <= last(b)
            	Interval(first(a),first(b))
			else  # first(a) ≤ first(b) ≤ last(b) ≤last(a)
				Interval(first(a),first(b))∪Interval(last(b),last(a))
			end
        end
    else #first(a)>= first(b)
        if first(a)>=last(b)
            a
        elseif last(a) <= last(b)
            EmptyDomain()
        else #first(b) < first(a) < last(b) < last(a)
            Interval(last(b),last(a))
        end
    end
end

## sort
Base.isless{T1<:Real,T2<:Real}(d1::Interval{T1},d2::Interval{T2})=d1≤first(d2)&&d1≤last(d2)
Base.isless{T<:Real}(d1::Interval{T},x::Real)=first(d1)≤x && last(d1)≤x
Base.isless{T<:Real}(x::Real,d1::Interval{T})=x≤first(d1) && x≤last(d1)


text:
part1(input) = optimize_equipment(input, part = 1)
part2(input) = optimize_equipment(input, part = 2)

function optimize_equipment(input; part)
    boss_stats = parse.(Int, last.(split.(eachline(input), ": ")))
    hero_hp = 100
    weapons = [(  8, 4, 0),
               ( 10, 5, 0),
               ( 25, 6, 0),
               ( 40, 7, 0),
               ( 74, 8, 0)]
    armors = [(  0, 0, 0),
              ( 13, 0, 1),
              ( 31, 0, 2),
              ( 53, 0, 3),
              ( 75, 0, 4),
              (102, 0, 5)]
    rings = [(  0, 0, 0),
             ( 25, 1, 0),
             ( 50, 2, 0),
             (100, 3, 0),
             ( 20, 0, 1),
             ( 40, 0, 2),
             ( 80, 0, 3)]

    best_cost = part == 1 ? typemax(Int) : 0
    for weapon in weapons, armor in armors, ring1 in rings, ring2 in rings
        0 < ring2[1] <= ring1[1] && continue
        ring2[1] < ring1[1] && continue
        cost, hero_damage, hero_armor = weapon .+ armor .+ ring1 .+ ring2
        part == 1 && cost >= best_cost && continue
        part == 2 && cost <= best_cost && continue
        if fight(hero_hp, hero_damage, hero_armor, boss_stats...) ⊻ (part == 2)
            best_cost = cost
        end
    end
    return best_cost
end

function fight(hp1, damage1, armor1, hp2, damage2, armor2)
    while true
        hp2 -= max(1, damage1 - armor2)
        hp2 <= 0 && return true
        hp1 -= max(1, damage2 - armor1)
        hp1 <= 0 && return false
    end
end


text:
module AllenNeuropixels
using PyCall
using Preferences
using Conda
using DataFrames
using Dates
using TimeZones
using CSV
using Requires

const allensdk = PyNULL()
const brain_observatory = PyNULL()
const ecephys = PyNULL()
const ecephys_project_cache = PyNULL()
const brain_observatory_cache = PyNULL()
const stimulus_info = PyNULL()
const mouse_connectivity_cache = PyNULL()
const ontologies_api = PyNULL()
const reference_space_cache = PyNULL()
const reference_space =PyNULL()
export allensdk, brain_observatory, ecephys, ecephys_project_cache, mouse_connectivity_cache, ontologies_api, reference_space_cache, reference_space

function __init__()
    Conda.pip_interop(true)
    Conda.update() # You might need to delete the .julia/conda folder and rebuild python; allensdk has some tricky compatibility requirements.
    Conda.pip("install", "allensdk")
    copy!(allensdk, pyimport("allensdk"))
    copy!(brain_observatory, pyimport("allensdk.brain_observatory"))
    copy!(stimulus_info, pyimport("allensdk.brain_observatory.stimulus_info"))
    copy!(ecephys, pyimport("allensdk.brain_observatory.ecephys"))
    copy!(ecephys_project_cache, pyimport("allensdk.brain_observatory.ecephys.ecephys_project_cache"))
    copy!(brain_observatory_cache, pyimport("allensdk.core.brain_observatory_cache"))
    copy!(mouse_connectivity_cache, pyimport("allensdk.core.mouse_connectivity_cache"))
    copy!(ontologies_api, pyimport("allensdk.api.queries.ontologies_api"))
    copy!(reference_space_cache, pyimport("allensdk.core.reference_space_cache"))
    copy!(reference_space, pyimport("allensdk.core.reference_space"))

    ecephys_project_cache.EcephysProjectCache.from_warehouse(manifest=ecephysmanifest)

    @require WGLMakie="276b4fcb-3e11-5398-bf8b-a0c2d153d008" @eval using .Plots
    @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" @eval using .Plots
end

function setdatadir(datadir::String)
    @set_preferences!("datadir" => datadir)
    @info("New default datadir set; restart your Julia session for this change to take effect")
end
const datadir = replace(@load_preference("datadir", joinpath(pkgdir(AllenNeuropixels), "data/")), "\\"=>"/")
const ecephysmanifest = replace(joinpath(datadir, "Ecephys", "manifest.json"), "\\"=>"/")
const brainobservatorymanifest = replace(joinpath(datadir, "BrainObservatory", "manifest.json"), "\\"=>"/")
const mouseconnectivitymanifest = replace(joinpath(datadir, "MouseConnectivity", "manifest.json"), "\\"=>"/")
const referencespacemanifest = replace(joinpath(datadir, "ReferenceSpace", "manifest.json"), "\\"=>"/")

function loaddataframe(file, dir=datadir)
    CSV.File(abspath(dir, file)) |> DataFrame
end
export convertdataframe

include("./EcephysCache.jl")
include("./BrainObservatory.jl")
include("./LFP.jl")
include("./MouseConnectivityCache.jl")
include("./Plots/Plots.jl")
include("./Ontologies.jl")
include("./ReferenceSpace.jl")

end


text:
module EMLbatch
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("batch/batch_basic.jl")
	include("batch/batch_isothermic.jl")
end


text:
export BreadthFirstPlanner

"Uninformed breadth-first search planner."
@kwdef mutable struct BreadthFirstPlanner <: Planner
    max_nodes::Int = typemax(Int)
    save_search::Bool = false # Flag to save search info
end

function solve(planner::BreadthFirstPlanner,
               domain::Domain, state::State, spec::Specification)
    @unpack max_nodes, save_search = planner
    # Initialize backpointers and queue
    node_id = hash(state)
    search_tree = Dict{UInt,PathNode}(node_id => PathNode(node_id, state, 0))
    queue = [node_id]
    # Run the search
    status, node_id = search!(planner, domain, spec, search_tree, queue)
    # Reconstruct plan and return solution
    if status != :failure
        plan, traj = reconstruct(node_id, search_tree)
        if save_search
            return PathSearchSolution(status, plan, traj, search_tree, queue)
        else
            return PathSearchSolution(status, plan, traj)
        end
    elseif save_search
        return PathSearchSolution(status, [], [], search_tree, queue)
    else
        return NullSolution()
    end
end

function search!(planner::BreadthFirstPlanner,
                 domain::Domain, spec::Specification,
                 search_tree::Dict{UInt,PathNode}, queue::Vector{UInt})
    count = 1
    while length(queue) > 0
        # Pop state off the queue
        node_id = popfirst!(queue)
        node = search_tree[node_id]
        # Return if max nodes are reached or goals are satisfied
        if is_goal(spec, domain, node.state)
            return :success, node_id
        elseif count >= planner.max_nodes
            return :max_nodes, node_id
        end
        count += 1
        # Expand current node
        expand!(planner, node, search_tree, queue, domain, spec)
    end
    return :failure, nothing
end

function expand!(planner::BreadthFirstPlanner, node::PathNode,
                 search_tree::Dict{UInt,PathNode}, queue::Vector{UInt},
                 domain::Domain, spec::Specification)
    state = node.state
    # Iterate over available actions
    actions = available(state, domain)
    for act in actions
        # Execute actions on state
        next_state = transition(domain, state, act; check=false)
        next_id = hash(next_state)
        # Skip if state has already been encountered
        if haskey(search_tree, next_id) continue end
        # Check if next state satisfies trajectory constraints
        if is_violated(spec, domain, state) continue end
        # Update backpointer and add next state to queue
        search_tree[next_id] = PathNode(next_id, next_state, 0, node.id, act)
        push!(queue, next_id)
    end
end


text:
"""
    PoptartCtx

Poptart code generation context.
"""
mutable struct PoptartCtx
    arg_inputs::OrderedDict{Arg, Symbol}
    option_inputs::Dict{Option, Symbol}
    flag_inputs::Dict{Flag, Symbol}
    leaf_windows::Dict{LeafCommand, Symbol}
    app::Symbol
    warning::Symbol
    help::Symbol
    version::Symbol
end

PoptartCtx() = PoptartCtx(OrderedDict{Arg, Symbol}(),
                          Dict{Option, Symbol}(),
                          Dict{Flag, Symbol}(),
                          Dict{LeafCommand, Symbol}(),
                          gensym(:app),
                          gensym(:warning),
                          gensym(:help),
                          gensym(:version))

"""
    Empty arguments-inputs mapping of context

It should be called for each window/casted function
"""
function empty_inputs!(ctx::PoptartCtx)
    ctx.arg_inputs = OrderedDict{Arg, Symbol}()
    ctx.option_inputs = Dict{Option, Symbol}()
    ctx.flag_inputs = Dict{Flag, Symbol}()
    ctx
end

"""
    codegen(cmd)

Generate Julia AST from given command object `cmd`. This will wrap
all the generated AST in a function `command_main`.
"""
function codegen(cmd::AbstractCommand)
    defs = Dict{Symbol,Any}()
    defs[:name] = :poptart_main
    defs[:args] = []

    ctx = PoptartCtx()
    defs[:body] = quote
        $(codegen(ctx, cmd))
    end

    ret_app = gensym(:app)

    return quote
        import Poptart
        $(poptart_compat())

        $(combinedef(defs))
        $ret_app = poptart_main()

        $(xexit(ret_app;esc=false))
    end
end

function codegen(ctx::PoptartCtx, cmd::EntryCommand)
    quote
        $(codegen_app(ctx, cmd_name(cmd)))
        $(codegen_body(ctx, cmd.root))
        return $(ctx.app)
    end
end

function codegen_body(ctx::PoptartCtx, cmd::LeafCommand)
    ret = Expr(:block)

    push!(ret.args, codegen_window(ctx, cmd), codegen_description(ctx, cmd))

    for args in (cmd.args, cmd.options, cmd.flags)
        expr = codegen_controls(ctx, args; cmd=cmd)
        push!(ret.args, expr)
    end

    button_run = gensym(:button_run)
    button_cancel = gensym(:button_cancel)

    params = gensym(:params)
    kwparams = gensym(:kwparams)

    run_event = quote
        $(ctx.warning).text = ""
        try
            $(codegen_params(ctx, params, kwparams, cmd))
            $(codegen_call(ctx, params, kwparams, cmd))
        catch e
            $(ctx.warning).text = string(e)
        end
    end

    cancel_event = xpoptart_desktop(:pause, ctx.app)

    button_expr = quote
        $(XButton(button_run, run_event, title = "run"))
        $(XButton(button_cancel, cancel_event, title = "cancel"))
        $(XLabel(ctx.warning, text = ""))
        $(xwindow_add_item(ctx, cmd, button_run, XSameLine(), button_cancel, ctx.warning))
    end
    push!(ret.args, button_expr)

    return ret
end

function codegen_body(ctx::PoptartCtx, cmd::NodeCommand)

    subcmd = gensym(:subcmd)
    ret = Expr(:block)

    for (i, subcmd) in enumerate(cmd.subcmds)
        push!(ret.args, codegen_body(ctx, subcmd))
        empty_inputs!(ctx)
    end

    push!(ret.args, codegen_entrywindow(ctx, cmd))

    ret
end

function codegen_entrywindow(ctx::PoptartCtx, cmd::NodeCommand)
    ret = Expr(:block)

    buttons = Symbol[]
    for (sub_cmd, window) in ctx.leaf_windows
        btn = gensym(:button)
        push!(ret.args, XWindowButton(btn, cmd_name(sub_cmd), window))
        push!(buttons, btn)
    end

    entry_window = gensym(:entrywindow)
    push!(ret.args, XWindow(entry_window, items=buttons, title=cmd_name(cmd)))
    push!(ret.args, :(push!($(ctx.app).windows, $entry_window)))
    ret
end

function codegen_app(ctx::PoptartCtx, appname::AbstractString="Comonicon")
    XApplication(ctx.app, title=appname, windows=[])
end

function codegen_window(ctx::PoptartCtx, cmd::LeafCommand)
    cmd_window = gensym(:window)
    ctx.leaf_windows[cmd] = cmd_window
    quote
        $(XWindow(cmd_window, title=cmd_name(cmd)))
        push!($(ctx.app).windows, $cmd_window)
    end
end

function codegen_description(ctx::PoptartCtx, cmd::LeafCommand)
    xwindow_add_item(ctx, cmd, XLabel(text = repr(cmd.doc)), XSeparator())
end

function codegen_params(ctx::PoptartCtx, params::Symbol, kwparams::Symbol, cmd::LeafCommand)
    hasparameters(cmd) || return
    args = gensym(:args)
    arg = gensym(:arg)

    ret = quote
        $params = []
        $(xget_args(ctx, args, ctx.arg_inputs))

        for $arg in $args
            if $arg === ""
                break
            end
            push!($params, $arg)
        end
        $kwparams = []
        $args = $(xget_kwargs(ctx, ctx.option_inputs))
        for $arg in $args
            if $arg.second === ""
                continue
            end
            push!($kwparams, $arg)
        end
        $args = $(xget_kwargs(ctx, ctx.flag_inputs))
        for $arg in $args
            push!($kwparams, $arg)
        end
    end
    ret
end

function xget_vararg(::PoptartCtx, args::Symbol, inputgroup::Symbol)
    input = gensym(:input)
    value = gensym(:value)
    quote
        for $input in $inputgroup.items
            $value = $input.buf
            push!($args, $value)
        end
    end
end

function xget_args(ctx::PoptartCtx, args::Symbol, arg_inputs::AbstractDict)
    ret = quote
        $args = []
    end
    params = ret.args[2].args[2]
    for (arg, input) in arg_inputs
        if arg.vararg
            push!(ret.args, xget_vararg(ctx, args, input))
        else
            value = inputvalue(ctx, arg, input)
            push!(params.args, value)
        end

    end
    ret
end

function xget_args(ctx::PoptartCtx, arg_inputs)
    ret = Expr(:vect)
    for (arg, input) in arg_inputs
        if arg.vararg
            continue
        end
        value = inputvalue(ctx, arg, input)
        push!(ret.args, value)

    end
    ret
end

function xget_kwargs(ctx::PoptartCtx, arg_inputs::AbstractDict)
    ret = :(Dict())
    for (arg, input) in arg_inputs
        value = inputvalue(ctx, arg, input)
        push!(ret.args, :($(QuoteNode(cmd_sym(arg))) => $value))
    end
    ret
end

function codegen_call(::PoptartCtx, params::Symbol, kwparams::Symbol, cmd::LeafCommand)
    ex_call = Expr(:call, cmd.entry)

    return :($(cmd.entry)($params...; $kwparams...))
end

function xwarn(ctx::PoptartCtx, message::AbstractString)
    :($(ctx.warning).Label = $message)
end

inputvalue(ctx::PoptartCtx, opt::Option, input::Symbol) = inputvalue(ctx, opt.arg, input)

function inputvalue(::PoptartCtx, arg::Arg, input::Symbol)
    if arg.type == Any
        value = :($input.buf)
    else
        value = :(parse($(arg.type), $input.buf))
    end
    value
end

function inputvalue(::PoptartCtx, flag::Flag, input::Symbol)
    :($input.value)
end

function Base.push!(ctx::PoptartCtx, arg_input::Pair{Arg, Symbol})
    push!(ctx.arg_inputs, arg_input)
end

function Base.push!(ctx::PoptartCtx, arg_input::Pair{Option, Symbol})
    push!(ctx.option_inputs, arg_input)
end

function Base.push!(ctx::PoptartCtx, arg_input::Pair{Flag, Symbol})
    push!(ctx.flag_inputs, arg_input)
end

function process_default(arg::Arg)
    process_default(arg.default)
end

function process_default(opt::Option)
    process_default(opt.arg.default)
end

function process_default(::Nothing)
    (buf="", tip="")
end

function process_default(value::Number)
    (buf=string(value), tip="")
end

function process_default(value::Expr)
    (buf="", tip="\nDefault value is $value")
end

function process_default(value::String)
    (buf="", tip="\nDefault value is $value")
end

function process_default(value)
    process_default(string(value))
end

function process_label(opt::Option)
    process_label(opt.arg, opt.doc)
end

function process_label(arg::Arg)
    process_label(arg, arg.doc)
end

function process_label(flag::Flag)
    label = flag.name * "::Bool"
    doc = flag.doc
    arg_docstring = string(doc)
    if arg_docstring != ""
        label *= "\n" * arg_docstring
    end
    label
end

function process_label(arg::Arg, doc::CommandDoc)
    label = arg.name
    if arg.type != Any
        label *= "::" * string(arg.type)
    end
    if arg.require
        label *= " *"
    end
    arg_docstring = string(doc)
    if arg_docstring != ""
        label *= "\n" * arg_docstring
    end
    label
end

"""
    Generate inputs for `args`
"""
function codegen_controls(ctx::PoptartCtx, args; cmd::LeafCommand)
    genexpr = Expr(:block)
    group = gensym(:group)
    push!(genexpr.args, quote
        $(XGroup(group, items = [XNewLine()]))
        $(xwindow_add_item(ctx, cmd, group))
    end
    )
    for arg in args
        expr, input = codegen_control(ctx, group, arg)
        push!(genexpr.args, expr)
    end
    return genexpr
end

function codegen_control(ctx::PoptartCtx, group::Symbol, arg)
    input_symbol = gensym(:input)
    expr = codegen_control(ctx, input_symbol, group, arg)
    push!(ctx, arg=>input_symbol)
    return expr, input_symbol
end

function is_vararg(arg::Arg)
    arg.vararg
end

function is_vararg(opt::Option)
    opt.arg.vararg
end

function codegen_control(ctx::PoptartCtx, input::Symbol, group::Symbol, arg::Union{Arg, Option})
    label = process_label(arg)

    buf, tip = process_default(arg)
    label *= tip
    if is_vararg(arg)
        codegen_control_vararg(ctx, input; name=arg.name, label=label, group=group)
    else
        codegen_control(ctx, input, buf; name=arg.name, label=label, group=group)
    end
end

function codegen_control_vararg(ctx::PoptartCtx, inputgroup::Symbol; name::AbstractString, label::AbstractString, group::Symbol)
    varaddbutton = gensym(:varaddbutton)
    inputcount = gensym(:inputcount)

    addbutton_event = :(
        try
            $inputcount[] += 1
            $(xpush_item(inputgroup, XInputText(label=:(repr($inputcount[])), buf="")))
        catch e
            @error "error from addbutton" e
        end
    )

    quote
        $inputcount = Ref(0)
        $(XButton(varaddbutton, addbutton_event, title="+"))
        $(XGroup(inputgroup, items=[]))
        $(xpush_item(group, XLabel(text=label), XSameLine(), varaddbutton, inputgroup))
    end
end

function codegen_control(ctx::PoptartCtx, input::Symbol, group::Symbol, flag::Flag)
    label = process_label(flag)

    codegen_control(ctx, input, false; name=flag.name, label=label, group=group)
end

function codegen_control(::PoptartCtx, input::Symbol, buf::AbstractString; name::AbstractString, label::AbstractString, group::Symbol)
    quote
        $(XInputText(input, label=name, buf=buf))
        $(xpush_item(group, XLabel(text=label), input))
    end
end

function codegen_control(::PoptartCtx, input::Symbol, value::Bool; name::AbstractString, label::AbstractString, group::Symbol)
    quote
        $(XCheckbox(input, label=name, value=value))
        $(xpush_item(group, XLabel(text=label), input))
    end
end

"""
Get thw window symbol corresponding to the cmd
"""
function xcmd_window(ctx::PoptartCtx, cmd::LeafCommand)
    ctx.leaf_windows[cmd]
end

function xwindow_add_item(ctx::PoptartCtx, cmd::LeafCommand, items...)
    xpush_item(xcmd_window(ctx, cmd), items...)
end

"""
Generate a button that opens a window
"""
function XWindowButton(btn::Symbol, window_name::AbstractString, window::Symbol)
    openwindow = :(open_window($window))
    XButton(btn, openwindow; title=window_name)
end

"""
temporary

Will be removed once poptart is updated
"""
function poptart_compat()
    quote
        function open_window(window::Poptart.Desktop.Window)
            window_isopen(window) && return

            if !window.props[:show_window_closing_widget]
                window.props[:isopen] = C_NULL
            else
                window.props[:isopen] = Ref(true)
            end
            return
        end

        function window_isopen(window::Poptart.Desktop.Window)
            isopen_prop = window.props[:isopen]
            window_isopen(isopen_prop)
        end

        window_isopen(::Ptr{Nothing}) = true
        window_isopen(x::Ref{Bool}) = x[]
    end
end


text:
# ## Many

"""
    Many(block) <: WrapperBlock

`Many` indicates that you can variable number of instances for
`block`. Consider a bounding box detection task where there may be any
number of targets in an image and this number varies for different
samples. The blocks `(Image{2}(), BoundingBox{2}()` imply that there is exactly
one bounding box for every image, which is not the case. Instead you
would want to use `(Image{2}(), Many(BoundingBox{2}())`.
"""
struct Many{B<:AbstractBlock} <: WrapperBlock
    block::B
end

FastAI.checkblock(many::Many, obss) =
    all(checkblock(wrapped(many), obs) for obs in obss)
FastAI.mockblock(many::Many) = [mockblock(wrapped(many)) for _ in 1:rand(1:3)]

function FastAI.encode(enc::Encoding, ctx, many::Many, obss)
    return map(obss) do obs
        encode(enc, ctx, wrapped(many), obs)
    end
end

function FastAI.decode(enc::Encoding, ctx, many::Many, obss)
    return map(obss) do obs
        decode(enc, ctx, wrapped(many), obs)
    end
end

# ## Tests

@testset "Many [block]" begin
    enc = ImagePreprocessing()
    block = Many(Image{2}())
    FastAI.testencoding(enc, block)
end


text:
module TestInplaceMutationWithReferenceables
using LiterateTest.AssertAsTest: @assert
test() = include(joinpath(@__DIR__, "../examples/inplace_mutation_with_referenceables.jl"))
end  # module


text:
"""
    labels(x::ComponentVector)

Get string labels for for each index of a `ComponentVector`. Useful for automatic plot legend labelling.

# Examples
```jldoctest
julia> x = ComponentArray(a=5, b=[(a=(a=20,b=1), b=0), (a=(a=33,b=1), b=0)], c=(a=(a=2, b=[1,2]), b=[1. 2.; 5 6]))
ComponentVector{Float64}(a = 5.0, b = [(a = (a = 20.0, b = 1.0), b = 0.0), (a = (a = 33.0, b = 1.0), b = 0.0)], c = (a = (a = 2.0, b = [1.0, 2.0]), b = [1.0 2.0; 5.0 6.0]))

julia> labels(x)
14-element Array{String,1}:
 "a"
 "b[1].a.a"
 "b[1].a.b"
 "b[1].b"
 "b[2].a.a"
 "b[2].a.b"
 "b[2].b"
 "c.a.a"
 "c.a.b[1]"
 "c.a.b[2]"
 "c.b[1,1]"
 "c.b[2,1]"
 "c.b[1,2]"
 "c.b[2,2]"
```
see also `label2index`
"""
labels(x::ComponentVector) = map(x->x[firstindex(x)+1:end], _labels(x))
labels(x) = map(x->x[firstindex(x):end], _labels(x))

function _labels(x::ComponentVector)
    vcat((".$(key)" .* _labels(x[key]) for key in keys(x))...)
end
_labels(x::AbstractArray{<:ComponentArray}) = vcat(("[$i]" .* _labels(x[i]) for i in eachindex(x))...)
_labels(x::LazyArray) = vcat(("[$i]" .* _labels(x[i]) for i in eachindex(x))...)
_labels(x::AbstractArray) = vcat(("[" * join(i.I, ",") * "]" for i in CartesianIndices(x))...)
_labels(x) = ""

"""
    label2index(x::ComponentVector, str::AbstractString)
    label2index(label_array, str::AbstractString)

Convert labels made by `labels` function to an array of flat indices of a `ComponentVector`.

# Examples
```jldoctest
julia> x = ComponentArray(a=5, b=[(a=(a=20,b=1), b=0), (a=(a=33,b=1), b=0)], c=(a=(a=2, b=[1,2]), b=[1. 2.; 5 6]))
ComponentVector{Float64}(a = 5.0, b = [(a = (a = 20.0, b = 1.0), b = 0.0), (a = (a = 33.0, b = 1.0), b = 0.0)], c = (a = (a = 2.0, b = [1.0, 2.0]), b = [1.0 2.0; 5.0 6.0]))

julia> labels(x)
14-element Array{String,1}:
 "a"
 "b[1].a.a"
 "b[1].a.b"
 "b[1].b"
 "b[2].a.a"
 "b[2].a.b"
 "b[2].b"
 "c.a.a"
 "c.a.b[1]"
 "c.a.b[2]"
 "c.b[1,1]"
 "c.b[2,1]"
 "c.b[1,2]"
 "c.b[2,2]"

julia> label2index(x, "c.a")
3-element Array{Int64,1}:
  8
  9
 10

julia> label2index(x, "b[1]")
3-element Array{Int64,1}:
 2
 3
 4
```

see also `labels`
"""
label2index(x::ComponentVector, str) = label2index(labels(x), str)
function label2index(labs, str)
    idx = findall(startswith.(labs, Regex("\\Q$str\\E(?:(\\.|\\[))"))) #str * r"(\.|\[)"))
    if !isempty(idx)
        return idx
    else
        return [findfirst(l -> l .== str, labs)]
    end
end


text:
using DiffEqFlux
using Plots

function autocatlysis!(du, u, (a, k₁, k₋₁), t)
    du .= k₁ * a * u - k₋₁ * u.^2
end

p = (10.0, 0.1, 0.01)

u0 = [0.1]
tspan = (0.0, 10.0)
a, k₁, k₋₁ = p
prob = ODEProblem(autocatlysis!,u0,tspan,p)
sol = solve(prob)
plot(sol)


text:
#=
Parallel WE estimation of the probability of for a diffusion with X(0) = x₀
satisfying X(T) ∈ B for the Muller potential.
=#

using Distributed
using Statistics
using Printf

# addprocs(4);

@everywhere include("muller_setup.jl");
@everywhere using WeightedEnsemble

# number of coarse steps in WE
n_we_steps = 10;
# number of time steps during mutation step
nΔt_coarse = ceil(Int, nΔt/n_we_steps);
# number of samples in coarse matrix
n_samples_per_bin = 10^2;
# ensemble size
n_particles = 10^2;

# define bin structure
@everywhere xc = LinRange(-1.5,1,7)
@everywhere yc = LinRange(-0.5,2,7)
@everywhere voronoi_pts = Array{Float64,1}[];
@everywhere for x in xc, y in yc
    if(V([x,y])<250)
        push!(voronoi_pts, [x,y])
    end
end
@everywhere B₀, bin_id, rebin! = setup_Voronoi_bins(voronoi_pts);

# define the mutation mapping
opts = MDOptions(n_iters=nΔt_coarse, n_save_iters = nΔt_coarse)
@everywhere function mutation(x)
    Xvals, _ = sample_trajectory(x, sampler, options=opts);
    return Xvals[end]
end
@everywhere mutation! = x-> sample_trajectory!(x, sampler, options=opts);

Random.seed!(100);
x0_vals = copy(voronoi_pts);
n_bins = length(B₀);
K̃ = WeightedEnsemble.pbuild_coarse_transition_matrix(mutation!, bin_id, x0_vals, n_bins, n_samples_per_bin);

# define coarse observable as a bin function
f̃ = f.(voronoi_pts);
_,v²_vectors = WeightedEnsemble.build_coarse_vectors(n_we_steps,K̃,float.(f̃));
v² = (x,t)-> v²_vectors[t+1][bin_id(x)]
selection! = (E, B, t)-> optimal_selection!(E, B, v², t)

we_sampler = DistributedWEsampler(mutation, selection!, rebin!);

# set up ensemble
E₀ = Dirac_to_Ensemble(x₀, n_particles);
rebin!(E₀, B₀, 0);

# run
E = deepcopy(E₀);
B = deepcopy(B₀);
Random.seed!(200)
prun_we!(E, B, we_sampler, n_we_steps);
p_est = f.(E.ξ) ⋅ E.ω
@printf("WE Estimate = %g\n", p_est)


text:
using Documenter
using Molly

makedocs(
    sitename = "Molly.jl",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home"                      => "index.md",
        "Documentation"             => "docs.md",
        "Differentiable simulation" => "differentiable.md",
        "Examples"                  => "examples.md",
        "API"                       => "api.md",
    ]
)

deploydocs(
    repo="github.com/JuliaMolSim/Molly.jl.git",
    push_preview = true
)


text:
# this is a little Monte Carlo exercise that illustrates that
# the OLS estimator is biased and inconsistent when errors are
# correlated with regressors, but that the IV estimator is consistent
using Econometrics, Plots, LinearAlgebra
function main()
reps = 1000 # number of Monte Carlo reps.
betaols = zeros(reps,2)
betaiv = zeros(reps,2)
n = 1000 # sample size

# covariance of X, W, e
cor_X_W = 0.2  # experiment with lowering or raising this: quality of instrument
cor_X_e = 0.5  # try setting this to zero to observe inefficiency of IV
sig = [
      1.0       cor_X_W     cor_X_e;
      cor_X_W   1.0         0.0;
     cor_X_e    0.0         1.0]
truebeta = [1, 2] # true beta
p = cholesky(sig).U
for i = 1:reps
	XWE = randn(n,3)*p
	e = XWE[:,3:3]
	w = [ones(n) XWE[:,2:2]]
	x = [ones(n) XWE[:,1:1]]
	y = x*truebeta + e
	# OLS
	betaols[i,:] = (x\y)'
	# IV
	betaiv[i,:] = (inv(w'*x)*w'*y)'
end

p1 = npdensity(betaols[:,2])
plot!(p1,title="OLS")
p2 = npdensity(betaiv[:,2])
plot!(p2,title="IV")
p = plot(p1,p2,layout=(2,1))
#savefig("olsiv.png")
gui()
println("true betas are ", truebeta)
println("OLS results")
dstats(betaols, short=true)
println("IV results")
dstats(betaiv, short=true)
return p
end
main()

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
using LiteQTL
using DelimitedFiles
using CSV
using DataFrames

function main()
    geno_file = joinpath(pathof(LiteQTL),"../..", "data", "processed", "spleen-bxd-genoprob.csv")
    pheno_file = joinpath(pathof(LiteQTL), "../..", "data","processed", "spleen-pheno-nomissing.csv")
    gmap_file = joinpath(pathof(LiteQTL), "../..","data","processed", "gmap.csv")

    export_matrix = false
    output_file = "output.csv"
    rqtl_file = joinpath(pathof(LiteQTL), "../..", "data", "UTHSC_SPL_RMA_1210.zip")

    LiteQTL.set_blas_threads(16);
    # Read in data.
    G = get_geno_data(geno_file, Float64)
    Y = get_pheno_data(pheno_file, Float64, transposed=false)
    # getting geno and pheno file size.
    n = size(Y,1)
    m = size(Y,2)
    p = size(G,2)
    println("******* Indivuduals n: $n, Traits m: $m, Markers p: $p ****************");
    # cpu_timing = benchmark(5, scan, Y, G, n; export_matrix);
    # println("CPU timing: $(cpu_timing[3])")

    # running analysis without covariates.
    @time lodc = scan(Y, G, export_matrix=export_matrix, maf_threshold=0.00, usegpu=false, lod_or_pval="lod");
    # lodg = scan(Y, G; usegpu=true)

    if !export_matrix
        gmap = CSV.read(gmap_file, DataFrame)
        idx = trunc.(Int, lodc[:,1])
        gmap_lod = hcat(gmap[idx,:], DataFrame(lodc, [:idx, :maxlod]))
    end

    # write output to file
    # writedlm(joinpath(Base.@__DIR__, "..", "data", "results", output_file), lod, ',')

    # TODO: generate plot?
    # return lod

end

lod = main()


text:
using JLD
using DiffPrecTest
using DataFrames, DataFramesMeta, CSV
using LaTeXStrings, LaTeXTabulars

cd("exper/sim1")

pArr = [100, 200, 500]
elemArr = [(5,5), (8, 7), (50, 25)]
methodArr = ["Sym-N", "Asym-N", "YinXia", "Sym-B", "Asym-B", "O-Sym-N", "O-Asym-N", "O-Sym-B", "O-Asym-B"]

df = DataFrame(p = Int[], row = Int[], col = Int[], method=String[], bias=Float64[], coverage=Float64[], lenCI=Float64[])

for ip=1:2
  for iElem=1:3
    res = load("results/sim1_res_$(ip)_$(iElem).jld", "results")

    for j in 1:9
      global pArr
      global elemArr
      global methodArr
      global df

      sr = computeSimulationResult([res[j, i] for i=1:1000], 0.)
      push!(df, [pArr[ip], elemArr[iElem][1], elemArr[iElem][2], methodArr[j], sr.bias, sr.coverage, sr.lenCoverage])
    end
  end
end

CSV.write("sim1.csv", df)

### create a table

latex_tabular("sim1_table.tex",
  Tabular("lllccc"),
  [
    Rule(:top),
    ["", "", "", "Coverage", "Length", L"{\rm Bias} \times 10^3"],
    Rule(:mid),
    vcat([
      hcat(
        (a = fill("", 9, 1);
        b = latexstring("p = $(pArr[ip])");
        a[1] = "\\multirow{9}{*}{$(b)}";
        a
        ),
        vcat([
          hcat(
            (a = fill("", 3, 1);
            b = latexstring("\\Delta_{$(elemArr[iElem][1]),$(elemArr[iElem][2])} = 0");
            a[1] = "\\multirow{3}{*}{$(b)}";
            a
            ),
            Matrix{Any}(
              @linq df |>
              where( (:method .== "Sym-N") .| (:method .== "YinXia") .| (:method .== "O-Sym-N"),
                    :p .== pArr[ip],
                    :row .== elemArr[iElem][1], :col .== elemArr[iElem][2]) |>  select(:method,
                    coverage = convert.(Int, :coverage * 1000),
                    lenCI = round.(:lenCI, digits=3),
                    bias = round.(:bias * 1000, digits=1))
            )
          )
          for iElem=1:3
          ]...
        )
      )
      for ip=1:2
    ]...),
    Rule(:bottom)
  ]
)


text:
export readdirpath

@doc """ return full path of files under dir
""" ->
function readdirpath(dir;recursive=false, sorted=false)
    if recursive
        #files = filterreaddirpath(dir)
        #return  vcat()
        error("readdirpath recursively is not implemented yet!")
    end
    file_paths = map(x->joinpath(dir,x ), readdir(dir))
    if sorted
        return sort(file_paths)
    end
    return file_paths
end


text:
include(joinpath(pwd(), "examples/direct_policy_optimization/simulate.jl"))

# Unpack trajectories
x̄, ū = unpack(z̄_nom, prob_nominal)
x, u = unpack(z, prob_dpo.prob.prob.nom)
Θ = get_policy(z, prob_dpo)

# Simulation setup
model_sim = model_slosh = RocketSlosh{RK3, FixedTime}(n_slosh, m, n_slosh, g, m1 - m2, l1, J, m2, l2, l3)

x1_sim = copy(x1_slosh)
T_sim = 10 * T

using Random
Random.seed!(1)

W = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(1.0e-5 * ones(model_sim.n)))
w = rand(W, T_sim)

W0 = Distributions.MvNormal(zeros(model_sim.n),
	Diagonal(1.0e-5 * ones(model_sim.n)))
w0 = rand(W0, 1)

z0_sim = vec(copy(x1_sim) + w0)

tf_nom = sum([ū[t][end] for t = 1:T-1])
t_nom = range(0, stop = tf_nom, length = T)
t_sim_nom = range(0, stop = tf_nom, length = T_sim)

tf_dpo = sum([u[t][end] for t = 1:T-1])
t_dpo = range(0, stop = tf_dpo, length = T)
t_sim_dpo = range(0, stop = tf_dpo, length = T_sim)

dt_sim_nom = tf_nom / (T_sim - 1)
dt_sim_dpo = tf_dpo / (T_sim - 1)

# Simulate
z_lqr, u_lqr, J_lqr, Jx_lqr, Ju_lqr = simulate(
	model_sim,
	policy, K,
    x̄, ū,
	Q, R,
	T_sim, ū[1][end],
	z0_sim,
	w,
	_norm = 2,
	ul = ul[1], uu = uu[2],
	u_idx = (1:model_sl.m - 1))

z_dpo, u_dpo, J_dpo, Jx_dpo, Ju_dpo = simulate(
	model_sim,
	policy, Θ,
    x, u,
	Q, R,
	T_sim, u[1][end],
	z0_sim, w,
	_norm = 2,
	ul = ul[1], uu = uu[1],
	u_idx = (1:model_sl.m - 1))

# state tracking
Jx_lqr
Jx_dpo

# control tracking
Ju_lqr
Ju_dpo

# objective value
J_lqr
J_dpo


text:
__precompile__(false)
module ElTopo

using Cxx
using Parameters
using GeometryTypes

include("../deps/deps.jl")

"""
    SurfTrack([...])

Constructs a SurfTrack instance which holds parameters for ElTopo surface tracker stabilization.
"""
@with_kw struct SurfTrack
    "Elements closer than this are considered \"near\" (or proximate)"
    proximity_epsilon::Float64 = 1e-4
    friction_coefficient::Float64 = 0.0
    min_triangle_area::Float64 = 1e-7

    " Collision epsilon to use during mesh improvment operations (i.e. if any mesh elements are closer than this, the operation is aborted). NOTE: This should be greater than collision_epsilon, to prevent improvement operations from moving elements into
     a collision configuration."
    improve_collision_epsilon::Float64 = 2e-6

    "Whether to set the min and max edge lengths as fractions of the initial average edge length"
    use_fraction::Bool = false

    # If use_fraction is true, these are taken to be fractions of the average edge length of the new surface. If use_fraction is false, these are absolute.
    min_edge_length::Float64 = 0.05
    max_edge_length::Float64 = 0.2
    max_volume_change::Float64 = 0.1

    "In-triangle minimum angle to enforce"
    min_triangle_angle::Float64 = 0
    "In-triangle maximum angle to enforce"
    max_triangle_angle::Float64 = 180

    use_curvature_when_splitting::Bool = false
    use_curvature_when_collapsing::Bool = false

    "Clamp curvature scaling with minimum bound"
    min_curvature_multiplier::Float64 = 1
    "Clamp curvature scaling with maximum bound"
    max_curvature_multiplier::Float64 = 1

    allow_vertex_movement::Bool = false

    "Minimum edge length improvement in order to flip an edge"
    edge_flip_min_length_change::Float64 = 0.05

    "Elements within this distance will trigger a merge attempt"
    merge_proximity_epsilon::Float64 = 1e-5

    "Whether to enforce collision-free surfaces (including during mesh maintenance operations)"
    collision_safety::Bool = true

    "Whether to allow changes in topology"
    allow_topology_changes::Bool = false

    "Wether to allow non-manifold (edges incident on more than two triangles)"
    allow_non_manifold::Bool = false

    "Whether to allow mesh improvement"
    perform_improvement::Bool = true
end

Libdl.dlopen(eltopo,Libdl.RTLD_GLOBAL)

addHeaderDir(dirname(@__FILE__)*"/../deps/usr/include/",kind=C_System)

cxxinclude("vector")
cxxinclude("subdivisionscheme.h")
cxxinclude("surftrack.h")

function stabilize(points,faces,p::SurfTrack)
    pointsout = Point{3,Float64}[]
    facesout = Face{3,UInt64}[]

    icxx"""
         std::vector<Vec3d> vs;
         std::vector<Vec3st> ts;
         std::vector<double> ms;

    $:(
    for i in 1:length(points)
        p1,p2,p3 = points[i]
        icxx"vs.push_back(Vec3d($p1,$p2,$p3));"
        icxx"ms.push_back(0.5);"
    end
    );

    $:(
    for i in 1:length(faces)
        v1,v2,v3 = faces[i] .- 1
        icxx"ts.push_back(Vec3st($v1,$v2,$v3));"
    end
    );

      SurfTrackInitializationParameters parameters;
      parameters.m_proximity_epsilon = $(p.proximity_epsilon);
      parameters.m_friction_coefficient = $(p.friction_coefficient);
      parameters.m_min_triangle_area = $(p.min_triangle_area);
      parameters.m_improve_collision_epsilon = $(p.improve_collision_epsilon);
      parameters.m_use_fraction = $(p.use_fraction);
      parameters.m_min_edge_length = $(p.min_edge_length);
      parameters.m_max_edge_length = $(p.max_edge_length);
      parameters.m_max_volume_change = $(p.max_volume_change);
      parameters.m_min_triangle_angle = $(p.min_triangle_angle);
      parameters.m_max_triangle_angle = $(p.max_triangle_angle);
      parameters.m_use_curvature_when_splitting = $(p.use_curvature_when_splitting);
      parameters.m_use_curvature_when_collapsing = $(p.use_curvature_when_collapsing);
      parameters.m_min_curvature_multiplier = $(p.min_curvature_multiplier);
      parameters.m_max_curvature_multiplier = $(p.max_curvature_multiplier);
      parameters.m_allow_vertex_movement = $(p.allow_vertex_movement);
      parameters.m_edge_flip_min_length_change = $(p.edge_flip_min_length_change);
      parameters.m_merge_proximity_epsilon = $(p.merge_proximity_epsilon);
      parameters.m_collision_safety = $(p.collision_safety);
      parameters.m_allow_topology_changes = $(p.allow_topology_changes);
      parameters.m_allow_non_manifold = $(p.allow_non_manifold);
      parameters.m_perform_improvement = $(p.perform_improvement);
      parameters.m_subdivision_scheme = new ButterflyScheme();

      SurfTrack surface_tracker(vs,ts,ms,parameters);

    surface_tracker.m_verbose = false;

    surface_tracker.improve_mesh();
    surface_tracker.topology_changes();
    surface_tracker.defrag_mesh();

    for ( int i = 0; i < surface_tracker.get_num_vertices(); ++i )
        $:(
           p1 = icxx"return surface_tracker.get_position(i)[0];";
           p2 = icxx"return surface_tracker.get_position(i)[1];";
           p3 = icxx"return surface_tracker.get_position(i)[2];";
           push!(pointsout,Point(p1,p2,p3));
        );

    for ( int i = 0; i < surface_tracker.m_mesh.num_triangles(); ++i )
        $:(
           v1 = icxx"return surface_tracker.m_mesh.get_triangle(i)[0];";
           v2 = icxx"return surface_tracker.m_mesh.get_triangle(i)[1];";
           v3 = icxx"return surface_tracker.m_mesh.get_triangle(i)[2];";
           push!(facesout,Face(v1+1,v2+1,v3+1));
        );
    """

    return pointsout,facesout
end

"""
    stabilize(msh::HomogenousMesh,p::SurfTrack)

Computes stabilized mesh for a given initial triangulation `msh` with stabilization parameters `p`.
"""
function stabilize(msh::HomogenousMesh,parameters::SurfTrack)
    newvertices, newfaces = stabilize(msh.vertices,msh.faces,parameters)
    return HomogenousMesh(newvertices,newfaces)
end

export stabilize, SurfTrack

end


text:
using Knet
using ArgParse
using AutoGrad
using GZip
using Images
using ImageMagick

function main(args)
    s = ArgParseSettings()
    s.description = string(
        "General Adversarial Networks implementation in Knet on MNIST data.",
        " by Ilker Kesen [ikesen16_at_ku.edu.tr], 2017.")

    @add_arg_table s begin
        ("--outdir"; default=nothing; help="generations save dir")
        ("--nogpu"; action=:store_true)
        ("--dnet"; nargs='+'; arg_type=Int; default=[128])
        ("--gnet"; nargs='+'; arg_type=Int; default=[128])
        ("--zdim"; arg_type=Int64; default=100)
        ("--winit"; arg_type=Float32; default=Float32(0.01))
        ("--epochs"; arg_type=Int; default=100)
        ("--batchsize"; arg_type=Int; default=100)
        ("--lr"; arg_type=Float32; default=Float32(0.001))
        ("--adam"; action=:store_true; help="adam optimizer")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gcheck"; arg_type=Int; default=0; help="gradient checking")
        ("--gridsize"; arg_type=Int64; nargs=2; default=[8,8])
        ("--gridscale"; arg_type=Float64; default=2.0)
        ("--activations"; nargs=3; default=["tanh","tanh","tanh"])
        ("--dropouts"; arg_type=Float64; nargs=3; default=[0.5,0.5,0.5])
    end

    # parse args
    @printf("\nScript started. [%s]\n", now()); flush(STDOUT)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true); println(o); flush(STDOUT)
    o[:seed] > 0 && srand(o[:seed])
    o[:activations] = map(x->eval(parse(x)), o[:activations])
    xscale = o[:activations][end] == tanh ? (255/2) : 255
    xnorm  = o[:activations][end] == tanh ? 1 : 0

    # load data
    (xtrn,xtst,ytrn,ytst)=loaddata()
    trn = minibatch(xtrn, ytrn, o[:batchsize]; xscale=xscale, xnorm=xnorm)
    tst = minibatch(xtst, ytst, o[:batchsize]; xscale=xscale, xnorm=xnorm)

    # get parameters
    atype = !o[:nogpu] ? KnetArray{Float32} : Array{Float32}
    xdim = size(trn[1][1],1)
    wd, wg = initweights(atype, xdim, o[:zdim], o[:dnet], o[:gnet])

    # gradient check
    if o[:gcheck] > 0
        x = convert(atype, trn[1][1])
        z = sample_noise(atype, size(trn[1][1],2), o[:zdim])
        gradcheck(loss, wd, wg, x, z; gcheck=o[:gcheck])
        gradcheck(loss, wg, wd, z; gcheck=o[:gcheck])
    end

    if o[:outdir] != nothing
        o[:outdir] = abspath(o[:outdir])
        !isdir(o[:outdir]) && mkdir(o[:outdir])
    end

    # initialize optimization params, using ADAM
    optd = initopt(wd,o[:lr])
    optg = initopt(wg,o[:lr])
    sample(x) = sample_noise(atype, size(x,2), o[:zdim])

    # performance with random initialization
    loss1, loss2 = test(wd,wg,trn,o)
    @printf("\nepoch: %d, losses: %g/%g\n", 0, loss1, loss2)
    if o[:outdir] != nothing
        out = generate(wg,o)
        png = makegrid(out; gridsize=o[:gridsize], scale=o[:gridscale])
        filename = @sprintf("%04d.png",0)
        save(joinpath(o[:outdir],filename), png)
    end

    # training, loss1 => discriminator, loss2 => generator
    for epoch = 1:o[:epochs]
        shuffle!(trn)
        loss1 = loss2 = 0
        for i = 1:length(trn)
            x = convert(atype, trn[i][1])
            z = convert(atype, sample(x))

            losses = train!(wd,wg,x,z,o,optd,optg)
            loss1 += losses[1]; loss2 += losses[2]
        end

        loss1, loss2 = test(wd,wg,trn,o)
        @printf("epoch: %d, losses: %g/%g\n", epoch, loss1, loss2)
        flush(STDOUT)
        if o[:outdir] != nothing
            out = generate(wg,o)
            png = makegrid(out; gridsize=o[:gridsize], scale=o[:gridscale])
            filename = @sprintf("%04d.png",epoch)
            save(joinpath(o[:outdir],filename), png)
        end
    end
end

function dropout(x,d)
    if d > 0
        return x .* (rand!(similar(AutoGrad.getval(x))) .> d) * (1/(1-d))
    else
        return x
    end
end

function D(w,x; f1=tanh, pdrop=0.0)
    for k = 1:2:length(w)-2
        x = f1(w[k] * x .+ w[k+1])
        x = dropout(x, pdrop)
    end
    x = w[end-1] * x .+ w[end]
end

function G(w,z; f2=tanh, f3=tanh, pdrop=0.5)
    x = z
    for k = 1:2:length(w)-2
        x = f2(w[k] * x .+ w[k+1])
        x = dropout(x, pdrop)
    end
    x = f3(w[end-1] * x .+ w[end])
end

function initweights(atype, xdim, zdim, dnet, gnet)
    # discriminator weights
    wd = Array(Any, 2*(1+length(dnet)))
    x = xdim
    for (k,y) in enumerate([dnet..., 2])
        wd[2k-1] = xavier(y,x)
        wd[2k] = zeros(y,1)
        x = y
    end

    # generator weights
    wg = Array(Any, 2*(1+length(gnet)))
    x = zdim
    for (k,y) in enumerate([gnet..., xdim])
        wg[2k-1] = xavier(y,x)
        wg[2k] = zeros(y,1)
        x = y
    end

    map(x->convert(atype, x), wd), map(x->convert(atype, x), wg)
end

# loss for discriminator network
function dloss(wd,wg,x,z,labels,fs,pdrops,values=[])
    real = x
    fake = G(wg,z; pdrop=pdrops[2], f2=fs[2], f3=fs[3])
    data = hcat(real,fake)
    pred = logp(D(wd,data; pdrop=pdrops[1], f1=fs[1]), 1)
    val  = -sum(pred .* labels) / size(data,2)
    push!(values, val)
    return val
end

# loss for generator network
function gloss(wg,wd,z,labels,fs,pdrops,values=[])
    data = G(wg,z; pdrop=pdrops[2], f2=fs[2], f3=fs[3])
    pred = logp(D(wd,data; pdrop=pdrops[1], f1=fs[1]), 1)
    val  = -sum(pred .* labels) / size(data,2)
    push!(values, val)
    return val
end

dlossgradient = grad(dloss)
glossgradient = grad(gloss)

function sample_noise(atype, batchsize, dimension, scale=1)
    return convert(atype, scale*randn(dimension,batchsize))
end

initopt(w,lr) = map(x->Adam(;lr=lr), w)

function train!(wd,wg,x,z,o,optd,optg)
    val = []
    atype = typeof(AutoGrad.getval(wd[1]))
    yx = zeros(2,size(x,2)); yx[1,:] = 1; yx = convert(atype, yx)
    yz = zeros(2,size(z,2)); yz[2,:] = 1; yz = convert(atype, yz)
    g = dlossgradient(wd,wg,x,z,hcat(yx,yz),o[:activations],o[:dropouts],val)
    for k = 1:length(wd); update!(wd[k], g[k], optd[k]); end
    g = glossgradient(wg,wd,z,yx,o[:activations],o[:dropouts],val)
    for k = 1:length(wg); update!(wg[k], g[k], optg[k]); end
    return val
end

function test(wd,wg,data,o)
    atype = !o[:nogpu] ? KnetArray{Float32} : Array{Float32}
    sample(x) = sample_noise(atype, size(x,2), o[:zdim])
    loss1 = loss2 = 0
    for (x,y) in data
        x = convert(atype, x)
        z = convert(atype, sample(x))
        yx = zeros(2,size(x,2)); yx[1,:] = 1; yx = convert(atype, yx)
        yz = zeros(2,size(z,2)); yz[2,:] = 1; yz = convert(atype, yz)
        loss1 += dloss(wd,wg,x,z,hcat(yx,yz),o[:activations],o[:dropouts])
        loss2 += gloss(wg,wd,z,yx,o[:activations],o[:dropouts])
    end
    return (loss1/length(data),loss2/length(data))
end

function generate(wg,o)
    atype = !o[:nogpu] ? KnetArray{Float32} : Array{Float32}
    sample(n) = sample_noise(atype, n, o[:zdim])
    ninstances = o[:gridsize][1] * o[:gridsize][2]
    z = convert(atype, sample(ninstances))
    y = G(wg,z; pdrop=o[:dropouts][end])
    y = convert(Array{Float64}, y)
    y = (y+1)/2
    return min(1,max(0,y))
end

function makegrid(y; gridsize=[4,4], scale=2.0, shape=(28,28))
    y = reshape(y, shape..., size(y,2))
    y = map(x->y[:,:,x]', [1:size(y,3)...])
    shp = map(x->Int(round(x*scale)), shape)
    y = map(x->Images.imresize(x,shp), y)
    gridx, gridy = gridsize
    outdims = (gridx*shp[1]+gridx+1,gridy*shp[2]+gridy+1)
    out = zeros(outdims...)
    for k = 1:gridx+1; out[(k-1)*(shp[1]+1)+1,:] = 1.0; end
    for k = 1:gridy+1; out[:,(k-1)*(shp[2]+1)+1] = 1.0; end

    x0 = y0 = 2
    for k = 1:length(y)
        x1 = x0+shp[1]-1
        y1 = y0+shp[2]-1
        out[x0:x1,y0:y1] = y[k]

        y0 = y1+2
        if k % gridy == 0
            x0 = x1+2
            y0 = 2
        else
            y0 = y1+2
        end
    end

    return convert(Array{Float64,2}, map(x->isnan(x)?0:x, out))
end

function minibatch(
    x, y, batchsize;
    atype=Array{Float32}, xrows=784, yrows=10, xscale=255/2, xnorm=1)
    xbatch(a)=convert(atype, reshape(a./xscale-xnorm, xrows, div(length(a),xrows)))
    ybatch(a)=(a[a.==0]=10; convert(atype, sparse(convert(Vector{Int},a),1:length(a),one(eltype(a)),yrows,length(a))))
    xcols = div(length(x),xrows)
    xcols == length(y) || throw(DimensionMismatch())
    data = Any[]
    for i=1:batchsize:xcols-batchsize+1
        j=i+batchsize-1
        push!(data, (xbatch(x[1+(i-1)*xrows:j*xrows]), ybatch(y[i:j])))
    end
    return data
end

function loaddata()
    info("Loading MNIST...")
    gzload("train-images-idx3-ubyte.gz")[17:end],
    gzload("t10k-images-idx3-ubyte.gz")[17:end],
    gzload("train-labels-idx1-ubyte.gz")[9:end],
    gzload("t10k-labels-idx1-ubyte.gz")[9:end]
end

function gzload(file; path=Knet.dir("data",file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = read(f)
    close(f)
    return(a)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/adversarial

max_issues_repo_name: muhammetfaik/GAN

max_forks_repo_name: muhammetfaik/GAN


text:
"""
    setup_dump_tendencies_diagnostics(
        ::ClimateMachineConfigType,
        interval::String,
        out_prefix::String;
        writer = NetCDFWriter(),
        interpol = nothing,
    )

Create the "DumpTendencies" `DiagnosticsGroup` which contains all the
tendencies for the simulation's `BalanceLaw`. These are collected on
an interpolated grid (`interpol` _must_ be specified).

!!! warn
    This diagnostics group can produce a lot of output and is
    very expensive.

!!! warn
    These diagnostics are intended for physics-debugging, and not numerics debugging, because
    `flux` and `source` are called on interpolated data, which may not exactly match the
    non-interpolated evaluation of fluxes and sources that the solver uses.
"""
function setup_dump_tendencies_diagnostics(
    ::ClimateMachineConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "DumpTendencies",
        Diagnostics.dump_tendencies_init,
        Diagnostics.dump_tendencies_fini,
        Diagnostics.dump_tendencies_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

# helpers for PV and tendency names
params(::T) where {T} = T.parameters
unval(::Val{i}) where {i} = i
param_suffix(p, ::Val{n}) where {n} = string("_", unval(p[1]))
param_suffix(p, ::Val{0}) = ""
param_suffix(p) = param_suffix(p, Val(length(p)))
prog_name(pv::PV) where {PV} = string(nameof(PV), param_suffix(params(pv)))
tend_name(t) = String(nameof(typeof(t)))

function precompute_args(
    bl,
    state,
    aux,
    diffusive,
    hyperdiffusive,
    t,
    direction,
)
    _args_fx1 = (; state, aux, t, direction)
    _args_fx2 = (; state, aux, t, diffusive, hyperdiffusive)
    _args_src = (; state, aux, t, direction, diffusive)

    cache_fx1 = precompute(bl, _args_fx1, Flux{FirstOrder}())
    cache_fx2 = precompute(bl, _args_fx2, Flux{SecondOrder}())
    cache_src = precompute(bl, _args_src, Source())

    args_fx1 = merge(_args_fx1, (; precomputed = cache_fx1))
    args_fx2 = merge(_args_fx2, (; precomputed = cache_fx2))
    args_src = merge(_args_src, (; precomputed = cache_src))

    return (args_fx1, args_fx2, args_src)
end

function dump_tendencies_init(dgngrp, t)
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q)
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if mpirank == 0
        if array_device(Q) isa CPU
            prognostic_data = Q.realdata
            auxiliary_data = dg.state_auxiliary.realdata
            gradient_flux_data = dg.state_gradient_flux.realdata
            hyperdiffusive_data = dg.states_higher_order[2].realdata
        else
            prognostic_data = Array(Q.realdata)
            auxiliary_data = Array(dg.state_auxiliary.realdata)
            gradient_flux_data = Array(dg.state_gradient_flux.realdata)
            hyperdiffusive_data = Array(dg.states_higher_order[2].realdata)
        end

        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)
        dim_names = tuple(collect(keys(dims))...)

        state, aux, diffusive, hyperdiffusive = (
            extract_state(dg, prognostic_data, 1, 1, Prognostic()),
            extract_state(dg, auxiliary_data, 1, 1, Auxiliary()),
            extract_state(dg, gradient_flux_data, 1, 1, GradientFlux()),
            extract_state(dg, hyperdiffusive_data, 1, 1, Hyperdiffusive()),
        )
        direction = dg.direction
        args_fx1, args_fx2, args_src = precompute_args(
            bl,
            state,
            aux,
            diffusive,
            hyperdiffusive,
            t,
            direction,
        )

        # set up the variables we're going to be writing
        vars = OrderedDict()
        for (tend_fun, tend_type, tend_args) in (
            (flux, Flux{FirstOrder}(), args_fx1),
            (flux, Flux{SecondOrder}(), args_fx2),
            (source, Source(), args_src),
        )
            for pv in prognostic_vars(bl)
                progname = prog_name(pv)
                for tend in eq_tends(pv, bl, tend_type)
                    flat_vals = flattened_tuple(
                        FlattenArr(),
                        tend_fun(pv, tend, bl, tend_args),
                    )
                    for i in 1:length(flat_vals)
                        varname = "$(progname)_$(tend_name(tend))_$(i)"
                        vars[varname] = (dim_names, FT, Dict())
                    end
                end
            end
        end

        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        noov = Settings.no_overwrite
        init_data(dgngrp.writer, dfilename, noov, dims, vars)
    end

    return nothing
end

function dump_tendencies_collect(dgngrp, t)
    interpol = dgngrp.interpol
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q.data)
    bl = dg.balance_law
    mpirank = MPI.Comm_rank(mpicomm)

    istate = similar(Q.data, interpol.Npl, number_states(bl, Prognostic()))
    interpolate_local!(interpol, Q.data, istate)
    iaux = similar(Q.data, interpol.Npl, number_states(bl, Auxiliary()))
    interpolate_local!(interpol, dg.state_auxiliary.data, iaux)
    igf = similar(Q.data, interpol.Npl, number_states(bl, GradientFlux()))
    interpolate_local!(interpol, dg.state_gradient_flux.data, igf)
    ihd = similar(Q.data, interpol.Npl, number_states(bl, Hyperdiffusive()))
    interpolate_local!(interpol, dg.states_higher_order[2].data, ihd)

    if interpol isa InterpolationCubedSphere
        i_ρu = varsindex(vars_state(bl, Prognostic(), FT), :ρu)
        project_cubed_sphere!(interpol, istate, tuple(collect(i_ρu)...))
    end

    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    all_aux_data = accumulate_interpolated_data(mpicomm, interpol, iaux)
    all_gf_data = accumulate_interpolated_data(mpicomm, interpol, igf)
    all_hd_data = accumulate_interpolated_data(mpicomm, interpol, ihd)

    if mpirank == 0
        varvals = OrderedDict()
        dims = dimensions(dgngrp.interpol)
        direction = dg.direction
        nx, ny, nz = [length(dim[1]) for dim in values(dims)]
        @traverse_interpolated_grid nx ny nz begin
            vars_view(st, data) =
                Vars{vars_state(bl, st, FT)}(view(data, lo, la, le, :))
            statei = vars_view(Prognostic(), all_state_data)
            auxi = vars_view(Auxiliary(), all_aux_data)
            diffusivei = vars_view(GradientFlux(), all_gf_data)
            hyperdiffusivei = vars_view(Hyperdiffusive(), all_hd_data)

            args_fx1, args_fx2, args_src = precompute_args(
                bl,
                statei,
                auxi,
                diffusivei,
                hyperdiffusivei,
                t,
                direction,
            )

            for (tend_fun, tend_type, tend_args) in (
                (flux, Flux{FirstOrder}(), args_fx1),
                (flux, Flux{SecondOrder}(), args_fx2),
                (source, Source(), args_src),
            )
                for pv in prognostic_vars(bl)
                    progname = prog_name(pv)
                    for tend in eq_tends(pv, bl, tend_type)
                        flat_vals = flattened_tuple(
                            FlattenArr(),
                            tend_fun(pv, tend, bl, tend_args),
                        )
                        for i in 1:length(flat_vals)
                            varname = "$(progname)_$(tend_name(tend))_$(i)"
                            vals = get!(
                                varvals,
                                varname,
                                Array{FT}(undef, nx, ny, nz),
                            )
                            vals[lo, la, le] = flat_vals[i]
                        end
                    end
                end
            end

        end

        append_data(dgngrp.writer, varvals, t)
    end

    MPI.Barrier(mpicomm)
    return nothing
end

function dump_tendencies_fini(dgngrp, t) end