text:
main_paths = ["."]
for main_path in main_paths
    format(
        main_path;
        whitespace_ops_in_indices = true,
        remove_extra_newlines = true,
        verbose = true,
        always_for_in = true,
        whitespace_typedefs = true,
    )
end

ext: jl


text:
module Persa

# package code goes here
include("preference.jl")
include("rating.jl")
include("dataset.jl")
include("learn.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TetGen_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TetGen")
JLLWrappers.@generate_main_file("TetGen", UUID("b47fdcd6-d2c1-58e9-bbba-c1cee8d8c179"))
end  # module TetGen_jll


text:
function wordcount(sentence)

end


text:
myParse(x) = parse.(Int, x)
energyLevels = Matrix(hcat(myParse.(collect.(readlines("input.txt")))...)')

nRows, nCols = size(energyLevels)
nLevels = length(energyLevels)

function coordfromIndex(index)
    row = ((index - 1) % nRows) + 1
    col = ((index - 1) ÷ nRows) + 1
    return row, col
end

function indexfromCoord(row, col)
    return nRows * (col - 1) + row
end

function surroundingIndices(index)
    row, col = coordfromIndex(index)
    arr = [(row - 1, col), (row, col - 1), (row + 1, col), (row, col + 1), (row + 1, col +1), (row - 1, col - 1), (row - 1, col + 1), (row + 1, col - 1)]
    return [indexfromCoord(x, y) for (x, y) in arr if x in 1:nRows && y in 1:nCols]
end

function simulatePt1(steps, elevels)

    function flashNode(index, flashed)
        flashed[index] = true
        for i in surroundingIndices(index)
            elevels[i] += 1
            if elevels[i] > 9 && !flashed[i]
                flashNode(i, flashed)
            end
        end
    end

    function simulateFlashes()
        flashed = [false for i in 1:nLevels]
        for i in 1:nLevels
            if elevels[i] > 9 && !flashed[i]
                flashNode(i, flashed)
            end
        end

        for (i, hasFlashed) in enumerate(flashed)
            if hasFlashed
                elevels[i] = 0
            end
        end

        return count(flashed)
    end

    flashes = 0

    for _ in 1:steps
        elevels .= elevels .+ 1
        flashes += simulateFlashes()
    end

    return flashes
end

simulate(100, copy(energyLevels))

function simulatePt2(elevels)

    function flashNode(index, flashed)
        flashed[index] = true
        for i in surroundingIndices(index)
            elevels[i] += 1
            if elevels[i] > 9 && !flashed[i]
                flashNode(i, flashed)
            end
        end
    end

    function simulateFlashes()
        flashed = [false for i in 1:nLevels]
        for i in 1:nLevels
            if elevels[i] > 9 && !flashed[i]
                flashNode(i, flashed)
            end
        end

        for (i, hasFlashed) in enumerate(flashed)
            if hasFlashed
                elevels[i] = 0
            end
        end

        return count(flashed) == length(flashed)
    end

    step = 1

    while true
        elevels .= elevels .+ 1
        if simulateFlashes() == true
            return step
        end
        step += 1
    end
end

simulatePt2(energyLevels)

max_stars_repo_path: Day 11/main.jl

max_issues_repo_path: Day 11/main.jl

max_forks_repo_path: Day 11/main.jl


text:
using ConditionalMoments
using Documenter

DocMeta.setdocmeta!(ConditionalMoments, :DocTestSetup, :(using ConditionalMoments); recursive=true)

makedocs(;
    modules=[ConditionalMoments],
    authors="William Davis",
    repo="https://github.com/williamjsdavis/ConditionalMoments.jl/blob/{commit}{path}#{line}",
    sitename="ConditionalMoments.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://williamjsdavis.github.io/ConditionalMoments.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/williamjsdavis/ConditionalMoments.jl",
)


text:
using Formatting
using CFTime, Dates
using DataStructures
using NCDatasets

function checkData(
    hist_dir :: String,
    casename :: String,
    varnames :: Array{String, 1},
    year_rng :: Array{Int64, 1};
    verbose :: Bool = false,
)

    expect_filenames = DataStructures.OrderedDict()

    println("I will also check a month prior to the begin year because we are using 12/31 of previous year as initial profile.")

    for y=(year_rng[1]-1):(year_rng[2]+1), m=1:12

        if (y == year_rng[1] - 1 && m != 12) || (y == year_rng[2]+1 && m != 1)
            continue
        end

        key = format("{:04d}-{:02d}", y, m)

        filename = format("$(casename).pop.h.daily.{:04d}-{:02d}-01.nc", y, m)
        full_path = joinpath(hist_dir, filename)
        verbose && println("Checking $(key)")

        if ! isfile(full_path)
            throw(ErrorException("Error: File $(filename) does not exist in folder $(hist_dir)."))
        end

        ds = Dataset(full_path, "r")
        dom = daysinmonth(DateTimeNoLeap, y, m)
        if ds.dim["time"] != dom
            throw(ErrorException("Error: File $(filename) does not have $(dom) records."))
        end

        for varname in varnames
            if ! haskey(ds, varname)
                throw(ErrorException("Error: File $(filename) does not have variable `$(varname)`."))
            end
        end
        close(ds)

        expect_filenames[key] = Dict(
            :fp => full_path,
            :fn => filename,
            :y => y,
            :m => m,
        )

    end

    return expect_filenames

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
abstract type AppResource end

struct AppRelativeResource <: AppResource
    namespace ::String
    relative_path ::String
    version ::String
    ts ::Int64
end
struct AppExternalResource <: AppResource
    url ::String
end
struct AppCustomResource <: AppResource
    params ::Dict{String, String}
end
struct AppAssetResource <: AppResource
    path ::String
    ts ::Int64
end

struct NamespaceFiles
    base_path ::String
    files ::Set{String}
    NamespaceFiles(base_path) = new(base_path, Set{String}())
end

struct ApplicationResources
    files ::Dict{String, NamespaceFiles}
    css ::Vector{AppResource}
    js ::Vector{AppResource}
    favicon ::Union{Nothing, AppAssetResource}
    ApplicationResources(files, css, js, favicon) = new(files, css, js, favicon)
end

function ApplicationResources(app::DashApp, registry::ResourcesRegistry)
    css = AppResource[]
    js = AppResource[]
    favicon::Union{Nothing, AppResource} = nothing
    files = Dict{String, NamespaceFiles}()

    serve_locally = get_setting(app, :serve_locally)
    assets_external_path = get_setting(app, :assets_external_path)
    dev = get_devsetting(app, :serve_dev_bundles)
    eager_loading = get_setting(app, :eager_loading)

    append_pkg = function(pkg)
        append!(css,
            _convert_resource_pkg(pkg, :css, dev = dev, serve_locally = serve_locally, eager_loading = eager_loading)
        )
        append!(js,
            _convert_resource_pkg(pkg, :js, dev = dev, serve_locally = serve_locally, eager_loading = eager_loading)
        )
    end

    asset_resource = serve_locally || isnothing(assets_external_path)  ?
                     (url, modified) -> AppAssetResource(url, modified) :
                     (url, modified) -> AppExternalResource(assets_external_path*url)

    append_pkg(get_dash_dependencies(registry, get_devsetting(app, :props_check)))

    append!(css,
        _convert_external.(get_setting(app, :external_stylesheets))
    )
    append!(js,
        _convert_external.(get_setting(app, :external_scripts))
    )

    if get_setting(app, :include_assets_files)
        walk_assets(app) do type, url, path, modified
            if type == :js
                push!(js, asset_resource(url, modified))
            elseif type == :css
                push!(css, asset_resource(url, modified))
            elseif type == :favicon
                favicon = AppAssetResource(url, modified)
            end
        end
    end

    append_pkg.(get_componens_pkgs(registry))
    append_pkg(get_dash_renderer_pkg(registry))

    fill_files(files, get_dash_dependencies(registry, get_devsetting(app, :props_check)), dev = dev, serve_locally = serve_locally)
    fill_files.(Ref(files), get_componens_pkgs(registry), dev = dev, serve_locally = serve_locally)
    fill_files(files, get_dash_renderer_pkg(registry), dev = dev, serve_locally = serve_locally)
    ApplicationResources(files, css, js, favicon)
end

function fill_files(dest::Dict{String, NamespaceFiles}, pkg::ResourcePkg; dev, serve_locally)
    !serve_locally && return

    if !haskey(dest, pkg.namespace)
        dest[pkg.namespace] = NamespaceFiles(pkg.path)
    end
    namespace_files = get!(dest, pkg.namespace, NamespaceFiles(pkg.path))
        full_path = (p)->joinpath(pkg.path, lstrip(p, ['/','\\']))
    for resource in pkg.resources
        paths = dev && has_dev_path(resource) ? get_dev_path(resource) : get_relative_path(resource)
        for path in paths
            push!.(
                Ref(namespace_files.files),
                lstrip.(path, Ref(['/', '\\']))
            )
        end
    end
end

function walk_assets(callback, app::DashApp)
    assets_ignore = get_setting(app, :assets_ignore)
    assets_regex = Regex(assets_ignore)

    assets_filter = isempty(assets_ignore) ?
        (f) -> true :
        (f) -> !occursin(assets_regex, f)

    assets_path = get_assets_path(app)
    if get_setting(app, :include_assets_files) && isdir(assets_path)
        for (base, dirs, files) in walkdir(assets_path)
            if !isempty(files)
                relative = ""
                if base != assets_path
                    relative = join(
                            splitpath(
                                relpath(base, assets_path)
                            ), "/"
                        ) * "/"
                end
                for file in Iterators.filter(assets_filter, files)
                    file_type::Symbol = endswith(file, ".js") ? :js :
                                (
                                    endswith(file, ".css") ? :css :
                                    (
                                        file == "favicon.ico" ? :favicon : :none
                                    )
                                )

                    if file_type != :none
                        full_path = joinpath(base, file)
                        callback(file_type,
                                relative*file,
                                full_path,
                                trunc(Int64, stat(full_path).mtime)
                                )
                    end
                end

            end

        end
    end
end

_convert_external(v::String) = AppExternalResource(v)
_convert_external(v::Dict{String, String}) = AppCustomResource(v)

function _convert_resource(resource::Resource; namespace, version, ts, dev, serve_locally)::Vector{AppResource}
    if !serve_locally && has_external_url(resource)
        return AppExternalResource.(get_external_url(resource))
    end
    relative_path = dev && has_dev_path(resource) ?
                    get_dev_path(resource) :
                    (has_relative_path(resource) ? get_relative_path(resource) : nothing)
    if !isnothing(relative_path)
        return AppRelativeResource.(namespace, relative_path, version, ts)
    end
    if serve_locally
        #TODO Warning
        #FIXME Why warning? It's strange logic in python§
        return AppResource[]
    end
    error("$(resource) does not have a relative_package_path, absolute_path, or an external_url.")
end

function _convert_resource_pkg(pkg::ResourcePkg, type::Symbol; dev, serve_locally, eager_loading)
    result = AppResource[]
    iterator = Iterators.filter(pkg.resources) do r
        get_type(r) == type && !isdynamic(r, eager_loading)
    end
    ts = ispath(pkg.path) ? trunc(Int64, stat(pkg.path).mtime) : 0
    for resource in iterator
        append!(
            result,
            _convert_resource(
                    resource,
                    namespace = pkg.namespace,
                    version = pkg.version,
                    ts = ts,
                    dev = dev, serve_locally = serve_locally
                )
            )
    end
    return result
end


text:
export start_game, server_handler

function agent_ip(args::Dict)
    string("http://", args["agentIP"], ":", args["agentPort"])
end

function breezy_ip(args::Dict)
    string("http://", args["breezyIP"], ":", args["breezyPort"])
end

function start_game(args::Dict)
    # build url to dota 2 breezy server
    start_url = string(breezy_ip(args), "/run/")
    start_data = Dict("agent"=>"SuReLI agent", "size"=>1)
    # create a run config for this agent
    HTTP.post(start_url, ["Content-Type" => "application/json"], JSON.json(start_data))
end

function server_handler(request::HTTP.Request, args::Dict, get_action::Function)
    # path is either an array containing "update" or nothing
    if occursin("update", request.target)
        println("Game done.")
        content = JSON.parse(String(request.body))

        # webhook to start new game in existing set of games
        if "webhook" in keys(content)
            # A webhook was sent to the agent to start a new game in the current set of games.
            webhook_url = string(breezy_ip(args), content["webhook"])
            # call webhook to trigger new game
            println("Starting a new game among the set")
            HTTP.get(webhook_url)
        # otherwise start new set of games, or end session
        else
            start_game(args)
        end
        # send response to our server
        return HTTP.Response(200, JSON.json(Dict("done"=>true)))
    else
        println("getting features")
        # relay path gives features from current game to agent
        features = Array{Float64}(JSON.parse(String(request.body)))
        action = get_action(features)
        return HTTP.Response(200, JSON.json(Dict("actionCode"=>action)))
    end
end


text:
"""
    LBA(;τ, A, k, ν, σ=1.0)

A model object for the linear ballistic accumulator.

# Fields

- `ν`: a vector of drift rates
- `A`: max start point
- `k`: A + k = b, where b is the decision threshold
- `σ=1`: drift rate standard deviation
- `τ`: a encoding-response offset

# Example

````julia
using SequentialSamplingModels
dist = LBA(ν=[3.0,2.0], A = .8, k = .2, τ = .3)
choice,rt = rand(dist, 10)
like = pdf.(dist, choice, rt)
loglike = logpdf.(dist, choice, rt)
````

# References

Brown, S. D., & Heathcote, A. (2008). The simplest complete model of choice response time: Linear ballistic accumulation. Cognitive psychology, 57(3), 153-178.
"""
mutable struct LBA{T1,T2,T3,T4} <: SequentialSamplingModel
    ν::T1
    A::T2
    k::T3
    τ::T4
    σ::Float64
end

Base.broadcastable(x::LBA) = Ref(x)

LBA(;τ, A, k, ν, σ=1.0) = LBA(ν, A, k, τ, σ)

function select_winner(dt)
    if any(x -> x > 0, dt)
        mi,mv = 0,Inf
        for (i,t) in enumerate(dt)
            if (t > 0) && (t < mv)
                mi = i
                mv = t
            end
        end
    else
        return 1,-1.0
    end
    return mi,mv
end

function sample_drift_rates(ν, σ)
    negative = true
    v = similar(ν)
    while negative
        v = [rand(Normal(d,σ)) for d in ν]
        any(x -> x > 0, v) ? negative = false : nothing
    end
    return v
end

function rand(d::LBA)
    @unpack τ,A,k,ν,σ = d
    b = A + k
    N = length(ν)
    v = sample_drift_rates(ν, σ)
    a = rand(Uniform(0,A), N)
    dt = @. (b - a)/v
    choice,mn = select_winner(dt)
    rt = τ .+ mn
    return choice,rt
end

function rand(d::LBA, N::Int)
    choice = fill(0, N)
    rt = fill(0.0, N)
    for i in 1:N
        choice[i],rt[i] = rand(d)
    end
    return (choice=choice,rt=rt)
end

logpdf(d::LBA, choice, rt) = log(pdf(d, choice, rt))

function logpdf(d::LBA, data::T) where {T<:NamedTuple}
    return sum(logpdf.(d, data...))
end

function logpdf(dist::LBA, data::Array{<:Tuple,1})
    LL = 0.0
    for d in data
        LL += logpdf(dist, d...)
    end
    return LL
end

function pdf(d::LBA, c, rt)
    @unpack τ,A,k,ν,σ = d
    b = A + k; den = 1.0
    rt < τ ? (return 1e-10) : nothing
    for (i,v) in enumerate(ν)
        if c == i
            den *= dens(d, v, rt)
        else
            den *= (1-cummulative(d, v, rt))
        end
    end
    pneg = pnegative(d)
    den = den/(1-pneg)
    den = max(den,1e-10)
    isnan(den) ? (return 0.0) : (return den)
end

logpdf(d::LBA, data::Tuple) = logpdf(d, data...)

function dens(d::LBA, v, rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt - τ; b = A + k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    dens = (1/A)*(-v*cdf(Normal(0,1), n1) + σ*pdf(Normal(0,1), n1) +
        v*cdf(Normal(0,1), n2) - σ*pdf(Normal(0,1) ,n2))
    return dens
end

function cummulative(d::LBA, v, rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt - τ; b = A + k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    cm = 1 + ((b-A-dt*v)/A)*cdf(Normal(0,1), n1) -
        ((b-dt*v)/A)*cdf(Normal(0,1), n2) + ((dt*σ)/A)*pdf(Normal(0,1), n1) -
        ((dt*σ)/A)*pdf(Normal(0,1), n2)
    return cm
end

function pnegative(d::LBA)
    @unpack ν,σ=d
    p=1.0
    for v in ν
        p*= cdf(Normal(0,1), -v/σ)
    end
    return p
end


text:
function _source_f(l::AbstractVector, r::AbstractVector, ϵ::Real)

    L = dot(r, l)

    non_normal = r - L * l

    norm_non_normal = sqrt(norm(non_normal) ^ 2 + ϵ ^ 2)
    nr = sqrt(norm(r) ^ 2 + ϵ ^ 2)

    (l - non_normal * (L / norm_non_normal ^ 2)) / nr

end

"""
Source line influence coefficient
"""
function source_line(colpt::AbstractVector, r1::AbstractVector, r2::AbstractVector, ϵ::Real)

    l = r2 - r1
    nl = norm(l) + 1e-10
    @. l /= nl

    (_source_f(l, colpt - r2, ϵ) - _source_f(l, colpt - r1, ϵ)) / (4 * pi)

end

function _gradPhi(r⃗::AbstractVector, ξ⃗::AbstractVector, η⃗::AbstractVector, ϵ::Real)

    rxi = dot(r⃗, ξ⃗)
    reta = dot(r⃗, η⃗)

    r = sqrt(norm(r⃗) ^ 2 + ϵ ^ 2)

    den = (r ^ 2 - rxi ^ 2 + ϵ ^ 2) * r

    (
        (rxi * η⃗ + reta * ξ⃗) * den - (den * r⃗ / r ^ 2 + 2 * (r⃗ - rxi * ξ⃗) * r) * rxi * reta
    ) / den ^ 2

end

"""
Doublet line influence coefficient
"""
function doublet_line(colpt::AbstractVector, r1::AbstractVector, r2::AbstractVector, η⃗::AbstractVector, ϵ::Real)

    l = r2 - r1
    nl = norm(l) + 1e-10
    @. l /= nl

    (_gradPhi(colpt - r2, l, η⃗, ϵ) - _gradPhi(colpt - r1, l, η⃗, ϵ)) / (4 * pi)

end

"""
Point doublet influence coefficient
"""
function point_doublet_infl(colpt::AbstractVector, pt::AbstractVector, ξ::AbstractVector, ϵ::Real)

    r = colpt .- pt
    nr = sqrt(norm(r) ^ 2 + ϵ ^ 2)

    z = ξ ⋅ r
    nprojr = r .- z .* ξ

    (3 * z) .* nprojr ./ (nr ^ 5 * 4 * π) .- ξ .* ((nr ^ 2 - 3 * z ^ 2) / (4 * π * nr ^ 5))

end

"""
Point source influence coefficient
"""
function point_source_infl(colpt::AbstractVector, pt::AbstractVector, ϵ::Real)

    r = colpt .- pt

    r ./ (4 * π * sqrt(norm(r) + ϵ ^ 2) ^ 3)

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
# Autogenerated wrapper script for TinyXML_jll for aarch64-apple-darwin
export libtinyxml

JLLWrappers.@generate_wrapper_header("TinyXML")
JLLWrappers.@declare_library_product(libtinyxml, "@rpath/libtinyxml.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtinyxml,
        "lib/libtinyxml.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, LinearAlgebra
using QBase

@testset "./src/types/states.jl" begin

@testset "is_density_matrix()" begin
    @testset "invalid matrices" begin
        @test !is_density_matrix([1 0;0 1])
        @test !is_density_matrix([0.5 1;-1 0.5])
        @test !is_density_matrix([0.5 1;1 0.5])
        @test !is_density_matrix([1.5 0;0 -0.5])
        @test !is_density_matrix([0.5 im;-im 0.5])
        @test !is_density_matrix([1 4;4 1])
        @test !is_density_matrix([0.5 0;0 0.5; 0 0])
    end
    @testset "valid matrices" begin
        @test is_density_matrix([0.5 0;0 0.5]::Matrix{Float64})
        @test is_density_matrix([1 0;0 0]::Matrix{Int64})
        @test is_density_matrix([0 0;0 1])
        @test is_density_matrix([0.5 0.5;0.5 0.5])
        @test is_density_matrix([0.5 0.5im;-0.5im 0.5]::Matrix{Complex{Float64}})
        @test is_density_matrix(State([1 0;0 0]))
    end
    @testset "using atol" begin
        ϵ = 1e-6

        ρ_almost_tr1 = [0.5 0;0 0.5+1e-6]
        @test !is_density_matrix(ρ_almost_tr1)
        @test is_density_matrix(ρ_almost_tr1, atol=1e-5)

        ρ_almost_herm = [1+1e-6 0;0 -1e-6]
        @test !is_density_matrix(ρ_almost_herm)
        @test is_density_matrix(ρ_almost_herm, atol=1e-5)
    end
end

@testset "_density_matrix_error()" begin
    @testset "non-hermitian error" begin
        try
            QBase._density_matrix_error([0.5 0;0.5 0.5])
            @test false
        catch err
            @test err isa DomainError
            @test err.val == [0.5 0;0.5 0.5]
            @test err.msg == "Density matrix `ρ` is not hermitian."
        end
    end

    @testset "non-trace-one error" begin
        try
            QBase._density_matrix_error([0.5 0;0 0.6])
            @test false
        catch err
            @test err isa DomainError
            @test err.val == [0.5 0;0 0.6]
            @test err.msg == "Density matrix `ρ` is not trace-one."
        end
    end

    @testset "non-positive semi-definite error" begin
        try
            QBase._density_matrix_error([1.5 0;0 -0.5])
            @test false
        catch err
            @test err isa DomainError
            @test err.val == [1.5 0;0 -0.5]
            @test err.msg == "Density matrix `ρ` is not positive semi-definite."
        end
    end
end

@testset "State()" begin
    @testset "type inheritance" begin
        ρ_int = State([1 0;0 0])
        @test ρ_int isa State{Int64}
        @test ρ_int isa Operator
        @test ρ_int == [1 0;0 0]

        ρ_float = State([1. 0.;0. 0.])
        @test ρ_float isa State{Float64}
        @test ρ_float isa Operator
        @test ρ_float == [1 0;0 0]

        ρ_complex = State([0.5 0.5im;-0.5im 0.5])
        @test ρ_complex isa State{Complex{Float64}}
        @test ρ_complex isa Operator
        @test ρ_complex == [0.5 0.5im;-0.5im 0.5]

        ρ_rational = State([1//2 1//2; 1//2 1//2])
        @test ρ_rational isa State{Rational{Int64}}
        @test ρ_rational isa Operator
        @test ρ_rational == [0.5 0.5;0.5 0.5]
    end

    @testset "adjoint matrix input" begin
        ρ = State([0.5 -0.5im;0.5im 0.5]')
        @test ρ isa State{Complex{Float64}}
        @test ρ == [0.5 -0.5im;0.5im 0.5]
    end

    @testset "using atol" begin
        ϵ = 1e-5

        @test_throws DomainError State([1 0;0 ϵ])
        @test State([1 0;0 ϵ], atol=1e-4) isa State{Float64}
    end
end

@testset "kron() for states" begin
    ρ0_int = State([1 0;0 0])
    ρ0_float = State([1. 0;0 0])

    ρ1_int = State([0 0;0 1])
    ρ1_float = State([0 0;0 1.])

    ρ01 = kron(ρ0_int,ρ1_int)
    @test ρ01 isa State{Int64}
    @test ρ01 == [0 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0]

    ρ11 = kron(ρ1_int, ρ1_float)
    @test ρ11 isa State{Float64}
    @test ρ11 == [0 0 0 0;0 0 0 0;0 0 0 0;0 0 0 1]

    ρ100 = kron(ρ1_float, ρ0_float, ρ0_int)
    @test ρ100 isa State{Float64}
    @test ρ100 ==  [
        0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0;
        0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0;
    ]

    @testset "using atol" begin
        ϵ = 1e-7

        ρ = State([0.5 0.5;0.5 0.5 + ϵ])

        @test_throws DomainError kron(ρ,ρ)

        ρρ = kron(ρ, ρ, atol=1e-6)
        @test ρρ isa State{Float64}
        @test isapprox(ρρ, ones(Float64,4,4)/4, atol=1e-6)
        @test ρρ.atol == 1e-6
    end
end

@testset "partial_trace() for states" begin
    @testset "maximally entangled bell state" begin
        ρ_AB = State([1 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 1]/2)

        ρ_A = partial_trace(ρ_AB, [2,2], 2)
        @test ρ_A isa State{Float64}
        @test ρ_A == [1 0;0 1]/2

        ρ_B = partial_trace(ρ_AB, [2,2], 1)
        @test ρ_B isa State{Float64}
        @test ρ_B == [1 0;0 1]/2
    end

    @testset "separable qubit states" begin
        ρ_x0 = State([1 1;1 1]/2)
        ρ_0 = State([1 0;0 0])

        ρ_AB = kron(ρ_x0, ρ_0)
        @test ρ_AB isa State{Float64}

        ρ_A = partial_trace(ρ_AB, [2,2], 2)
        @test ρ_A isa State{Float64}
        @test ρ_A == ρ_x0

        ρ_B = partial_trace(ρ_AB, [2,2], 1)
        @test ρ_B isa State{Float64}
        @test ρ_B == ρ_0
    end

    @testset "using atol" begin
        ρ_AB = State([1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1], atol= 3)

        @test_throws DomainError partial_trace(ρ_AB, [2,2], 2)

        ρ_A = partial_trace(ρ_AB, [2,2], 2, atol=3)
        @test ρ_A isa State{Int64}
        @test ρ_A == [2 2;2 2]
        @test ρ_A.atol == 3
    end
end

@testset "rank() for State types" begin
    @test 1 == rank(pure_state([1,-im]/sqrt(2)))
    @test 2 == rank(State([1 0;0 1]/2))
    @test 1 == rank(State([0.9 0;0 0.1], atol=0.11))
    @test 2 == rank(State([0.9 0;0 0.1]))
end

@testset "eigvals() for State types" begin
    priors = Probabilities([1/3,1/3,1/3])

    qb1 = bloch_qubit_state(2*π/3,0)
    qb2 = bloch_qubit_state(5*π/12,π)
    qb3 = bloch_qubit_state(0,0)
    ρ_states = [qb1,qb2,qb3]

    ρ_mix = mixed_state(priors, ρ_states)
    @test isapprox(eigvals(ρ_mix), [0.372, 0.628], atol=1e-3)
end

@testset "is_pure()" begin
    @test is_pure(State([1 0;0 0]))
    @test is_pure([0 0;0 1])

    @test !is_pure(State([0.5 0;0 0.5]))
    @test_throws DomainError is_pure(State([1 0;0 1]))

    @test !is_pure(State([1-1e-5 0;0 1e-5]))
    @test is_pure(State([1-1e-5 0;0 1e-5], atol=1e-4))

end

@testset "is_mixed()" begin
    @test !is_mixed(State([1 0;0 0]))
    @test !is_mixed([0 0;0 1])

    @test is_mixed(State([0.5 0;0 0.5]))
    @test_throws DomainError is_mixed(State([1 0;0 1]))

    @test is_mixed(State([1-1e-5 0;0 1e-5]))
    @test !is_mixed(State([1-1e-5 0;0 1e-5], atol=1e-4))
end

end


text:
include("common.jl")
@testset "HistoryFunction" begin
    # check constant extrapolation with problem with vanishing delays at t = 0
    @testset "vanishing delays" begin
        prob = DDEProblem((du,u,h,p,t) -> -h(p, t/2)[1], [1.0], (p, t) -> [1.0], (0.0, 10.0))
        solve(prob, MethodOfSteps(RK4()))
    end

    @testset "general" begin
        # naive history functions
        h_notinplace(p, t; idxs=nothing) = idxs === nothing ? [t, -t] : [t, -t][idxs]

        function h_inplace(val, p, t; idxs=nothing)
            if idxs === nothing
                val[1] = t
                val[2] = -t
            else
                val .= [t; -t][idxs]
            end
        end

        @testset "agrees (h=$h)" for h in (h_notinplace, h_inplace)
            @test DelayDiffEq.agrees(h, zeros(2), nothing, 0)
            @test !DelayDiffEq.agrees(h, ones(2), nothing, 1)
        end

        # ODE integrator
        prob = ODEProblem((du,u,p,t)->@.(du=p*u), ones(2), (0.0, 1.0),1.01)
        integrator = init(prob, Tsit5())

        # combined history function
        history_notinplace = DelayDiffEq.HistoryFunction(h_notinplace,
                                                         integrator.sol,
                                                         integrator)
        history_inplace = DelayDiffEq.HistoryFunction(h_inplace,
                                                      integrator.sol,
                                                      integrator)

        # test evaluation of history function
        @testset "evaluation" for idxs in (nothing, [2]) # (idxs=$idxs)
            # expected value
            trueval = h_notinplace(nothing, -1; idxs = idxs)

            # out-of-place
            @test history_notinplace(nothing, -1, Val{0}; idxs = idxs) == trueval

            # in-place
            val = zero(trueval)
            history_inplace(val, nothing, -1; idxs = idxs)
            @test val == trueval

            val = zero(trueval)
            history_inplace(val, nothing, -1, Val{0}; idxs = idxs)
            @test val == trueval
        end

        # test constant extrapolation
        @testset "constant extrapolation" for #  (deriv=$deriv, idxs=$idxs)
            deriv in (Val{0}, Val{1}), idxs in (nothing, [2])
            # expected value
            trueval = deriv == Val{0} ?
                (idxs === nothing ? integrator.u : integrator.u[[2]]) :
                (idxs === nothing ? zeros(2) : [0.0])

            # out-of-place
            integrator.isout = false
            @test history_notinplace(nothing, 1, deriv; idxs = idxs) == trueval &&
                integrator.isout

            # in-place
            integrator.isout = false
            @test history_inplace(nothing, nothing, 1, deriv; idxs = idxs) == trueval &&
                integrator.isout

            integrator.isout = false
            val = 1 .- trueval # ensures that val ≠ trueval
            history_inplace(val, nothing, 1, deriv; idxs = idxs)
            @test val == trueval && integrator.isout
        end

        # add step to integrator
        @testset "update integrator" begin
            OrdinaryDiffEq.loopheader!(integrator)
            OrdinaryDiffEq.perform_step!(integrator, integrator.cache)
            integrator.t = integrator.dt
            @test 0.01 < integrator.t < 1
            @test integrator.sol.t[end] == 0
        end

        # test integrator interpolation
        @testset "integrator interpolation" for # (deriv=$deriv, idxs=$idxs)
            deriv in (Val{0}, Val{1}), idxs in (nothing, [2])
            # expected value
            trueval = OrdinaryDiffEq.current_interpolant(0.01, integrator, idxs, deriv)

            # out-of-place
            integrator.isout = false
            @test history_notinplace(nothing, 0.01, deriv; idxs = idxs) == trueval &&
                integrator.isout

            # in-place
            integrator.isout = false
            val = zero(trueval)
            history_inplace(val, nothing, 0.01, deriv; idxs = idxs)
            @test val == trueval && integrator.isout
        end

        # add step to solution
        @testset "update solution" begin
            integrator.t = 0
            OrdinaryDiffEq.loopfooter!(integrator)
            @test integrator.t == integrator.sol.t[end]
        end

        # test solution interpolation
        @testset "solution interpolation" for # (deriv=$deriv, idxs=$idxs)
            deriv in (Val{0}, Val{1}), idxs in (nothing, [2])
            # expected value
            trueval = integrator.sol.interp(0.01, idxs, deriv, integrator.p)

            # out-of-place
            @test history_notinplace(nothing, 0.01, deriv; idxs = idxs) == trueval &&
                !integrator.isout

            # in-place
            val = zero(trueval)
            history_inplace(val, nothing, 0.01, deriv; idxs = idxs)
            @test val == trueval && !integrator.isout
        end

        # test integrator extrapolation
        @testset "integrator extrapolation" for # (deriv=$deriv, idxs=$idxs)
            deriv in (Val{0}, Val{1}), idxs in (0, [2])
            idxs == 0 && (idxs = nothing)
            # expected value
            trueval = OrdinaryDiffEq.current_interpolant(1, integrator, idxs, deriv)

            # out-of-place
            integrator.isout = false
            @test history_notinplace(nothing, 1, deriv; idxs = idxs) == trueval &&
                integrator.isout

            # in-place
            integrator.isout = false
            val = zero(trueval)
            history_inplace(val, nothing, 1, deriv; idxs = idxs)
            @test val == trueval && integrator.isout
        end
    end
end


text:
using Oceananigans.AbstractOperations: UnaryOperation, Derivative, BinaryOperation, MultiaryOperation
using Oceananigans.Fields: PressureField
using Oceananigans.Buoyancy: BuoyancyField

function simple_binary_operation(op, a, b, num1, num2)
    a_b = op(a, b)
    interior(a) .= num1
    interior(b) .= num2
    return a_b[2, 2, 2] == op(num1, num2)
end

function three_field_addition(a, b, c, num1, num2)
    a_b_c = a + b + c
    interior(a) .= num1
    interior(b) .= num2
    interior(c) .= num2
    return a_b_c[2, 2, 2] == num1 + num2 + num2
end

function x_derivative(a)
    dx_a = ∂x(a)

    for k in 1:3
        interior(a)[:, 1, k] .= [1, 2, 3]
        interior(a)[:, 2, k] .= [1, 2, 3]
        interior(a)[:, 3, k] .= [1, 2, 3]
    end

    return dx_a[2, 2, 2] == 1
end

function y_derivative(a)
    dy_a = ∂y(a)

    for k in 1:3
        interior(a)[1, :, k] .= [1, 3, 5]
        interior(a)[2, :, k] .= [1, 3, 5]
        interior(a)[3, :, k] .= [1, 3, 5]
    end

    return dy_a[2, 2, 2] == 2
end

function z_derivative(a)
    dz_a = ∂z(a)

    for k in 1:3
        interior(a)[1, k, :] .= [1, 4, 7]
        interior(a)[2, k, :] .= [1, 4, 7]
        interior(a)[3, k, :] .= [1, 4, 7]
    end

    return dz_a[2, 2, 2] == 3
end

function x_derivative_cell(FT, arch)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3, 3, 3))
    a = Field(Cell, Cell, Cell, arch, grid, nothing)
    dx_a = ∂x(a)

    for k in 1:3
        interior(a)[:, 1, k] .= [1, 4, 4]
        interior(a)[:, 2, k] .= [1, 4, 4]
        interior(a)[:, 3, k] .= [1, 4, 4]
    end

    return dx_a[2, 2, 2] == 3
end

function times_x_derivative(a, b, location, i, j, k, answer)
    a∇b = @at location b * ∂x(a)
    return a∇b[i, j, k] == answer
end

function compute_derivative(model, ∂)
    T, S = model.tracers
    S.data.parent .= π

    @compute ∂S = ComputedField(∂(S))
    result = Array(interior(∂S))

    return all(result .≈ zero(eltype(model.grid)))
end

function compute_unary(unary, model)
    set!(model; S=π)
    T, S = model.tracers

    @compute uS = ComputedField(unary(S), data=model.pressures.pHY′.data)
    result = Array(interior(uS))

    return all(result .≈ unary(eltype(model.grid)(π)))
end

function compute_plus(model)
    set!(model; S=π, T=42)
    T, S = model.tracers

    @compute ST = ComputedField(S + T, data=model.pressures.pHY′.data)

    result = Array(interior(ST))

    return all(result .≈ eltype(model.grid)(π+42))
end

function compute_many_plus(model)
    set!(model; u=2, S=π, T=42)
    T, S = model.tracers
    u, v, w = model.velocities

    @compute uTS = ComputedField(@at((Cell, Cell, Cell), u + T + S))
    result = Array(interior(uTS))

    return all(result .≈ eltype(model.grid)(2+π+42))
end

function compute_minus(model)
    set!(model; S=π, T=42)
    T, S = model.tracers

    @compute ST = ComputedField(S - T, data=model.pressures.pHY′.data)
    result = Array(interior(ST))

    return all(result .≈ eltype(model.grid)(π-42))
end

function compute_times(model)
    set!(model; S=π, T=42)
    T, S = model.tracers

    @compute ST = ComputedField(S * T, data=model.pressures.pHY′.data)
    result = Array(interior(ST))

    return all(result .≈ eltype(model.grid)(π*42))
end

function compute_kinetic_energy(model)
    u, v, w = model.velocities
    set!(u, 1)
    set!(v, 2)
    set!(w, 3)

    kinetic_energy_operation = @at (Cell, Cell, Cell) (u^2 + v^2 + w^2) / 2
    @compute kinetic_energy = ComputedField(kinetic_energy_operation, data=model.pressures.pHY′.data)
    result = Array(interior(kinetic_energy))

    return all(result .≈ 7)
end

function horizontal_average_of_plus(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    S₀(x, y, z) = sin(π * z)
    T₀(x, y, z) = 42 * z
    set!(model; S=S₀, T=T₀)
    T, S = model.tracers

    @compute ST = AveragedField(S + T, dims=(1, 2))

    zC = znodes(Cell, model.grid)
    correct_profile = @. sin(π * zC) + 42 * zC

    result = Array(interior(ST))[:]

    return all(result .≈ correct_profile)
end

function zonal_average_of_plus(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    S₀(x, y, z) = sin(π*z) * sin(π*y)
    T₀(x, y, z) = 42*z + y^2
    set!(model; S=S₀, T=T₀)
    T, S = model.tracers

    @compute ST = AveragedField(S + T, dims=1)

    yC = ynodes(Cell, model.grid, reshape=true)
    zC = znodes(Cell, model.grid, reshape=true)
    correct_slice = @. sin(π * zC) * sin(π * yC) + 42*zC + yC^2

    result = Array(interior(ST))

    return all(result[1, :, :] .≈ correct_slice[1, :, :])
end

function volume_average_of_times(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    S₀(x, y, z) = 1 + sin(2π*x)
    T₀(x, y, z) = y
    set!(model; S=S₀, T=T₀)
    T, S = model.tracers

    @compute ST = AveragedField(S * T, dims=(1, 2, 3))

    result = Array(interior(ST))

    return all(result[1, 1, 1] .≈ 0.5)
end

function horizontal_average_of_minus(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    S₀(x, y, z) = sin(π*z)
    T₀(x, y, z) = 42*z
    set!(model; S=S₀, T=T₀)
    T, S = model.tracers

    @compute ST = AveragedField(S - T, dims=(1, 2))

    zC = znodes(Cell, model.grid)
    correct_profile = @. sin(π * zC) - 42 * zC

    result = Array(interior(ST))

    return all(result[1, 1, 1:Nz] .≈ correct_profile)
end

function horizontal_average_of_times(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    S₀(x, y, z) = sin(π*z)
    T₀(x, y, z) = 42*z
    set!(model; S=S₀, T=T₀)
    T, S = model.tracers

    @compute ST = AveragedField(S * T, dims=(1, 2))

    zC = znodes(Cell, model.grid)
    correct_profile = @. sin(π * zC) * 42 * zC

    result = Array(interior(ST))

    return all(result[1, 1, 1:Nz] .≈ correct_profile)
end

function multiplication_and_derivative_ccf(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    w₀(x, y, z) = sin(π * z)
    T₀(x, y, z) = 42 * z
    set!(model; w=w₀, T=T₀)

    w = model.velocities.w
    T = model.tracers.T

    @compute wT = AveragedField(w * ∂z(T), dims=(1, 2))

    zF = znodes(Face, model.grid)
    correct_profile = @. 42 * sin(π * zF)

    result = Array(interior(wT))

    # Omit both halos and boundary points
    return all(result[1, 1, 2:Nz] .≈ correct_profile[2:Nz])
end

const C = Cell
const F = Face

function multiplication_and_derivative_ccc(model)
    Ny, Nz = model.grid.Ny, model.grid.Nz

    w₀(x, y, z) = sin(π*z)
    T₀(x, y, z) = 42*z
    set!(model; w=w₀, T=T₀)

    w = model.velocities.w
    T = model.tracers.T

    wT_ccc = @at (C, C, C) w * ∂z(T)
    @compute wT_ccc_avg = AveragedField(wT_ccc, dims=(1, 2))

    zF = znodes(Face, model.grid)
    sinusoid = sin.(π * zF)
    interped_sin = [(sinusoid[k] + sinusoid[k+1]) / 2 for k in 1:model.grid.Nz]
    correct_profile = interped_sin .* 42

    result = Array(interior(wT_ccc_avg))

    # Omit boundary-adjacent points from comparison
    return all(result[1, 1, 2:Nz-1] .≈ correct_profile[2:Nz-1])
end

function computation_including_boundaries(FT, arch)
    topo = (Periodic, Bounded, Bounded)
    grid = RegularCartesianGrid(FT, topology=topo, size=(13, 17, 19), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, float_type=FT, grid=grid)

    u, v, w = model.velocities
    @. u.data = 1 + rand()
    @. v.data = 2 + rand()
    @. w.data = 3 + rand()

    op = @at (Cell, Face, Face) u * v * w
    @compute uvw = ComputedField(op)

    return all(interior(uvw) .!= 0)
end

function operations_with_computed_field(model)
    u, v, w = model.velocities
    uv = ComputedField(u * v)
    @compute uvw = ComputedField(uv * w)
    return true
end

function operations_with_averaged_field(model)
    u, v, w = model.velocities
    UV = AveragedField(u * v, dims=(1, 2))
    wUV = ComputedField(w * UV)
    compute!(wUV)
    return true
end

function pressure_field(model)
    p = PressureField(model)
    u, v, w = model.velocities
    @compute up = ComputedField(u * p)
    return true
end

function computations_with_buoyancy_field(FT, arch, buoyancy)
    grid = RegularCartesianGrid(FT, size=(1, 1, 1), extent=(1, 1, 1))
    tracers = buoyancy isa BuoyancyTracer ? :b : (:T, :S)
    model = IncompressibleModel(architecture=arch, float_type=FT, grid=grid,
                                tracers=tracers, buoyancy=buoyancy)

    b = BuoyancyField(model)
    u, v, w = model.velocities

    compute!(b)

    ub = ComputedField(b * u)
    vb = ComputedField(b * v)
    wb = ComputedField(b * w)

    compute!(ub)
    compute!(vb)
    compute!(wb)

    return true # test that it doesn't error
end

function computations_with_averaged_fields(model)
    u, v, w, T, S = fields(model)

    set!(model, u = (x, y, z) -> z, v = 2, w = 3)

    # Two ways to compute turbulent kinetic energy
    U = AveragedField(u, dims=(1, 2))
    V = AveragedField(v, dims=(1, 2))

    tke_op = @at (Cell, Cell, Cell) ((u - U)^2  + (v - V)^2 + w^2) / 2
    tke = ComputedField(tke_op)
    compute!(tke)

    return all(interior(tke)[2:3, 2:3, 2:3] .== 9/2)
end

function computations_with_averaged_field_derivative(model)

    set!(model, u = (x, y, z) -> z, v = 2, w = 3)

    u, v, w, T, S = fields(model)

    # Two ways to compute turbulent kinetic energy
    U = AveragedField(u, dims=(1, 2))
    V = AveragedField(v, dims=(1, 2))

    # This tests a vertical derivative of an AveragedField
    shear_production_op = @at (Cell, Cell, Cell) u * w * ∂z(U)
    shear = ComputedField(shear_production_op)
    compute!(shear)

    set!(model, T = (x, y, z) -> 3 * z)

    return all(interior(shear)[2:3, 2:3, 2:3] .== interior(T)[2:3, 2:3, 2:3])
end

function computations_with_computed_fields(model)
    u, v, w, T, S = fields(model)

    set!(model, u = (x, y, z) -> z, v = 2, w = 3)

    # Two ways to compute turbulent kinetic energy
    U = AveragedField(u, dims=(1, 2))
    V = AveragedField(v, dims=(1, 2))

    u′ = ComputedField(u - U)
    v′ = ComputedField(v - V)

    tke_op = @at (Cell, Cell, Cell) (u′^2  + v′^2 + w^2) / 2
    tke = ComputedField(tke_op)
    compute!(tke)

    return all(interior(tke)[2:3, 2:3, 2:3] .== 9/2)
end

@testset "Abstract operations" begin
    @info "Testing abstract operations..."

    for FT in float_types
        arch = CPU()
        grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3, 3, 3))
        u, v, w = VelocityFields(arch, grid)
        c = Field(Cell, Cell, Cell, arch, grid, nothing)

        @testset "Unary operations and derivatives [$FT]" begin
            for ψ in (u, v, w, c)
                for op_symbol in Oceananigans.AbstractOperations.unary_operators
                    op = eval(op_symbol)
                    @test typeof(op(ψ)[2, 2, 2]) <: Number
                end

                for d_symbol in Oceananigans.AbstractOperations.derivative_operators
                    d = eval(d_symbol)
                    @test typeof(d(ψ)[2, 2, 2]) <: Number
                end
            end
        end

        @testset "Binary operations [$FT]" begin
            generic_function(x, y, z) = x + y + z
            for (ψ, ϕ) in ((u, v), (u, w), (v, w), (u, c), (generic_function, c), (u, generic_function))
                for op_symbol in Oceananigans.AbstractOperations.binary_operators
                    op = eval(op_symbol)
                    @test typeof(op(ψ, ϕ)[2, 2, 2]) <: Number
                end
            end
        end

        @testset "Multiary operations [$FT]" begin
            generic_function(x, y, z) = x + y + z
            for (ψ, ϕ, σ) in ((u, v, w), (u, v, c), (u, v, generic_function))
                for op_symbol in Oceananigans.AbstractOperations.multiary_operators
                    op = eval(op_symbol)
                    @test typeof(op((Cell, Cell, Cell), ψ, ϕ, σ)[2, 2, 2]) <: Number
                end
            end
        end
    end

    @testset "Fidelity of simple binary operations" begin
        arch = CPU()
        @info "  Testing simple binary operations..."
        for FT in float_types
            num1 = FT(π)
            num2 = FT(42)
            grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3, 3, 3))

            u, v, w = VelocityFields(arch, grid)
            T, S = TracerFields((:T, :S), arch, grid)

            for op in (+, *, -, /)
                @test simple_binary_operation(op, u, v, num1, num2)
                @test simple_binary_operation(op, u, w, num1, num2)
                @test simple_binary_operation(op, u, T, num1, num2)
                @test simple_binary_operation(op, T, S, num1, num2)
            end
            @test three_field_addition(u, v, w, num1, num2)
        end
    end

    @testset "Derivatives" begin
        arch = CPU()
        @info "  Testing derivatives..."
        for FT in float_types
            grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(3, 3, 3),
                                        topology=(Periodic, Periodic, Periodic))

            u, v, w = VelocityFields(arch, grid)
            T, S = TracerFields((:T, :S), arch, grid)
            for a in (u, v, w, T)
                @test x_derivative(a)
                @test y_derivative(a)
                @test z_derivative(a)
            end
            @test x_derivative_cell(FT, arch)
        end
    end

    @testset "Combined binary operations and derivatives" begin
        @info "  Testing combined binary operations and derivatives..."
        arch = CPU()
        Nx = 3 # Δx=1, xC = 0.5, 1.5, 2.5
        for FT in float_types
            grid = RegularCartesianGrid(FT, size=(Nx, Nx, Nx), extent=(Nx, Nx, Nx))
            a, b = (Field(Cell, Cell, Cell, arch, grid, nothing) for i in 1:2)

            set!(b, 2)
            set!(a, (x, y, z) -> x < 2 ? 3x : 6)

            #                            0   0.5   1   1.5   2   2.5   3
            # x -▶                  ∘ ~~~|--- * ---|--- * ---|--- * ---|~~~ ∘
            #        i Face:    0        1         2        3          4
            #        i Cell:        0         1         2         3         4

            #              a = [    0,       1.5,      4.5,       6,        0    ]
            #              b = [    0,        2,        2,        2,        0    ]
            #          ∂x(a) = [        1.5,       3,       1.5,      -6         ]

            # x -▶                  ∘ ~~~|--- * ---|--- * ---|--- * ---|~~~ ∘
            #        i Face:    0        1         2         3         4
            #        i Cell:        0         1         2         3         4

            # ccc: b * ∂x(a) = [             4.5,      4.5      -4.5,            ]
            # fcc: b * ∂x(a) = [         3,        6,        3,       -6         ]

            @test times_x_derivative(a, b, (C, C, C), 1, 2, 2, 4.5)
            @test times_x_derivative(a, b, (F, C, C), 1, 2, 2, 3)

            @test times_x_derivative(a, b, (C, C, C), 2, 2, 2, 4.5)
            @test times_x_derivative(a, b, (F, C, C), 2, 2, 2, 6)

            @test times_x_derivative(a, b, (C, C, C), 3, 2, 2, -4.5)
            @test times_x_derivative(a, b, (F, C, C), 3, 2, 2, 3)
        end
    end

    for arch in archs
        @testset "AbstractOperations and ComputedFields [$(typeof(arch))]" begin

            @info "  Testing combined binary operations and derivatives..."

            for FT in (Float64,) #float_types

                grid = RegularCartesianGrid(FT, size=(4, 4, 4), extent=(1, 1, 1),
                                            topology=(Periodic, Periodic, Bounded))

                buoyancy = SeawaterBuoyancy(gravitational_acceleration = 1,
                                                     equation_of_state = LinearEquationOfState(α=1, β=1))

                model = IncompressibleModel(architecture = arch,
                                              float_type = FT,
                                                    grid = grid,
                                                buoyancy = buoyancy)

                @testset "Construction of abstract operations [$FT, $(typeof(arch))]" begin
                    @info "    Testing construction of abstract operations [$FT, $(typeof(arch))]..."

                    u, v, w, T, S = fields(model)

                    @test_throws ArgumentError @at (Nothing, Nothing, Cell) T * S

                    for ϕ in (u, v, w, T, S)
                        for op in (sin, cos, sqrt, exp, tanh)
                            @test op(ϕ) isa UnaryOperation
                        end

                        for ∂ in (∂x, ∂y, ∂z)
                            @test ∂(ϕ) isa Derivative
                        end

                        @test u ^ 2 isa BinaryOperation
                        @test u * 2 isa BinaryOperation
                        @test u + 2 isa BinaryOperation
                        @test u - 2 isa BinaryOperation
                        @test u / 2 isa BinaryOperation

                        for ψ in (u, v, w, T, S)
                            @test ψ ^ ϕ isa BinaryOperation
                            @test ψ * ϕ isa BinaryOperation
                            @test ψ + ϕ isa BinaryOperation
                            @test ψ - ϕ isa BinaryOperation
                            @test ψ / ϕ isa BinaryOperation

                            for χ in (u, v, w, T, S)
                                @test ψ * ϕ * χ isa MultiaryOperation
                                @test ψ + ϕ + χ isa MultiaryOperation
                            end
                        end
                    end
                end

                @info "    Testing computation of abstract operations [$FT, $(typeof(arch))]..."

                @testset "Derivative computations [$FT, $(typeof(arch))]" begin
                    @info "      Testing compute! derivatives..."
                    @test compute_derivative(model, ∂x)
                    @test compute_derivative(model, ∂y)
                    @test compute_derivative(model, ∂z)
                end

                @testset "Unary computations [$FT, $(typeof(arch))]" begin
                    @info "      Testing compute! unary operations..."
                    for unary in Oceananigans.AbstractOperations.unary_operators
                        @test compute_unary(eval(unary), model)
                    end
                end

                @testset "Binary computations [$FT, $(typeof(arch))]" begin
                    @info "      Testing compute! binary operations..."
                    @test compute_plus(model)
                    @test compute_minus(model)
                    @test compute_times(model)
                end

                @testset "Multiary computations [$FT, $(typeof(arch))]" begin
                    @info "      Testing compute! multiary operations..."
                    @test compute_many_plus(model)

                    @info "      Testing compute! kinetic energy..."
                    @test compute_kinetic_energy(model)
                end

                @testset "Operations with ComputedField and PressureField [$FT, $(typeof(arch))]" begin
                    @info "      Testing operations with ComputedField..."
                    @test operations_with_computed_field(model)

                    @info "      Testing PressureField..."
                    @test pressure_field(model)
                end

                @testset "Horizontal averages of operations [$FT, $(typeof(arch))]" begin
                    @info "      Testing horizontal averges..."
                    @test horizontal_average_of_plus(model)
                    @test horizontal_average_of_minus(model)
                    @test horizontal_average_of_times(model)

                    @test multiplication_and_derivative_ccf(model)
                    @test multiplication_and_derivative_ccc(model)
                end

                @testset "Zonal averages of operations [$FT, $(typeof(arch))]" begin
                    @info "      Testing zonal averges..."
                    @test zonal_average_of_plus(model)
                end

                @testset "Volume averages of operations [$FT, $(typeof(arch))]" begin
                    @info "      Testing volume averges..."
                    @test volume_average_of_times(model)
                end

                @testset "Operations with AveragedField [$FT, $(typeof(arch))]" begin
                    @info "      Testing operations with AveragedField..."

                    T, S = model.tracers

                    TS = AveragedField(T * S, dims=(1, 2))

                    @test_throws ArgumentError @at (Nothing, Nothing, Cell) T * S
                    @test_throws ArgumentError TS * S

                    @test operations_with_averaged_field(model)
                end

                @testset "Compute! on faces along bounded dimensions" begin
                    @info "      Testing compute! on faces along bounded dimensions..."
                    @test computation_including_boundaries(FT, arch)
                end

                EquationsOfState = (LinearEquationOfState, SeawaterPolynomials.RoquetEquationOfState,
                                    SeawaterPolynomials.TEOS10EquationOfState)

                buoyancies = (BuoyancyTracer(), SeawaterBuoyancy(FT),
                              (SeawaterBuoyancy(FT, equation_of_state=eos(FT)) for eos in EquationsOfState)...)

                for buoyancy in buoyancies
                    @testset "Computations with BuoyancyFields [$FT, $(typeof(arch)), $(typeof(buoyancy).name.wrapper)]" begin
                        @info "      Testing computations with BuoyancyField " *
                              "[$FT, $(typeof(arch)), $(typeof(buoyancy).name.wrapper)]..."

                        @test computations_with_buoyancy_field(FT, arch, buoyancy)
                    end
                end

                @testset "Computations with AveragedFields [$FT, $(typeof(arch))]" begin
                    @info "      Testing computations with AveragedField [$FT, $(typeof(arch))]..."

                    # These don't work on the GPU right now
                    if arch isa CPU
                        @test computations_with_averaged_fields(model)
                        @test computations_with_averaged_field_derivative(model)
                    else
                        @test_skip computations_with_averaged_fields(model)
                        @test_skip computations_with_averaged_field_derivative(model)
                    end
                end

                @testset "Computations with ComputedFields [$FT, $(typeof(arch))]" begin
                    @info "      Testing computations with ComputedField [$FT, $(typeof(arch))]..."

                    # These don't work on the GPU right now
                    if arch isa CPU
                        @test computations_with_computed_fields(model)
                    else
                        @test_skip computations_with_computed_fields(model)
                    end
                end

                @testset "Conditional computation of ComputedField and BuoyancyField [$FT, $(typeof(arch))]" begin
                    @info "      Testing conditional computation of ComputedField and BuoyancyField " *
                          "[$FT, $(typeof(arch))]..."

                    set!(model, u=2, v=0, w=0, T=3, S=0)
                    u, v, w, T, S = fields(model)

                    uT = ComputedField(u * T)

                    α = model.buoyancy.equation_of_state.α
                    g = model.buoyancy.gravitational_acceleration
                    b = BuoyancyField(model)

                    compute!(uT, 1.0)
                    compute!(b, 1.0)
                    @test all(interior(uT) .== 6)
                    @test all(interior(b) .== g * α * 3)

                    set!(model, u=2, T=4)
                    compute!(uT, 1.0)
                    compute!(b, 1.0)
                    @test all(interior(uT) .== 6)
                    @test all(interior(b) .== g * α * 3)

                    compute!(uT, 2.0)
                    compute!(b, 2.0)
                    @test all(interior(uT) .== 8)
                    @test all(interior(b) .== g * α * 4)
                end
            end
        end
    end
end


text:
using LinearAlgebra, ForwardDiff

function row_col!(row,col,r,c)
    for cc in c
        for rr in r
            push!(row,convert(Int,rr))
            push!(col,convert(Int,cc))
        end
    end
    return row, col
end

v = [10.0; 1.0; 0.]
n = 0.0
μ = 0.5

P = [1. 0. 0.;
	 0. 1. 0.]

_norm(x) = sqrt(x'*x)
function vec_norm(x)
	if norm(x) == 0.0
		return ones(length(x))./norm(ones(length(x)))
	else
		x./_norm(x)
	end
end

function d_vec_norm(x)
	if norm(x) == 0.0
		y = 1.0*ones(length(x))#./norm(ones(length(x)))
		return (I - y*y'/(y'*y))/_norm(y)
	else
		(I - x*x'/(x'*x))/_norm(x)
	end
end
b = 1.0e-3*ones(2)
b = zeros(2)
d_vec_norm(b)
_norm(b)
vec_norm(b)
d_vec_norm(b)
norm(vec(d_vec_norm(b)) - vec(ForwardDiff.jacobian(vec_norm,norm(b) == 0.0 ? 1.0*ones(2) : b)))

ψ = 0.1
P*v + ψ*vec_norm(b)

##
include(joinpath(pwd(),"src/Ipopt.jl"))

nx = 4
nc = 2 + 1 + 1
function obj(x)
	b = x[1:2]
	ψ = x[3]
	s = x[4]

	return 100.0*s
end

function con!(c,x)
	b = x[1:2]
	ψ = x[3]
	s = x[4]

	c[1:2] = P*v + ψ*vec_norm(b)
	c[3] = μ*n - _norm(b)
	c[4] = s - min(ψ,(μ*n - _norm(b)))

	return nothing
end

function ∇con!(C,x)
	b = x[1:2]
	ψ = x[3]

	C[1:2,1:2] = ψ*d_vec_norm(b)
	C[1:2,3] = vec_norm(b)

	C[3,1:2] = -1.0*vec_norm(b)
	C[3,3] = 0.0

	if ψ < μ*n - _norm(b)
		C[4,3] = -1.0
	else
		C[4,1:2] = vec_norm(b)
	end
	C[4,4] = 1.0
	return nothing
end

_x = [1.0e-3; 1.0e-3; 1.0e-3; 1.0e-3]
obj(_x)
_c = zeros(nc)
con!(_c,_x)
_c
_C = zeros(nc,nx)
∇con!(_C,_x)

function eval_f(x::Vector{Float64})
  return obj(x)
end

function eval_g(x::Vector{Float64}, g::Vector{Float64})
	con!(g,x)

	return g
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
	grad_f .= ForwardDiff.gradient(obj,x)
end

function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
	if mode == :Structure
		row = []
		col = []
		row_col!(row,col,(1:nc),(1:nx))
		rows .= vec(row)
		cols .= vec(col)
	else
		∇con!(reshape(values,nc,nx),x)
	end
end

function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
	if mode == :Structure
		nothing
	else
		nothing
	end
end

x_L = [-1.0e8; -1.0e8; 0.0; 0.0]
x_U = 1.0e8*ones(nx)

g_L = zeros(nc)
g_U = [zeros(2);1.0e8;1.0e8]

prob = Ipopt.createProblem(nx, x_L, x_U, nc, g_L, g_U, nc*nx, nx*nx,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
Ipopt.addOption(prob,"hessian_approximation","limited-memory")
# Ipopt.addOption(prob,"mu_strategy","monotone")

Ipopt.addOption(prob,"constr_viol_tol",1.0e-3)
Ipopt.addOption(prob,"tol",1.0e-3)

prob.x .= 1.0*rand(nx)
@time solvestat = Ipopt.solveProblem(prob)

b_sol = prob.x[1:2]
ψ_sol = prob.x[3]
s_sol = prob.x[4]

norm(v[1:2])
vec_norm(v[1:2])'*vec_norm(b_sol)

con!(_c,prob.x)
_c
_C = zeros(nc,nx)
∇con!(_C,prob.x)
_C


text:
struct CompletePolyDR{S<:Grid,T<:Grid,nx} <: AbstractDecisionRule{S,T,nx}
    grid_exo::S
    grid_endo::T
    coefs::Vector{Matrix{Float64}}
    order::Int
end

#####
##### 1-argument decision rule
#####

function CompletePolyDR(
        grid_exo::EmptyGrid, grid_endo::Grid{ns},
        ::Union{Val{nx},Type{Val{nx}}}, order::Int=3
    ) where ns where nx
    coeffs = [Array{Float64}(undef, BM.n_complete(ns, order), nx)]
    CompletePolyDR{EmptyGrid,typeof(grid_endo),nx}(grid_exo, grid_endo, coeffs, order)
end

function set_values!(
        dr::CompletePolyDR{<:G}, values::Vector{Matrix{Float64}}
    ) where G <: Union{<:EmptyGrid,<:UnstructuredGrid}
    # TODO the following should be once for all with the result stored in dr
    B_grid = BM.complete_polynomial(nodes(Matrix,dr.grid_endo), dr.order)
    for i in 1:length(values)
        dr.coefs[i][:,:] = B_grid\values[i]
        # TODO: understand why the following seems to modify values[i]
        # q_B_grid = qr(B_grid, Val(true))
        # ldiv!(dr.coefs[i], q_B_grid, values[i])
    end
end

function set_values!(
        dr::CompletePolyDR{<:G,<:Grid,nx},
        values::Vector{<:Array{Value{nx}}}
    ) where G <: Union{EmptyGrid,UnstructuredGrid} where nx
    B_grid = BM.complete_polynomial(nodes(Matrix,dr.grid_endo), dr.order)
    # q_B_grid = qr(B_grid, Val(true))

    if length(values) != length(dr.coefs)
        msg = "The length of values ($(length(values))) is not the same "
        msg *= "as the length of the coefficient Vector ($(length(dr.coefs)))"
        error(msg)
    end

    for i in 1:length(values)
        N = length(values[i])
        # data = copy(reshape(reinterpret(Float64, vec(values[i])), (nx, N))')
        data = copy(from_LOP(values[i]))
        dr.coefs[i] = B_grid \ data
        # ldiv!(dr.coefs[i], q_B_grid, data)
    end
end

function evaluate(dr::CompletePolyDR{<:EmptyGrid}, z::AbstractMatrix{Float64})
    B = BM.complete_polynomial(z, dr.order)
    B*dr.coefs[1]
end

function evaluate(dr::CompletePolyDR{<:EmptyGrid,<:Grid{d}}, points::AbstractVector{Point{d}}) where d
    N = length(points)
    mat = copy(reshape(reinterpret(Float64, vec(points)), (d, N))')
    evaluate(dr, mat)
end

####
#### UnstructuredGrid × CartesianGrid 2 continous arguments d.r.
####

function CompletePolyDR(
        grid_exo::S, grid_endo::Grid{ns},
        ::Union{Val{nx},Type{Val{nx}}}, order::Int=3
    ) where S <: UnstructuredGrid where ns where nx
    n_coefs = BM.n_complete(ns, order)
    coeffs = [Array{Float64}(undef, n_coefs, nx) for i in 1:n_nodes(grid_exo)]
    CompletePolyDR{S,typeof(grid_endo),nx}(grid_exo, grid_endo, coeffs, order)
end

function evaluate(dr::CompletePolyDR{<:UnstructuredGrid}, i::Int, z::AbstractMatrix{Float64})
    @boundscheck begin
        n_funcs = length(dr.coefs)
        if i > n_funcs
            msg = "Only $n_funcs are known, but function $i was requested"
            throw(BoundsError(msg))
        end
    end
    B = BM.complete_polynomial(z, dr.order)
    B*dr.coefs[i]
end

function evaluate(dr::CompletePolyDR{<:UnstructuredGrid}, i::Int, z::AbstractVector{Point{d}}) where d
    N = length(z)
    mat = reshape(reinterpret(Float64, vec(z)), (d, N))'
    evaluate(dr, i, mat)
end


text:
# Douglas-Rachford line-search method

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct DRLS_iterable{
    R, C <: Union{R, Complex{R}}, Tx <: AbstractArray{C}, Tf, Tg, TH
}
    f::Tf
    g::Tg
    x0::Tx
    lambda::R
    gamma::R
    c::R
    N::Int
    H::TH
end

Base.IteratorSize(::Type{<:DRLS_iterable}) = Base.IsInfinite()

mutable struct DRLS_state{R, Tx, TH}
    x::Tx
    u::Tx
    v::Tx
    w::Tx
    res::Tx
    res_prev::Tx
    xbar::Tx
    xbar_prev::Tx
    d::Tx
    x_d::Tx
    gamma::R
    f_u::R
    g_v::R
    H::TH
    tau::Maybe{R}
end

function DRE(state::DRLS_state)
    return (
        state.f_u + state.g_v
        - real(dot(state.x - state.u, state.res))/state.gamma
        + 1/(2*state.gamma) * norm(state.res)^2
    )
end

function Base.iterate(iter::DRLS_iterable)
    x = iter.x0
    u, f_u = prox(iter.f, x, iter.gamma)
    w = 2 * u - x
    v, g_v = prox(iter.g, w, iter.gamma)
    res = u - v
    xbar = x - iter.lambda * res
    state = DRLS_state(
        x, u, v, w, res, zero(x), xbar, zero(x), zero(x), zero(x), iter.gamma,
        f_u, g_v, iter.H, nothing
    )
    return state, state
end

function Base.iterate(iter::DRLS_iterable{R}, state::DRLS_state) where R
    DRE_curr = DRE(state)

    mul!(state.d, iter.H, -state.res)
    state.x_d .= state.x .+ state.d
    copyto!(state.xbar_prev, state.xbar)
    copyto!(state.res_prev, state.res)
    state.tau = R(1)
    state.x .= state.x_d

    for k in 1:iter.N
        state.f_u = prox!(state.u, iter.f, state.x, iter.gamma)
        state.w .= 2 .* state.u .- state.x
        state.g_v = prox!(state.v, iter.g, state.w, iter.gamma)
        state.res .= state.u .- state.v

        if k == 1
            update!(iter.H, state.d, state.res - state.res_prev)
        end

        state.xbar .= state.x .- iter.lambda * state.res
        DRE_candidate = DRE(state)

        if DRE_candidate <= DRE_curr - iter.c / iter.gamma * norm(state.res)^2
            return state, state
        end

        state.tau = state.tau / 2
        state.x .= state.tau .* state.x_d .+ (1 - state.tau) .* state.xbar_prev
    end

    @warn "stepsize `tau` became too small ($(state.tau)), stopping the iterations"
    return nothing
end

# Solver

struct DRLS{R}
    alpha::R
    beta::R
    gamma::Maybe{R}
    lambda::R
    N::Int
    memory::Int
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function DRLS{R}(; alpha::R=R(0.95), beta::R=R(0.5),
        gamma::Maybe{R}=nothing, lambda::R=R(1), N::Int=20, memory::Int=5,
        maxit::Int=1000, tol::R=R(1e-8), verbose::Bool=false, freq::Int=10
    ) where R
        @assert 0 < alpha < 1
        @assert 0 < beta < 1
        @assert gamma === nothing || gamma > 0
        @assert 0 < lambda < 2
        @assert N > 0
        @assert memory >= 0
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(alpha, beta, gamma, lambda, N, memory, maxit, tol, verbose, freq)
    end
end

function (solver::DRLS{R})(
    x0::AbstractArray{C}; f=Zero(), g=Zero(), L::Maybe{R}=nothing
) where {R, C <: Union{R, Complex{R}}}
    stop(state::DRLS_state) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it, state.gamma / state.gamma, norm(state.res, Inf),
        (state.tau === nothing ? 0.0 : state.tau)
    )

    gamma = if solver.gamma === nothing && L !== nothing
        if ProximalOperators.is_convex(f)
            solver.alpha / L
        else
            solver.alpha * (2 - solver.lambda) / (2 * L)
        end
    else
        solver.gamma
    end

    m = if ProximalOperators.is_convex(f)
        max(gamma * L - solver.lambda / 2, 0)
    else
        1
    end
    C_gamma_lambda = (
        solver.lambda / ((1 + gamma * L) ^ 2) *
        ((2 - solver.lambda) / 2 - gamma * L * m)
    )
    c = solver.beta * C_gamma_lambda

    iter = DRLS_iterable(
        f, g, x0, solver.lambda, gamma, c, solver.N, LBFGS(x0, solver.memory)
    )
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.u, state_final.v, num_iters
end

DRLS(::Type{R}; kwargs...) where R = DRLS{R}(; kwargs...)
DRLS(; kwargs...) = DRLS(Float64; kwargs...)


text:
convex_hull(a::Matrix{<: Real}) = convex_hull(collect.(eachcol(a)))
function convex_hull(a::Array{Vector{T},1} where T<: Real)
    cw(a, b, c) = (a[1]*(b[2]-c[2])+b[1]*(c[2]-a[2])+c[1]*(a[2]-b[2]) < 0);
    ccw(a, b, c) = (a[1]*(b[2]-c[2])+b[1]*(c[2]-a[2])+c[1]*(a[2]-b[2]) > 0);

    if length(a) == 1
        return a
    end

    a = sort(a, lt = (a,b) -> (a[1] < b[1] || a[1] == b[1] && a[2] < b[2]))
    p1 = a[1];
    p2 = a[end];
    up = [p1]
    down = [p1]

    for i in 2:length(a)
        if (i==length(a) || cw(p1, a[i], p2))
            while (length(up) >= 2 && !cw(up[end-1], up[end], a[i]))
                up = up[1:end-1]
            end

            push!(up, a[i])
        end
        if (i==length(a) || ccw(p1, a[i], p2))
            while (length(down) >= 2 && !ccw(down[end-1], down[end], a[i]))
                down = down[1:end-1]
            end

            push!(down, a[i])
        end
    end

    return hcat(vcat(up, reverse(down[1:end-1]))...);
end

function area(polygon::Matrix{<:Real})
    area = 0.0;

    for i in 1:size(polygon, 2)
        j = i % size(polygon, 2) + 1;
        area += polygon[1, i] * polygon[2, j] - polygon[2, i] * polygon[1, j];
    end

   return abs(area) / 2;
end


text:
using DataFrames
using CSV
include("../src/TrueSkill.jl")
using .TrueSkill
global const ttt = TrueSkill
using Test

p_d_m = [log(0.5)]
p_d_m_trueskill = [log(0.5)]
p_d_m_hat = [log(0.5)]
loocv_hat = [log(0.5)]
composition = Vector{Vector{Vector{String}}}()
push!(composition, [["b"],["a"]])
h = ttt.History(composition)
for i in 1:50
    lc = ttt.learning_curves(h)
    push!(p_d_m, p_d_m[end]+log(ttt.Game([[ttt.Player(lc["a"][end][2])],[ttt.Player(lc["b"][end][2])]]).evidence))
    push!(composition, [["a"],["b"]])
    h = ttt.History(composition)
    push!(p_d_m_trueskill,ttt.log_evidence(h))
    ttt.convergence(h,iterations=10)
    push!(loocv_hat,ttt.log_evidence(h))
    push!(p_d_m_hat,ttt.log_evidence(h, online=true))

    lc = ttt.learning_curves(h)
    push!(p_d_m, p_d_m[end]+log(ttt.Game([[ttt.Player(lc["b"][end][2])],[ttt.Player(lc["b"][end][2])]]).evidence))
    push!(composition, [["b"],["a"]])
    h = ttt.History(composition)
    push!(p_d_m_trueskill,ttt.log_evidence(h))
    ttt.convergence(h,iterations=10)
    push!(loocv_hat,ttt.log_evidence(h))
    push!(p_d_m_hat,ttt.log_evidence(h, online=true))
end

@test sum(p_d_m_trueskill .-1e-4 .<= p_d_m) == length(p_d_m)

@test sum(p_d_m_hat .-1e-4 .<= p_d_m) == length(p_d_m)

#@test sum(loocv_hat .-1e-4  .<= p_d_m ) == length(p_d_m)

df = DataFrame()
df[!,:p_d_m] = p_d_m
df[!,:p_d_m_trueskill] = p_d_m_trueskill
df[!,:p_d_m_hat] = p_d_m_hat
df[!,:loocv_hat] = loocv_hat

CSV.write("output/prior_predictions.csv", df; header=true)

# using Plots
# plot(loocv_hat-p_d_m)
# plot!(p_d_m_hat-p_d_m)
# plot!(p_d_m_trueskill-p_d_m)


text:
using Conda
using Base.Test
using Compat

already_installed = "curl" in Conda._installed_packages()

@test Conda.exists("curl")
Conda.add("curl")

if is_unix()
    curl_path = joinpath(Conda.PREFIX, "bin", "curl-config")
end
if is_windows()
    using BinDeps
    manager = Conda.Manager(["curl"])
    curl_libpath = BinDeps.libdir(manager, "")
    curl_path = joinpath(curl_libpath, "curl.exe")
end

@test isfile(curl_path)

@test isfile(joinpath(Conda.BINDIR, basename(curl_path)))

Conda.rm("curl")
if is_unix()
    @test !isfile(curl_path)
end

if already_installed
    Conda.add("curl")
end

@test isfile(joinpath(Conda.SCRIPTDIR, "conda" * (is_windows() ? ".exe": "")))

@test isfile(joinpath(Conda.PYTHONDIR, "python" * (is_windows() ? ".exe": "")))

channels = Conda.channels()
@test (isempty(channels) || channels == ["defaults"])

Conda.add_channel("foo")
@test Conda.channels() == ["foo", "defaults"]

Conda.rm_channel("foo")
channels = Conda.channels()
@test (isempty(channels) || channels == ["defaults"])

if is_windows()
    Conda.add("qt")
    Conda.rm("qt")
end


text:
using GLMakie

Base.@kwdef mutable struct Lorenz
    dt::Float64 = 0.01
    σ::Float64 = 10
    ρ::Float64 = 28
    β::Float64 = 8/3
    x::Float64 = 1
    y::Float64 = 1
    z::Float64 = 1
end

function step!(l::Lorenz)
    dx = l.σ * (l.y - l.x)
    dy = l.x * (l.ρ - l.z) - l.y
    dz = l.x * l.y - l.β * l.z
    l.x += l.dt * dx
    l.y += l.dt * dy
    l.z += l.dt * dz
    Point3f(l.x, l.y, l.z)
end

attractor = Lorenz()

points = Observable(Point3f[])
colors = Observable(Int[])

set_theme!(theme_black())

fig, ax, l = lines(points, color = colors,
    colormap = :inferno, transparency = true,
    axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
        viewmode = :fit, limits = (-30, 30, -30, 30, 0, 50)))

frames = 1:120
for frame in frames
    for i in 1:50
        push!(points[], step!(attractor))
        push!(colors[], frame)
    end
    ax.azimuth[] = 1.7pi + 0.3 * sin(2pi * frame / 120)
    notify.((points, colors))
    l.colorrange = (0, frame)
    sleep(0.1)
    #println(frame)
end

max_stars_repo_name: m-geo/lorenz

max_issues_repo_name: m-geo/lorenz

max_forks_repo_name: m-geo/lorenz


text:
using Revise
using DifferentialEquations
using DynamicalSystems
using Basins
using Plots
using Printf

ds = Systems.henon(zeros(2); a = 1.4, b = 0.3)
integ_df  = integrator(ds)

xres=200
yres=200

# range for forced pend
xg = range(-2.,2.,length=xres)
yg = range(-2.,2.,length=yres)

# compute basin
@time basin = Basins.basins_map2D(xg, yg, integ_df)

# Basin entropy
@show Sb,Sbb = basin_entropy(bsn.basin; eps_x=20, eps_y=20)

# Wada merge Haussdorff distances
@time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn.basin)
epsilon = xg[2]-xg[1]
@show dmax = max_dist/epsilon
@show dmin = min_dist/epsilon

# Wada grid
W = detect_wada_grid_method(integ_df, bsn; max_iter=8)
@show W[:,end]

# Uncertainty exponent for these parameter and grid
bd = box_counting_dim(xg, yg, bsn.basin)
α = 2 - bd

println("---------------")
println("---------------")
println("Basin Report: ")
println("---------------")
println("---------------")

@printf("Basin entropy %.2f \n", Sb)
@printf("Boundary Basin Entropy: %.2f\n", Sbb)
@printf("Uncertainty exponent: α= %.2f\n", α )
@printf("Number of basins: %d\n", length(unique(bsn.basin)))
@printf("Merge Method: Max fattening parameter: %.2f\n", dmax)
@printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] )

plot(xg,yg,bsn.basin', seriestype=:heatmap)


text:
abstract type AbstractGraphLayer end

(l::AbstractGraphLayer)(x::AbstractMatrix) = l(l.fg, x)


text:
function delta(n)
# DELTA calculates the sun's declination angle
# y = DELTA(jd) gets the Julian day (jd) and calculate declination of the sun

d=(asin(-sin(23.45/180*pi).*cos(360/365.25*(n+10)/180*pi)));

return d

end


text:
using Test
using Qecsim.BasicModels
using Qecsim.Model
using Qecsim.PauliTools:to_bsf

@testset "BasicCode" begin
    pauli_stabilizers = ["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"]
    pauli_logical_xs = ["XXXXX"]
    pauli_logical_zs = ["ZZZZZ"]
    my_nkd = (5, 1, 3)
    my_label = "5-qubit"
    code = BasicCode(pauli_stabilizers, pauli_logical_xs, pauli_logical_zs, my_nkd,
                     my_label)
    @test nkd(code) == my_nkd
    @test label(code) == my_label
    @test stabilizers(code) == to_bsf(pauli_stabilizers)
    @test logical_xs(code) == to_bsf(pauli_logical_xs)
    @test logical_zs(code) == to_bsf(pauli_logical_zs)
    # defaults
    code = BasicCode(pauli_stabilizers, pauli_logical_xs, pauli_logical_zs)
    @test isequal(nkd(code), (5, 1, missing))  # isequal equates missing
    @test isa(label(code), String)
    # 5-qubit
    code = FiveQubitCode()
    @test nkd(code) == (5, 1, 3)
    @test label(code) == "5-qubit"
    @test validate(code) === nothing  # no error
    # Steane
    code = SteaneCode()
    @test nkd(code) == (7, 1, 3)
    @test label(code) == "Steane"
    @test validate(code) === nothing  # no error
end


text:
module HCEstimator

include("Tools.jl")
using .Tools
include("DistSystem.jl")
using .DistSystem
include("Estimator.jl")
using .Estimator

export DistSystem, build_model, Tools
import JuMP

function build_model(model::JuMP.Model, sys, obj=true)
    if isempty(sys.dgs)
        push!(sys.dgs, DistSystem.null_der(sys.buses))
    end
    model = nl_pf(model, sys)
    if obj
        model = add_objective(model, sys)
    end
    return model
end

end


text:
export maketankeditor, showtankeditor, tank_liststore, treeview

using SQLite
using DataFrames
using Gtk
using Rockets

function showtankeditor()
    db = opendb()
    list = tank_liststore(db)
    win = maketankeditor(list, db)
    showall(win)
    return win, db
end

function maketankeditor(tanklist::GtkListStore, db::DBInterface.Connection)
    win = GtkWindow("Tank Editor", 500, 250)

    toplayout = GtkBox(:h)  # :h makes a horizontal layout, :v a vertical layout
    set_gtk_property!(toplayout, :margin_left, 10)
    set_gtk_property!(toplayout, :margin_top, 10)
    set_gtk_property!(toplayout,:spacing,10)
    push!(win, toplayout)

    # Master View Containing list of tanks
    master_layout = GtkBox(:v)
    set_gtk_property!(master_layout,:spacing,10)
    push!(toplayout, master_layout)

    model = GtkTreeModel(tanklist)
    tree = treeview(model)

    push!(master_layout, tree)

    actions_layout = GtkBox(:h)
    set_gtk_property!(actions_layout,:spacing,10)
    push!(master_layout, actions_layout)

    addbtn = GtkButton("Add")
    removebtn = GtkButton("Remove")
    reloadbtn = GtkButton("Reload")

    push!(actions_layout, addbtn)
    push!(actions_layout, removebtn)
    push!(actions_layout, reloadbtn)

    # Detail View with info about each selected tank
    detail_layout = GtkGrid()

    namelabel = GtkLabel("Name:")
    totallabel = GtkLabel("Total mass:")
    drylabel = GtkLabel("Dry mass:")

    nameentry = GtkEntry()
    totalentry = GtkEntry()
    dryentry = GtkEntry()

    detail_layout[1, 1] = namelabel
    detail_layout[2, 1] = nameentry
    detail_layout[1, 2] = drylabel
    detail_layout[2, 2] = dryentry
    detail_layout[1, 3] = totallabel
    detail_layout[2, 3] = totalentry

    set_gtk_property!(detail_layout, :column_spacing, 10)
    set_gtk_property!(detail_layout, :row_spacing, 10)

    push!(toplayout, detail_layout)

    # Deal with selection and model update
    selection = GAccessor.selection(tree)
    chid = signal_connect(selection, "changed") do widget
        selection = GAccessor.selection(tree)
        if hasselection(selection)
            current = selected(selection)
            tankname = tanklist[current, 1]
            tank = findfirsttank(db, tankname)

            if !isnothing(tank)
                set_gtk_property!(nameentry, :text, tankname)
                set_gtk_property!(totalentry, :text, tank.total_mass)
                set_gtk_property!(dryentry, :text, tank.dry_mass)
            end
        end
    end

    # To disconnect this signal handler later
    # signal_handler_disconnect(selection, chid)

    # Adding row
    addid = signal_connect(addbtn, "clicked") do widget
        name = get_gtk_property(nameentry, :text, String)

        # Cannot get text as numbers. Need to make that conversion ourselves
        drytxt = get_gtk_property(dryentry, :text, String)
        totaltxt = get_gtk_property(totalentry, :text, String)

        dry = parse(Float64, drytxt)
        total = parse(Float64, totaltxt)

        add_tank!(db, name, Tank(dry, total))
        push!(tanklist, (name,))
    end

    # Remove row
    rmid = signal_connect(removebtn, "clicked") do widget
        selection = GAccessor.selection(tree)
        if hasselection(selection)
            current = selected(selection)
            tankname = tanklist[current, 1]
            tank = findfirsttank(db, tankname)

            if !isnothing(tank)
                # Put on this on main GUI task. I think.
                # GUI stuff cannot be run on other tasks
                Gtk.@sigatom begin
                    deleteat!(tanklist, current)

                    set_gtk_property!(nameentry, :text, "")
                    set_gtk_property!(totalentry, :text, "")
                    set_gtk_property!(dryentry, :text, "")
                end
            end
        end
    end

    # Reloading
    reloadid = signal_connect(reloadbtn, "clicked") do widget
        Gtk.@sigatom begin
            empty!(tanklist)
            for name in tanknames(db)
               push!(tanklist, (name,))
            end
        end
    end

    return win
end

function tank_liststore(db::DBInterface.Connection)
    list = GtkListStore(String)

    for name in tanknames(db)
       push!(list, (name,))
    end

    list
end

function treeview(model::GtkTreeModel)
    tree = GtkTreeView(model)

    txtcell = GtkCellRendererText()

    # Render first column as text
    col1 = GtkTreeViewColumn("Tank", txtcell, Dict([("text",0)]))

    push!(tree, col1)

    return tree
end

"Just meant for experimenting with list view"
function liststore_window_test()
    ls = GtkListStore(String, Int, Bool)
    push!(ls,("Peter",20,false))
    push!(ls,("Paul",30,false))
    push!(ls,("Mary",25,true))

    m  = GtkTreeModel(ls)
    tv = GtkTreeView(m)

    rTxt = GtkCellRendererText()
    rTog = GtkCellRendererToggle()

    c1 = GtkTreeViewColumn("Name", rTxt, Dict([("text",0)]))
    c2 = GtkTreeViewColumn("Age", rTxt, Dict([("text",1)]))
    c3 = GtkTreeViewColumn("Female", rTog, Dict([("active",2)]))

    push!(tv, c1, c2, c3)

    # Get selection
    selection = GAccessor.selection(tv)

    win = GtkWindow(tv, "List View")
    showall(win)
end


text:
@testset "Test Microscope mode" begin
    @test Microscopy.MicroscopeConfocal <: Microscopy.Microscope
    @test Microscopy.MicroscopeWidefield <: Microscopy.Microscope
end

max_stars_repo_name: JuliaMicroscopy/Microscopy

max_forks_repo_name: JuliaMicroscopy/Microscopy


text:
module AtomModule
export Atom, Element
import Base.==

using ..BaseModule

atomicNumbers = Dict{AbstractString,Integer}(
  "H" => 	1,
  "He" => 	2,
  "Li" => 	3,
  "Be" =>	4,
  "B" =>	5,
  "C" =>	6,
  "N" =>	7,
  "O" =>	8,
  "F" =>	9,
  "Ne" =>	10,
  "Na" =>	11,
  "Mg" =>	12,
  "Al" =>	13,
  "Si" =>	14,
  "P" =>	15,
  "S" =>	16,
  "Cl" =>	17,
  "Ar" =>	18,
  "K" =>	19,
  "Ca" =>	20,
  "Sc" =>	21,
  "Ti" =>	22,
  "V" =>	23,
  "Cr" =>	24,
  "Mn" =>	25,
  "Fe" =>	26,
  "Co" =>	27,
  "Ni" =>	28,
  "Cu" =>	29,
  "Zn" =>	30,
  "Ga" =>	31,
  "Ge" =>	32,
  "As" =>	33,
  "Se" =>	34,
  "Br" =>	35,
  "Kr" =>	36,
  "Rb" =>	37,
  "Sr" =>	38,
  "Y" =>	39,
  "Zr" =>	40,
  "Nb" =>	41,
  "Mo" =>	42,
  "Tc" =>	43,
  "Ru" =>	44,
  "Rh" =>	45,
  "Pd" =>	46,
  "Ag" =>	47,
  "Cd" =>	48,
  "In" =>	49,
  "Sn" =>	50,
  "Sb" =>	51,
  "Te" =>	52,
  "I" =>	53,
  "Xe" =>	54,
  "Cs" =>	55,
  "Ba" =>	56,
  "La" =>	57,
  "Ce" =>	58,
  "Pr" =>	59,
  "Nd" =>	60,
  "Pm" =>	61,
  "Sm" =>	62,
  "Eu" =>	63,
  "Gd" =>	64,
  "Tb" =>	65,
  "Dy" =>	66,
  "Ho" =>	67,
  "Er" =>	68,
  "Tm" =>	69,
  "Yb" =>	70,
  "Lu" =>	71,
  "Hf" =>	72,
  "Ta" =>	73,
  "W" =>	74,
  "Re" =>	75,
  "Os" =>	76,
  "Ir" =>	77,
  "Pt" =>	78,
  "Au" =>	79,
  "Hg" =>	80,
  "Tl" =>	81,
  "Pb" =>	82,
  "Bi" =>	83,
  "Po" =>	84,
  "At" =>	85,
  "Rn" =>	86,
  "Fr" =>	87,
  "Ra" =>	88,
  "Ac" =>	89,
  "Th" =>	90,
  "Pa" =>	91,
  "U" =>	92,
  "Np" =>	93,
  "Pu" =>	94,
  "Am" =>	95,
  "Cm" =>	96,
  "Bk" =>	97,
  "Cf" =>	98,
  "Es" =>	99,
  "Fm" =>	100,
  "Md" =>	101,
  "No" =>	102,
  "Lr" =>	103,
  "Rf" =>	104,
  "Db" =>	105,
  "Sg" =>	106,
  "Bh" =>	107,
  "Hs" =>	108,
  "Mt" =>	109,
  "Ds" =>	110,
  "Rg" =>	111,
  "Cp" =>	112,
  "Uut" =>	113,
  "Uuq" =>	114,
  "Uup" =>	115,
  "Uuh" =>	116,
  "Uus" =>	117,
  "Uuo" =>	118
)

immutable Element
	symbol::AbstractString
	atomicNumber::Int

	Element(sym::AbstractString) = new(sym,atomicNumbers[sym])
end

==(el1::Element,el2::Element) = el1.symbol==el2.symbol && el1.atomicNumber==el2.atomicNumber
isequal(el1::Element,el2::Element) = isequal(el1.symbol,el2.symbol) && isequal(el1.atomicNumber,el2.atomicNumber)

import Base.hash
hash(el::Element) = hash(el.symbol,hash(el.atomicNumber))

type Atom
	element::Element
	position::Position
end

end # module

max_stars_repo_name: tkelman/QuantumLab

max_issues_repo_name: tkelman/QuantumLab

max_forks_repo_name: tkelman/QuantumLab


text:
# Automatically generated using Clang.jl wrap_c

const NFRAMES = 3
const SECRET_NUM = -1234

@cenum(LAYER_TYPE,
    CONVOLUTIONAL = 0,
    DECONVOLUTIONAL = 1,
    CONNECTED = 2,
    MAXPOOL = 3,
    SOFTMAX = 4,
    DETECTION = 5,
    DROPOUT = 6,
    CROP = 7,
    ROUTE = 8,
    COST = 9,
    NORMALIZATION = 10,
    AVGPOOL = 11,
    LOCAL = 12,
    SHORTCUT = 13,
    ACTIVE = 14,
    RNN = 15,
    GRU = 16,
    LSTM = 17,
    CRNN = 18,
    BATCHNORM = 19,
    NETWORK = 20,
    XNOR = 21,
    REGION = 22,
    YOLO = 23,
    ISEG = 24,
    REORG = 25,
    REORG_OLD = 26,
    UPSAMPLE = 27,
    LOGXENT = 28,
    L2NORM = 29,
    BLANK = 30,
)
@cenum(ACTIVATION,
    LOGISTIC = 0,
    RELU = 1,
    RELIE = 2,
    LINEAR = 3,
    RAMP = 4,
    TANH = 5,
    PLSE = 6,
    LEAKY = 7,
    ELU = 8,
    LOGGY = 9,
    STAIR = 10,
    HARDTAN = 11,
    LHTAN = 12,
    SELU = 13,
)
@cenum(COST_TYPE,
    SSE = 0,
    MASKED = 1,
    L1 = 2,
    SEG = 3,
    SMOOTH = 4,
    WGAN = 5,
)

struct tree
    leaf::Ptr{Cint}
    n::Cint
    parent::Ptr{Cint}
    child::Ptr{Cint}
    group::Ptr{Cint}
    name::Ptr{Cstring}
    groups::Cint
    group_size::Ptr{Cint}
    group_offset::Ptr{Cint}
end

struct layer
    type::LAYER_TYPE
    activation::ACTIVATION
    cost_type::COST_TYPE
    forward::Ptr{Cvoid}
    backward::Ptr{Cvoid}
    update::Ptr{Cvoid}
    forward_gpu::Ptr{Cvoid}
    backward_gpu::Ptr{Cvoid}
    update_gpu::Ptr{Cvoid}
    batch_normalize::Cint
    shortcut::Cint
    batch::Cint
    forced::Cint
    flipped::Cint
    inputs::Cint
    outputs::Cint
    nweights::Cint
    nbiases::Cint
    extra::Cint
    truths::Cint
    h::Cint
    w::Cint
    c::Cint
    out_h::Cint
    out_w::Cint
    out_c::Cint
    n::Cint
    max_boxes::Cint
    groups::Cint
    size::Cint
    side::Cint
    stride::Cint
    reverse::Cint
    flatten::Cint
    spatial::Cint
    pad::Cint
    sqrt::Cint
    flip::Cint
    index::Cint
    binary::Cint
    xnor::Cint
    use_bin_output::Cint
    steps::Cint
    hidden::Cint
    truth::Cint
    smooth::Cfloat
    dot::Cfloat
    angle::Cfloat
    jitter::Cfloat
    saturation::Cfloat
    exposure::Cfloat
    shift::Cfloat
    ratio::Cfloat
    learning_rate_scale::Cfloat
    clip::Cfloat
    focal_loss::Cint
    noloss::Cint
    softmax::Cint
    classes::Cint
    coords::Cint
    background::Cint
    rescore::Cint
    objectness::Cint
    does_cost::Cint
    joint::Cint
    noadjust::Cint
    reorg::Cint
    log::Cint
    tanh::Cint
    mask::Ptr{Cint}
    total::Cint
    bflops::Cfloat
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    t::Cint
    alpha::Cfloat
    beta::Cfloat
    kappa::Cfloat
    coord_scale::Cfloat
    object_scale::Cfloat
    noobject_scale::Cfloat
    mask_scale::Cfloat
    class_scale::Cfloat
    bias_match::Cint
    random::Cint
    ignore_thresh::Cfloat
    truth_thresh::Cfloat
    thresh::Cfloat
    focus::Cfloat
    classfix::Cint
    absolute::Cint
    onlyforward::Cint
    stopbackward::Cint
    dontload::Cint
    dontsave::Cint
    dontloadscales::Cint
    numload::Cint
    temperature::Cfloat
    probability::Cfloat
    scale::Cfloat
    cweights::Cstring
    indexes::Ptr{Cint}
    input_layers::Ptr{Cint}
    input_sizes::Ptr{Cint}
    map::Ptr{Cint}
    counts::Ptr{Cint}
    sums::Ptr{Ptr{Cfloat}}
    rand::Ptr{Cfloat}
    cost::Ptr{Cfloat}
    state::Ptr{Cfloat}
    prev_state::Ptr{Cfloat}
    forgot_state::Ptr{Cfloat}
    forgot_delta::Ptr{Cfloat}
    state_delta::Ptr{Cfloat}
    combine_cpu::Ptr{Cfloat}
    combine_delta_cpu::Ptr{Cfloat}
    concat::Ptr{Cfloat}
    concat_delta::Ptr{Cfloat}
    binary_weights::Ptr{Cfloat}
    biases::Ptr{Cfloat}
    bias_updates::Ptr{Cfloat}
    scales::Ptr{Cfloat}
    scale_updates::Ptr{Cfloat}
    weights::Ptr{Cfloat}
    weight_updates::Ptr{Cfloat}
    align_bit_weights_gpu::Cstring
    mean_arr_gpu::Ptr{Cfloat}
    align_workspace_gpu::Ptr{Cfloat}
    transposed_align_workspace_gpu::Ptr{Cfloat}
    align_workspace_size::Cint
    align_bit_weights::Cstring
    mean_arr::Ptr{Cfloat}
    align_bit_weights_size::Cint
    lda_align::Cint
    new_lda::Cint
    bit_align::Cint
    col_image::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    output::Ptr{Cfloat}
    delta_pinned::Cint
    output_pinned::Cint
    loss::Ptr{Cfloat}
    squared::Ptr{Cfloat}
    norms::Ptr{Cfloat}
    spatial_mean::Ptr{Cfloat}
    mean::Ptr{Cfloat}
    variance::Ptr{Cfloat}
    mean_delta::Ptr{Cfloat}
    variance_delta::Ptr{Cfloat}
    rolling_mean::Ptr{Cfloat}
    rolling_variance::Ptr{Cfloat}
    x::Ptr{Cfloat}
    x_norm::Ptr{Cfloat}
    m::Ptr{Cfloat}
    v::Ptr{Cfloat}
    bias_m::Ptr{Cfloat}
    bias_v::Ptr{Cfloat}
    scale_m::Ptr{Cfloat}
    scale_v::Ptr{Cfloat}
    z_cpu::Ptr{Cfloat}
    r_cpu::Ptr{Cfloat}
    h_cpu::Ptr{Cfloat}
    prev_state_cpu::Ptr{Cfloat}
    temp_cpu::Ptr{Cfloat}
    temp2_cpu::Ptr{Cfloat}
    temp3_cpu::Ptr{Cfloat}
    dh_cpu::Ptr{Cfloat}
    hh_cpu::Ptr{Cfloat}
    prev_cell_cpu::Ptr{Cfloat}
    cell_cpu::Ptr{Cfloat}
    f_cpu::Ptr{Cfloat}
    i_cpu::Ptr{Cfloat}
    g_cpu::Ptr{Cfloat}
    o_cpu::Ptr{Cfloat}
    c_cpu::Ptr{Cfloat}
    dc_cpu::Ptr{Cfloat}
    binary_input::Ptr{Cfloat}
    bin_re_packed_input::Ptr{UInt32}
    t_bit_input::Cstring
    input_layer::Ptr{layer}
    self_layer::Ptr{layer}
    output_layer::Ptr{layer}
    reset_layer::Ptr{layer}
    update_layer::Ptr{layer}
    state_layer::Ptr{layer}
    input_gate_layer::Ptr{layer}
    state_gate_layer::Ptr{layer}
    input_save_layer::Ptr{layer}
    state_save_layer::Ptr{layer}
    input_state_layer::Ptr{layer}
    state_state_layer::Ptr{layer}
    input_z_layer::Ptr{layer}
    state_z_layer::Ptr{layer}
    input_r_layer::Ptr{layer}
    state_r_layer::Ptr{layer}
    input_h_layer::Ptr{layer}
    state_h_layer::Ptr{layer}
    wz::Ptr{layer}
    uz::Ptr{layer}
    wr::Ptr{layer}
    ur::Ptr{layer}
    wh::Ptr{layer}
    uh::Ptr{layer}
    uo::Ptr{layer}
    wo::Ptr{layer}
    uf::Ptr{layer}
    wf::Ptr{layer}
    ui::Ptr{layer}
    wi::Ptr{layer}
    ug::Ptr{layer}
    wg::Ptr{layer}
    softmax_tree::Ptr{tree}
    workspace_size::Cint
end

@cenum(learning_rate_policy,
    CONSTANT = 0,
    STEP = 1,
    EXP = 2,
    POLY = 3,
    STEPS = 4,
    SIG = 5,
    RANDOM = 6,
)

struct network
    n::Cint
    batch::Cint
    seen::Ptr{UInt64}
    t::Ptr{Cint}
    epoch::Cfloat
    subdivisions::Cint
    layers::Ptr{layer}
    output::Ptr{Cfloat}
    policy::learning_rate_policy
    learning_rate::Cfloat
    momentum::Cfloat
    decay::Cfloat
    gamma::Cfloat
    scale::Cfloat
    power::Cfloat
    time_steps::Cint
    step::Cint
    max_batches::Cint
    scales::Ptr{Cfloat}
    steps::Ptr{Cint}
    num_steps::Cint
    burn_in::Cint
    cudnn_half::Cint
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    inputs::Cint
    outputs::Cint
    truths::Cint
    notruth::Cint
    h::Cint
    w::Cint
    c::Cint
    max_crop::Cint
    min_crop::Cint
    max_ratio::Cfloat
    min_ratio::Cfloat
    center::Cint
    flip::Cint
    angle::Cfloat
    aspect::Cfloat
    exposure::Cfloat
    saturation::Cfloat
    hue::Cfloat
    random::Cint
    track::Cint
    augment_speed::Cint
    try_fix_nan::Cint
    gpu_index::Cint
    hierarchy::Ptr{tree}
    input::Ptr{Cfloat}
    truth::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    workspace::Ptr{Cfloat}
    train::Cint
    index::Cint
    cost::Ptr{Cfloat}
    clip::Cfloat
end

struct network_state
    truth::Ptr{Cfloat}
    input::Ptr{Cfloat}
    delta::Ptr{Cfloat}
    workspace::Ptr{Cfloat}
    train::Cint
    index::Cint
    net::network
end

struct image
    w::Cint
    h::Cint
    c::Cint
    data::Ptr{Cfloat}
end

struct box
    x::Cfloat
    y::Cfloat
    w::Cfloat
    h::Cfloat
end

struct detection
    bbox::box
    classes::Cint
    prob::Ptr{Cfloat}
    mask::Ptr{Cfloat}
    objectness::Cfloat
    sort_class::Cint
end

struct matrix
    rows::Cint
    cols::Cint
    vals::Ptr{Ptr{Cfloat}}
end

struct data
    w::Cint
    h::Cint
    X::matrix
    y::matrix
    shallow::Cint
    num_boxes::Ptr{Cint}
    boxes::Ptr{Ptr{box}}
end

@cenum(data_type,
    CLASSIFICATION_DATA = 0,
    DETECTION_DATA = 1,
    CAPTCHA_DATA = 2,
    REGION_DATA = 3,
    IMAGE_DATA = 4,
    COMPARE_DATA = 5,
    WRITING_DATA = 6,
    SWAG_DATA = 7,
    TAG_DATA = 8,
    OLD_CLASSIFICATION_DATA = 9,
    STUDY_DATA = 10,
    DET_DATA = 11,
    SUPER_DATA = 12,
    LETTERBOX_DATA = 13,
    REGRESSION_DATA = 14,
    SEGMENTATION_DATA = 15,
    INSTANCE_DATA = 16,
    ISEG_DATA = 17,
)

struct load_args
    threads::Cint
    paths::Ptr{Cstring}
    path::Cstring
    n::Cint
    m::Cint
    labels::Ptr{Cstring}
    h::Cint
    w::Cint
    c::Cint
    out_w::Cint
    out_h::Cint
    nh::Cint
    nw::Cint
    num_boxes::Cint
    min::Cint
    max::Cint
    size::Cint
    classes::Cint
    background::Cint
    scale::Cint
    center::Cint
    coords::Cint
    mini_batch::Cint
    track::Cint
    augment_speed::Cint
    jitter::Cfloat
    flip::Cint
    angle::Cfloat
    aspect::Cfloat
    saturation::Cfloat
    exposure::Cfloat
    hue::Cfloat
    d::Ptr{data}
    im::Ptr{image}
    resized::Ptr{image}
    type::data_type
    hierarchy::Ptr{tree}
end

struct metadata
    classes::Cint
    names::Ptr{Cstring}
end

@cenum(IMTYPE,
    PNG = 0,
    BMP = 1,
    TGA = 2,
    JPG = 3,
)
@cenum(BINARY_ACTIVATION,
    MULT = 0,
    ADD = 1,
    SUB = 2,
    DIV = 3,
)

struct update_args
    batch::Cint
    learning_rate::Cfloat
    momentum::Cfloat
    decay::Cfloat
    adam::Cint
    B1::Cfloat
    B2::Cfloat
    eps::Cfloat
    t::Cint
end

struct box_label
    id::Cint
    x::Cfloat
    y::Cfloat
    w::Cfloat
    h::Cfloat
    left::Cfloat
    right::Cfloat
    top::Cfloat
    bottom::Cfloat
end