text:
# # Dry Atmosphere GCM diagnostics
#
# This file computes selected diagnostics for the GCM and outputs them
# on the spherical interpolated diagnostic grid.
#
# Use it by calling `Diagnostics.setup_atmos_default_diagnostics()`.
#
# TODO:
# - enable zonal means and calculation of covariances using those means
#   - ds.T_zm = mean(.*1., ds.T; dims = 3)
#   - ds.u_zm = mean((ds.u); dims = 3 )
#   - v_zm = mean(ds.v; dims = 3)
#   - w_zm = mean(ds.w; dims = 3)
#   - ds.uvcovariance = (ds.u .- ds.u_zm) * (ds.v .- v_zm)
#   - ds.vTcovariance = (ds.v .- v_zm) * (ds.T .- ds.T_zm)
# - add more variables, including horiz streamfunction from laplacial of vorticity (LN)
# - density weighting
# - maybe change thermo/dyn separation to local/nonlocal vars?

import CUDA
using LinearAlgebra
using Printf
using Statistics

using ..Atmos
using ..Atmos: recover_thermo_state
using ..DGMethods.NumericalFluxes
using ..TurbulenceClosures: turbulence_tensors

include("diagnostic_fields.jl")
include("vorticity_balancelaw.jl")

mutable struct VorticityBLState <: DiagnosticsGroupParams
    bl::Union{Nothing, VorticityModel}
    dg::Union{Nothing, DGModel}
    state::Union{Nothing, MPIStateArray}
    dQ::Union{Nothing, MPIStateArray}

    VorticityBLState() = new(nothing, nothing, nothing, nothing)
end

"""
    setup_atmos_default_diagnostics(
        ::AtmosGCMConfigType,
        interval::String,
        out_prefix::String;
        writer::AbstractWriter,
        interpol = nothing,
    )

Create the "AtmosGCMDefault" `DiagnosticsGroup` which contains the following
diagnostic variables:
- u: zonal wind
- v: meridional wind
- w: vertical wind
- rho: air density
- temp: air temperature
- pres: air pressure
- thd: dry potential temperature
- et: total specific energy
- ei: specific internal energy
- ht: specific enthalpy based on total energy
- hi: specific enthalpy based on internal energy
- vort: vertical component of relative vorticity
- vort2: vertical component of relative vorticity from DGModel kernels via a mini balance law

When an `EquilMoist` moisture model is used, the following diagnostic
variables are also output:

- qt: mass fraction of total water in air
- ql: mass fraction of liquid water in air
- qv: mass fraction of water vapor in air
- qi: mass fraction of ice in air
- thv: virtual potential temperature
- thl: liquid-ice potential temperature

All these variables are output with `lat`, `long`, and `level` dimensions
of an interpolated grid (`interpol` _must_ be specified) as well as a
(unlimited) `time` dimension at the specified `interval`.
"""
function setup_atmos_default_diagnostics(
    ::AtmosGCMConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert !isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosGCMDefault",
        Diagnostics.atmos_gcm_default_init,
        Diagnostics.atmos_gcm_default_fini,
        Diagnostics.atmos_gcm_default_collect,
        interval,
        out_prefix,
        writer,
        interpol,
        VorticityBLState(),
    )
end

# Declare all (3D) variables for this diagnostics group
function vars_atmos_gcm_default_simple_3d(atmos::AtmosModel, FT)
    @vars begin
        u::FT
        v::FT
        w::FT
        rho::FT
        temp::FT
        pres::FT
        thd::FT                 # θ_dry
        et::FT                  # e_tot
        ei::FT                  # e_int
        ht::FT
        hi::FT
        vort::FT                # Ω₃
        vort2::FT               # Ω_bl₃

        moisture::vars_atmos_gcm_default_simple_3d(atmos.moisture, FT)
    end
end
vars_atmos_gcm_default_simple_3d(::MoistureModel, FT) = @vars()
function vars_atmos_gcm_default_simple_3d(m::EquilMoist, FT)
    @vars begin
        qt::FT                  # q_tot
        ql::FT                  # q_liq
        qv::FT                  # q_vap
        qi::FT                  # q_ice
        thv::FT                 # θ_vir
        thl::FT                 # θ_liq

    end
end
num_atmos_gcm_default_simple_3d_vars(m, FT) =
    varsize(vars_atmos_gcm_default_simple_3d(m, FT))
atmos_gcm_default_simple_3d_vars(m, array) =
    Vars{vars_atmos_gcm_default_simple_3d(m, eltype(array))}(array)

# Collect all (3D) variables for this diagnostics group
function atmos_gcm_default_simple_3d_vars!(
    atmos::AtmosModel,
    state_prognostic,
    thermo,
    dyni,
    dyn_bli,
    vars,
)
    vars.u = state_prognostic.ρu[1] / state_prognostic.ρ
    vars.v = state_prognostic.ρu[2] / state_prognostic.ρ
    vars.w = state_prognostic.ρu[3] / state_prognostic.ρ
    vars.rho = state_prognostic.ρ
    vars.et = state_prognostic.ρe / state_prognostic.ρ

    vars.temp = thermo.temp
    vars.pres = thermo.pres
    vars.thd = thermo.θ_dry
    vars.ei = thermo.e_int
    vars.ht = thermo.h_tot
    vars.hi = thermo.h_int

    vars.vort = dyni.Ω₃

    vars.vort2 = dyn_bli.Ω_bl₃

    atmos_gcm_default_simple_3d_vars!(
        atmos.moisture,
        state_prognostic,
        thermo,
        vars,
    )

    return nothing
end
function atmos_gcm_default_simple_3d_vars!(
    ::MoistureModel,
    state_prognostic,
    thermo,
    vars,
)
    return nothing
end
function atmos_gcm_default_simple_3d_vars!(
    moist::EquilMoist,
    state_prognostic,
    thermo,
    vars,
)
    vars.moisture.qt = state_prognostic.moisture.ρq_tot / state_prognostic.ρ
    vars.moisture.ql = thermo.moisture.q_liq
    vars.moisture.qv = thermo.moisture.q_vap
    vars.moisture.qi = thermo.moisture.q_ice
    vars.moisture.thv = thermo.moisture.θ_vir
    vars.moisture.thl = thermo.moisture.θ_liq_ice

    return nothing
end

# Dynamic variables
function vars_dyn(FT)
    @vars begin
        Ω₁::FT
        Ω₂::FT
        Ω₃::FT
    end
end
dyn_vars(array) = Vars{vars_dyn(eltype(array))}(array)

function vars_dyn_bl(FT)
    @vars begin
        Ω_bl₁::FT
        Ω_bl₂::FT
        Ω_bl₃::FT
    end
end
dyn_bl_vars(array) = Vars{vars_dyn_bl(eltype(array))}(array)

"""
    atmos_gcm_default_init(dgngrp, currtime)

Initialize the GCM default diagnostics group, establishing the output file's
dimensions and variables.
"""
function atmos_gcm_default_init(dgngrp::DiagnosticsGroup, currtime)
    dg = Settings.dg
    grid = dg.grid
    atmos = dg.balance_law
    FT = eltype(Settings.Q)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if !(dgngrp.interpol isa InterpolationCubedSphere)
        @warn """
            Diagnostics ($dgngrp.name): currently requires `InterpolationCubedSphere`!
            """
        return nothing
    end

    # set up the vorticity mini balance law
    vort_state = dgngrp.params
    vort_state.bl = VorticityModel()
    vort_state.dg = DGModel(
        vort_state.bl,
        grid,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )
    vort_state.state = init_ode_state(vort_state.dg, FT(0))
    vort_state.dQ = similar(
        vort_state.state;
        vars = @vars(Ω_bl::SVector{3, FT}),
        nstate = 3,
    )

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # adjust the level dimension for `planet_radius`
        level_val = dims["level"]
        dims["level"] = (
            level_val[1] .- FT(planet_radius(Settings.param_set)),
            level_val[2],
        )

        # set up the variables we're going to be writing
        vars = OrderedDict()
        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_gcm_default_simple_3d(atmos, FT)),
        )
        for varname in varnames
            var = Variables[varname]
            vars[varname] = (tuple(collect(keys(dims))...), FT, var.attrib)
        end

        # create the output file
        dprefix = @sprintf(
            "%s_%s_%s",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)

        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

"""
    atmos_gcm_default_collect(dgngrp, currtime)

    Master function that performs a global grid traversal to compute various
    diagnostics using the above functions.
"""
function atmos_gcm_default_collect(dgngrp::DiagnosticsGroup, currtime)
    interpol = dgngrp.interpol
    if !(interpol isa InterpolationCubedSphere)
        @warn """
            Diagnostics ($dgngrp.name): currently requires `InterpolationCubedSphere`!
            """
        return nothing
    end
    vort_state = dgngrp.params

    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    atmos = dg.balance_law
    grid = dg.grid
    grid_info = basic_grid_info(dg)
    topl_info = basic_topology_info(grid.topology)
    Nqk = grid_info.Nqk
    Nqh = grid_info.Nqh
    npoints = prod(grid_info.Nq)
    nrealelem = topl_info.nrealelem
    nvertelem = topl_info.nvertelem
    nhorzelem = topl_info.nhorzrealelem

    # get needed arrays onto the CPU
    device = array_device(Q)
    if device isa CPU
        ArrayType = Array
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
    else
        ArrayType = CUDA.CuArray
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
    end
    FT = eltype(state_data)

    # TODO: can this be done in one pass?
    #
    # Non-local vars, e.g. relative vorticity
    vgrad = VectorGradients(dg, Q)
    vort = Vorticity(dg, vgrad)

    # run the vorticity mini balance law
    ix_ρu = varsindex(vars(Q), :ρu)
    ix_ρ = varsindex(vars(Q), :ρ)
    ρ = Q.data[:, ix_ρ, :]
    u = Q.data[:, ix_ρu, :] ./ ρ

    vort_state.dg.state_auxiliary.data .= u
    vort_state.dg(vort_state.dQ, vort_state.state, nothing, FT(0))

    # Compute thermo variables element-wise
    thermo_array = Array{FT}(undef, npoints, num_thermo(atmos, FT), nrealelem)
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        aux = extract_state(dg, aux_data, ijk, e, Auxiliary())

        thermo = thermo_vars(atmos, view(thermo_array, ijk, :, e))
        compute_thermo!(atmos, state, aux, thermo)
    end

    # Interpolate the state, thermo, dgdiags and dyn vars to sphere (u and vorticity
    # need projection to zonal, merid). All this may happen on the GPU.
    istate =
        ArrayType{FT}(undef, interpol.Npl, number_states(atmos, Prognostic()))
    interpolate_local!(interpol, Q.realdata, istate)

    ithermo = ArrayType{FT}(undef, interpol.Npl, num_thermo(atmos, FT))
    interpolate_local!(interpol, ArrayType(thermo_array), ithermo)

    idyn = ArrayType{FT}(undef, interpol.Npl, size(vort.data, 2))
    interpolate_local!(interpol, vort.data, idyn)

    idyn_bl = ArrayType{FT}(undef, interpol.Npl, size(vort_state.dQ.data, 2))
    interpolate_local!(interpol, vort_state.dQ.data, idyn_bl)

    # TODO: get indices here without hard-coding them
    _ρu, _ρv, _ρw = 2, 3, 4
    project_cubed_sphere!(interpol, istate, (_ρu, _ρv, _ρw))
    _Ω₁, _Ω₂, _Ω₃ = 1, 2, 3
    project_cubed_sphere!(interpol, idyn, (_Ω₁, _Ω₂, _Ω₃))
    project_cubed_sphere!(interpol, idyn_bl, (_Ω₁, _Ω₂, _Ω₃))

    # FIXME: accumulating to rank 0 is not scalable
    all_state_data = accumulate_interpolated_data(mpicomm, interpol, istate)
    all_thermo_data = accumulate_interpolated_data(mpicomm, interpol, ithermo)
    all_dyn_data = accumulate_interpolated_data(mpicomm, interpol, idyn)
    all_dyn_bl_data = accumulate_interpolated_data(mpicomm, interpol, idyn_bl)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # set up the array for the diagnostic variables based on the interpolated grid
        nlong = length(dims["long"][1])
        nlat = length(dims["lat"][1])
        nlevel = length(dims["level"][1])

        simple_3d_vars_array = Array{FT}(
            undef,
            nlong,
            nlat,
            nlevel,
            num_atmos_gcm_default_simple_3d_vars(atmos, FT),
        )

        @traverse_interpolated_grid nlong nlat nlevel begin
            statei = Vars{vars_state(atmos, Prognostic(), FT)}(view(
                all_state_data,
                lo,
                la,
                le,
                :,
            ))
            thermoi = thermo_vars(atmos, view(all_thermo_data, lo, la, le, :))
            dyni = dyn_vars(view(all_dyn_data, lo, la, le, :))
            dyn_bli = dyn_bl_vars(view(all_dyn_bl_data, lo, la, le, :))
            simple_3d_vars = atmos_gcm_default_simple_3d_vars(
                atmos,
                view(simple_3d_vars_array, lo, la, le, :),
            )

            atmos_gcm_default_simple_3d_vars!(
                atmos,
                statei,
                thermoi,
                dyni,
                dyn_bli,
                simple_3d_vars,
            )
        end

        # assemble the diagnostics for writing
        varvals = OrderedDict()
        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_gcm_default_simple_3d(atmos, FT)),
        )
        for (vari, varname) in enumerate(varnames)
            varvals[varname] = simple_3d_vars_array[:, :, :, vari]
        end

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_gcm_default_fini(dgngrp::DiagnosticsGroup, currtime) end

max_forks_repo_name: jm-c/CLIMA


text:
__precompile__(true)
"""
    NCEI

The `NCEI` module provides an interface to the National Centers of Environmental
Information API. The Climate Data Online API is fully supported.

Please see the documenation for more information.
"""
module NCEI

using Dates: Day, Date, DateTime, TimeType, today, Year
using DataFrames: DataFrame
using HTTP: request
using LazyJSON: value # Needs to be run first in MacOS for some reason
using LazyJSON.PropertyDicts: get

import Base: names, parse

# Code from Reexport.jl
"""
    @reexport(ex::Expr)
Reexports the package (e.g., `@reexport using PkgA`).
"""
macro reexport(ex)
    isa(ex, Expr) && (ex.head == :module ||
                      ex.head == :using ||
                      (ex.head == :toplevel &&
                       all(e->isa(e, Expr) && e.head == :using, ex.args))) ||
        error("@reexport: syntax error")

    if ex.head == :module
        modules = convert(Vector{Any}, [ex.args[2]])
        ex = Expr(:toplevel, ex, :(using .$(ex.args[2])))
    elseif ex.head == :using && all(e->isa(e, Symbol), ex.args)
        modules = convert(Vector{Any}, [ex.args[end]])
    else
        modules = convert(Vector{Any}, [e.args[end] for e in ex.args])
    end

    esc(Expr(:toplevel, ex,
             [:(eval(Expr(:export, names($mod)...))) for mod in modules]...))
end

@reexport using DataFrames

include.(["Endpoint.jl", "API.jl"])

export
    Date,
    DateTime,
    cdo_datasets,
    cdo_datacategories,
    cdo_datatypes,
    cdo_locationcategories,
    cdo_locations,
    cdo_stations,
    cdo_data
end

ext: jl


text:
using SpineOpt

run_spineopt(ARGS...)


text:
function treepositions(hc::Hclust, useheight::Bool)
    order = StatsBase.indexmap(hc.order)
    nodepos = Dict(-i => (float(order[i]), 0.0) for i in hc.order)

    xs = Array{Float64}(undef, 4, size(hc.merges, 1))
    ys = Array{Float64}(undef, 4, size(hc.merges, 1))

    for i in 1:size(hc.merges, 1)
        x1, y1 = nodepos[hc.merges[i, 1]]
        x2, y2 = nodepos[hc.merges[i, 2]]

        xpos = (x1 + x2) / 2
        useheight ? h = hc.heights[i] : h = 1
        ypos = max(y1, y2) + h

        nodepos[i] = (xpos, ypos)
        xs[:, i] .= [x1, x1, x2, x2]
        ys[:, i] .= [y1, ypos, ypos, y2]
    end

    return xs, ys
end

@recipe function f(hc::Hclust; useheight=true)
    typeof(useheight) <: Bool || error("'useheight' argument must be true or false")

    legend := false
    xforeground_color_axis := :white
    xgrid := false
    xlims := (0.5, length(hc.order) + 0.5)

    linecolor --> :black
    xticks --> (1:nnodes(hc), string.(1:nnodes(hc))[hc.order])
    ylims --> (0, Inf)
    yshowaxis --> useheight

    treepositions(hc, useheight)
end


text:
using Formatting, FFTW, Jets, JetPackWaveFD, LinearAlgebra, SpecialFunctions, Statistics, Test, WaveFD

function make_op(modeltype, interpmethod, fs; v₀=1500, ϵ₀=0.2, η₀=0.4, comptype = Float32, st = 0.0, wavelet = WaveletCausalRicker(f=5.0))
    nsponge = 10
    pad = 20
    nx = 60+2*pad
    nz = 50+2*pad
    nt = 101
    dx,dz,dt = 25.0,25.0,0.002
    xmin,zmin,tmin = -pad*dx,-pad*dz,0.0
    xmax,zmax,tmax = xmin+dx*(nx-pad-1),zmin+dz*(nz-pad-1),tmin+dt*(nt-1)

    sx = dx*div(nx,2)
    sz = dz
    rx = dx*[pad:nx-pad-1;]
    rz = 2*dz*ones(length(rx))

    local ϵ,η
    if modeltype  == "v"
        ϵ = zeros(Float32,nz,nx)
        η = zeros(Float32,nz,nx)
    else
        ϵ = Array{Float32}(undef,0,0)
        η = Array{Float32}(undef,0,0)
    end

    b = ones(Float32,nz,nx)

    F = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(b = b,
        f = 0.85, ϵ = ϵ, η = η,
        sx = sx, sz = sz, st = st, rx = rx, rz = rz, dx = dx, dz = dz, z0 = zmin, x0 = xmin,
        dtrec = dt, dtmod = dt, ntrec = nt, nsponge = nsponge, comptype = comptype, compscale = 1e-4,
        freqQ = 5.0, qMin = 0.1, qInterior = 100.0, wavelet = wavelet, freesurface = fs,
        reportinterval = 0, interpmethod = interpmethod, isinterior = false)

    m₀ = zeros(domain(F))
    m₀[:,:,modelindex(F,"v")] .= v₀
    if modeltype == "vϵη"
        m₀[:,:,modelindex(F,"ϵ")] .= ϵ₀
        m₀[:,:,modelindex(F,"η")] .= η₀
    end

    m₀,F
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- close" begin
    m₀, F = make_op("vϵη", :hicks, false)
    d = F * m₀

    @test isfile(state(F).srcfieldfile*"-pspace")
    @test isfile(state(F).srcfieldfile*"-mspace")
    close(F)
    @test !(isfile(state(F).srcfieldfile*"-pspace"))
    @test !(isfile(state(F).srcfieldfile*"-mspace"))
    d = F * m₀
    @test isfile(state(F).srcfieldfile*"-pspace")
    @test isfile(state(F).srcfieldfile*"-mspace")
    close(F)
    @test !(isfile(state(F).srcfieldfile*"-pspace"))
    @test !(isfile(state(F).srcfieldfile*"-mspace"))

    close(F)
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- perfstat" begin
    m₀, F = make_op("vϵη", :hicks, false)
    d = F * m₀

    s = Jets.perfstat(F)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test s["%imaging"] ≈ 0.0

    J = jacobian!(F, m₀)
    d = J * m₀
    s = Jets.perfstat(J)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test isa(s["%imaging"], Float64)

    m = J' * d
    s = Jets.perfstat(J)
    @test isa(s["MCells/s"], Float64)
    @test isa(s["%io"], Float64)
    @test isa(s["%inject/extract"], Float64)
    @test isa(s["%imaging"], Float64)

    close(F)
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- non-causal source injection" begin
    wavelet = WaveletRicker(;f=5.0)
    mₒ, F₁ = make_op(:hicks, "v", false; wavelet, st=-1.0)
    mₒ, F₂ = make_op(:hicks, "v", false; wavelet, st=-2.0)

    d₁ = F₁ * mₒ
    d₂ = F₂ * mₒ
    @test d₁ ≈ d₂
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- linearization, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear, :hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    d = F * m₀

    mmask = zeros(domain(F))
    mmask[:,:,modelindex(F,"v")] .= 1
    if modeltype == "vϵη"
        mmask[:,:,modelindex(F,"ϵ")] .= 0.001
        mmask[:,:,modelindex(F,"η")] .= 0.001
    end

    μobs, μexp = linearization_test(F, m₀,
        mmask = mmask,
        μ = 100*sqrt.([1.0,1.0/2.0,1.0/4.0,1.0/8.0,1.0/16.0,1.0/32.0,1.0/64.0,1.0/128.0,1.0/256.0,1.0/512.0]))

    @show μobs, μexp
    @show minimum(abs, μobs - μexp)
    @test minimum(abs, μobs - μexp) < 0.01

    close(F)
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- dot product, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    d = F * m₀
    J = jacobian!(F, m₀)

    m = zeros(domain(J))
    m[:,:,modelindex(F,"v")] .= 100 .* (-1 .+ 2*rand(Float32,size(m)[1:2]))
    if modeltype == "vϵη"
        m[:,:,modelindex(F,"ϵ")] .= rand(Float32,size(m)[1:2])
        m[:,:,modelindex(F,"η")] .= rand(Float32,size(m)[1:2])
    end
    d = -1 .+ 2*rand(range(J))

    lhs, rhs = dot_product_test(J, m, d)
    err = abs(lhs - rhs) / abs(lhs + rhs)
    @show lhs,rhs,(lhs-rhs)/(lhs+rhs)
    @test lhs ≈ rhs rtol=1e-4
    @test err < 1e-4

    close(F)
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- linearity, model type is $modeltype, interpmethod=$interpmethod, fs=$fs" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks), fs in (true,false)
    m₀, F = make_op(modeltype, interpmethod, fs)
    d = F * m₀
    J = jacobian!(F, m₀)
    lhs,rhs = linearity_test(J)
    @test lhs ≈ rhs

    close(F)
end

# note the compression is exercised on the second pass of F * m₀
@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- serialization, modeltype=$modeltype, C=$C" for modeltype in ("vϵη", "v"), C in (Float32,UInt32)
    m₀, F = make_op(modeltype, :hicks, false, comptype = C)
    d₁ = F * m₀
    d₂ = F * m₀
    @test d₁ ≈ d₂

    close(F)
end

@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- srcillum, C=$C, modeltype=$modeltype" for C in (Float32, UInt32), modeltype in ("vϵη", "v")
    m₀, F = make_op(:hicks, modeltype, false, comptype = C)
    d₁ = F * m₀;
    J = jacobian(F, m₀);
    s1 = srcillum(F, m₀);
    s2 = srcillum(J);
    close(F)
    @test s1 ≈ s2
    @test maximum(s1) > 0
end

#=
@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:hicks,:linear)
    z,x,dz,dx,dtrec,dtmod,t,sz,sx,c = 0.5,1.0,0.02,0.02,0.004,0.0004,1.0,0.25,0.02,1.5
    nz,nx,ntrec = round(Int,z/dz)+1,round(Int,x/dx)+1,round(Int,t/dtrec)+1

    fmax = c/6/max(dx,dz)
    fdef = WaveletCausalRicker(f=fmax/2);

    function analytic_result()
        # non-negative frequency axis
        nt_pad = nextprod([2,3,5,7],4*ntrec)
        nw = div(nt_pad,2)+1
        w = range(0, stop=pi/dtrec, length=nw)

        # wavelet
        f = WaveFD.get(fdef, dtrec*collect(0:nt_pad-1))
        F = rfft(f)

        # Morse and Feshbach, 1953, p. 891
        u_ana = zeros(Float32,ntrec,nz,nx)
        U = zeros(Complex{Float64},nw)
        for ix = 1:nx, iz = 1:nz
            z = (iz-1)*dz
            x = (ix-1)*dx
            r = sqrt((z-sz)^2+(x-sx)^2)
            for iw = 2:nw
                if r*w[iw] > eps(Float32)
                    U[iw] = -im*pi*hankelh2(0,w[iw]/c*r)*F[iw]
                end
            end
            u_ana[:,iz,ix] = irfft(U,nt_pad)[1:ntrec]
            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100))
        end
        u_ana
    end
    write(stdout, "computing analytic result (be patient, can take several minutes)...\n")
    u_ana = analytic_result()
    write(stdout, "...done.\n")

    local ϵ,η
    if modeltype == "v"
        ϵ = zeros(Float32,nz,nx)
        η = zeros(Float32,nz,nx)
    else
        ϵ = Float32[]
        η = Float32[]
    end

    # Finite difference
    M = JopNlProp2DAcoVTIDenQ_DEO2_FDTD(
        b = (1.0f0/2.0f0)*ones(Float32,nz,nx),
        f = 0.75,
        ϵ = ϵ,
        η = η,
        sx = sx,
        sz = sz,
        rz = sz*ones(nx),
        rx = dx*collect(0:nx-1),
        interpmethod = interpmethod,
        dx = dx,
        dz = dz,
        dtrec = dtrec,
        dtmod = dtmod,
        ntrec = ntrec,
        wavelet = fdef,
        nsponge = 100,
        freqQ = 5.0,
        qMin = 0.1,
        qInterior = 100.0,
        freesurface = false,
        reportinterval = 0)

    m = zeros(domain(M))
    m[:,:,modelindex(M,"v")] .= c
    write(stdout, "computing numerical result (be patient)...\n")
    d_tmp = M*m
    write(stdout, "...done.\n")

    # Read FD field back from disk
    u_fd = zeros(Float32,ntrec,nz,nx)
    u_snap = zeros(Float32,nz,nx)
    u_snap_ginsu = zeros(Float32,size(state(M).ginsu)...)
    io = open("$(state(M).srcfieldfile)-pold")
    for it = 1:state(M).ntrec
        read!(io,u_snap_ginsu)
        super!(u_snap,state(M).ginsu,u_snap_ginsu)
        u_fd[it,:,:] = u_snap[:,:]
    end
    close(io)

    # compute Pearson correlation (https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient)
    xx = zeros(nz,nx)
    for ix = 1:nx, iz = 1:nz
        xx[iz,ix] = cor(u_fd[:,iz,ix],u_ana[:,iz,ix])
    end
    xx_mean = mean(xx)
    @test xx_mean > .9

    # measure geometric spreading
    amp_ana = zeros(nz,nx)
    amp_fd = zeros(nz,nx)
    for ix = 1:nx, iz = 1:nz
        amp_ana[iz,ix] = maximum(u_ana[:,iz,ix])
        amp_fd[iz,ix] = maximum(u_fd[:,iz,ix])
    end

    # don't trust the values at the injection point
    szi = round(Int,sz/dz)+1
    sxi = round(Int,sx/dx)+1
    amp_fd[szi,sxi] = 0.0
    amp_ana[szi,sxi] = 0.0

    # we only care about relative values
    amp_fd /= maximum(amp_fd)
    amp_ana /= maximum(amp_ana)
    @test cor(amp_fd[:],amp_ana[:]) > .99

    close(M)
end
=#

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
function check_error(err, error_msg)
    if err != 0
        error(error_msg * " with error: $(err)")
    end
end

const FREE_FONT_LIBRARY = FT_Library[C_NULL]

function ft_init()
    FREE_FONT_LIBRARY[1] != C_NULL && error("Freetype already initalized. init() called two times?")
    err = FT_Init_FreeType(FREE_FONT_LIBRARY)
    return err == 0
end

function ft_done()
    FREE_FONT_LIBRARY[1] == C_NULL && error("Library == CNULL. FreeTypeAbstraction.done() called before init(), or done called two times?")
    err = FT_Done_FreeType(FREE_FONT_LIBRARY[1])
    FREE_FONT_LIBRARY[1] = C_NULL
    return err == 0
end

function newface(facename, faceindex::Real=0, ftlib=FREE_FONT_LIBRARY)
    face = Ref{FT_Face}()
    err = FT_New_Face(ftlib[1], facename, Int32(faceindex), face)
    check_error(err, "Couldn't load font $facename")
    return face[]
end

struct FontExtent{T}
    vertical_bearing::Vec{2, T}
    horizontal_bearing::Vec{2, T}

    advance::Vec{2, T}
    scale::Vec{2, T}
end

BroadcastStyle(::Type{<: FontExtent}) = Style{FontExtent}()
BroadcastStyle(::Style{FontExtent}, x) = Style{FontExtent}()
BroadcastStyle(x, ::Style{FontExtent}) = Style{FontExtent}()

function broadcasted(op::Function, f::FontExtent, scaling::StaticVector)
    return FontExtent(
        op.(f.vertical_bearing, scaling[1]),
        op.(f.horizontal_bearing, scaling[2]),
        op.(f.advance, scaling),
        op.(f.scale, scaling),
    )
end

function broadcasted(op::Function, f::FontExtent)
    return FontExtent(
        op.(f.vertical_bearing),
        op.(f.horizontal_bearing),
        op.(f.advance),
        op.(f.scale),
    )
end

function broadcasted(op::Function, ::Type{T}, f::FontExtent) where T
    return FontExtent(
        map(x-> op(T, x), f.vertical_bearing),
        map(x-> op(T, x), f.horizontal_bearing),
        map(x-> op(T, x), f.advance),
        map(x-> op(T, x), f.scale),
    )
end

function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::T = 64.0) where T <: AbstractFloat
    return FontExtent(
        Vec{2, T}(fontmetric.vertBearingX, fontmetric.vertBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiBearingX, fontmetric.horiBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiAdvance, fontmetric.vertAdvance) ./ scale,
        Vec{2, T}(fontmetric.width, fontmetric.height) ./ scale
    )
end

function ==(x::FontExtent, y::FontExtent)
    return (x.vertical_bearing == y.vertical_bearing &&
            x.horizontal_bearing == y.horizontal_bearing &&
            x.advance == y.advance &&
            x.scale == y.scale)
end

function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::Integer)
    return FontExtent(
        div.(Vec{2, Int}(fontmetric.vertBearingX, fontmetric.vertBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiBearingX, fontmetric.horiBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiAdvance, fontmetric.vertAdvance), scale),
        div.(Vec{2, Int}(fontmetric.width, fontmetric.height), scale)
    )
end

function bearing(extent::FontExtent)
    return Vec2f0(extent.horizontal_bearing[1],
                  -(extent.scale[2] - extent.horizontal_bearing[2]))
end

function safe_free(face)
    ptr = getfield(face, :ft_ptr)
    if ptr != C_NULL && FREE_FONT_LIBRARY[1] != C_NULL
        FT_Done_Face(face)
    end
end

function boundingbox(extent::FontExtent)
    mini = bearing(extent)
    return Rect2D(mini, Vec2f0(extent.scale))
end

mutable struct FTFont
    ft_ptr::FreeType.FT_Face
    current_pixelsize::Base.RefValue{Int}
    use_cache::Bool
    cache::Dict{Tuple{Int, Char}, FontExtent{Float32}}
    function FTFont(ft_ptr::FreeType.FT_Face, pixel_size::Int=64, use_cache::Bool=true)
        cache = Dict{Tuple{Int, Char}, FontExtent{Float32}}()
        face = new(ft_ptr, Ref(pixel_size), use_cache, cache)
        finalizer(safe_free, face)
        FT_Set_Pixel_Sizes(face, pixel_size, 0);
        return face
    end
end

use_cache(face::FTFont) = getfield(face, :use_cache)
get_cache(face::FTFont) = getfield(face, :cache)

function FTFont(path::String)
    return FTFont(newface(path))
end

# C interop
function Base.cconvert(::Type{FreeType.FT_Face}, font::FTFont)
    return font
end

function Base.unsafe_convert(::Type{FreeType.FT_Face}, font::FTFont)
    return getfield(font, :ft_ptr)
end

function Base.propertynames(font::FTFont)
    return fieldnames(FreeType.FT_FaceRec)
end

function Base.getproperty(font::FTFont, fieldname::Symbol)
    fontrect = unsafe_load(getfield(font, :ft_ptr))
    field = getfield(fontrect, fieldname)
    if field isa Ptr{FT_String}
        field == C_NULL && return ""
        return unsafe_string(field)
    # Some fields segfault with unsafe_load...Lets find out which another day :D
    elseif field isa Ptr{FreeType.LibFreeType.FT_GlyphSlotRec}
        return unsafe_load(field)
    else
        return field
    end
end

get_pixelsize(face::FTFont) = getfield(face, :current_pixelsize)[]

function set_pixelsize(face::FTFont, size::Integer)
    get_pixelsize(face) == size && return size
    err = FT_Set_Pixel_Sizes(face, size, size)
    check_error(err, "Couldn't set pixelsize")
    getfield(face, :current_pixelsize)[] = size
    return size
end

function kerning(c1::Char, c2::Char, face::FTFont)
    i1 = FT_Get_Char_Index(face, c1)
    i2 = FT_Get_Char_Index(face, c2)
    kerning2d = Ref{FreeType.FT_Vector}()
    err = FT_Get_Kerning(face, i1, i2, FreeType.FT_KERNING_DEFAULT, kerning2d)
    # Can error if font has no kerning! Since that's somewhat expected, we just return 0
    err != 0 && return Vec2f0(0)
    # 64 since metrics are in 1/64 units (units to 26.6 fractional pixels)
    divisor = 64
    return Vec2f0(kerning2d[].x / divisor, kerning2d[].y / divisor)
end

function loadchar(face::FTFont, c::Char)
    err = FT_Load_Char(face, c, FT_LOAD_RENDER)
    check_error(err, "Could not load char to render.")
end

function get_extent(face::FTFont, char::Char)
    if use_cache(face)
        get!(get_cache(face), (get_pixelsize(face), char)) do
            return internal_get_extent(face, char)
        end
    else
        return internal_get_extent(face, char)
    end
end

function internal_get_extent(face::FTFont, char::Char)
    err = FT_Load_Char(face, char, FT_LOAD_DEFAULT)
    check_error(err, "Could not load char to get extend.")
    metrics = face.glyph.metrics
    # 64 since metrics are in 1/64 units (units to 26.6 fractional pixels)
    return FontExtent(metrics, Float32(64))
end


text:
include("SA.jl")
using HDF5

# n_lst = collect(8000:100:15000)
n_lst = collect(10000)
e_lst = 3*n_lst
reps = 5
t = zeros(length(n_lst),reps)

dt = replace(string(now()),"T","_")
dt = replace(dt,":","_")
dt = replace(dt,"-","_")
fname = "emp_results"
folder0 = string("empirical_results_",dt)
folder = ""
# junk run to get everyhting in memory, not sure how to fix this yet.
SA.optimize()
# now start the real testing
for i = 1:length(n_lst)

    ARGS = Array(UTF8String,7)
    ARGS[1] = "--simple_graph"
    ARGS[2] = @sprintf("-n=%d",n_lst[i])
    ARGS[3] = @sprintf("-e=%d",e_lst[i])
    ARGS[4] = "-i=100"
    ARGS[5] = "-r=0.5"
    ARGS[7] = "--no_console_output"
    for j = 1:reps
        folder = string(folder0,"__",i,"_",j)
        ARGS[6] = @sprintf("--out_folder=%s",folder)
        @printf("Run %d trial %d, n=%d, e=%d\n",i,j,n_lst[i],e_lst[i])
        t[i,j] = SA.optimize(ARGS)
    end
    h5open(string(joinpath(folder,string(fname,"_data")),".h5"), "w") do file
    write(file, "t", t)  # alternatively, say "@write file A"
    write(file, "n_lst", n_lst)
    write(file, "e_lst", e_lst)
    write(file, "ARGS", ARGS)
    end
end
PyPlot.clf()
PyPlot.plot(n_lst,t,"b*")
PyPlot.title("Runtime vs. n")
PyPlot.xlabel("n")
PyPlot.ylabel("Runtime")
fname = string(fname,"_runtime",".pdf")
PyPlot.savefig(joinpath(folder,fname))
# println(t)


text:
M7 = Posit{7,0}
M7_1 = Posit{7,1}
M7_2 = Posit{7,2}

@test M7_1(0b0000000)   == M7_1(0.0)
@test M7_1(0b0000001)   == M7_1(0.0009765625)
@test M7_1(0b000001_0)  == M7_1(0.00390625)
@test M7_1(0b000001_1)  == M7_1(0.0078125)
@test M7_1(0b00001_0_0) == M7_1(0.015625)
@test M7_1(0b00001_0_1) == M7_1(0.0234375)
@test M7_1(0b00001_1_0) == M7_1(0.03125)
@test M7_1(0b00001_1_1) == M7_1(0.046875)
@test M7_1(0b0001_0_00) == M7_1(0.0625)
@test M7_1(0b0001_0_01) == M7_1(0.078125)
@test M7_1(0b0001_0_10) == M7_1(0.09375)
@test M7_1(0b0001_0_11) == M7_1(0.109375)
@test M7_1(0b0001_1_00) == M7_1(0.125)
@test M7_1(0b0001_1_01) == M7_1(0.15625)
@test M7_1(0b0001_1_10) == M7_1(0.1875)
@test M7_1(0b0001_1_11) == M7_1(0.21875)
@test M7_1(0b001_0_000) == M7_1(0.25)
@test M7_1(0b001_0_001) == M7_1(0.28125)
@test M7_1(0b001_0_010) == M7_1(0.3125)
@test M7_1(0b001_0_011) == M7_1(0.34375)
@test M7_1(0b001_0_100) == M7_1(0.375)
@test M7_1(0b001_0_101) == M7_1(0.40625)
@test M7_1(0b001_0_110) == M7_1(0.4375)
@test M7_1(0b001_0_111) == M7_1(0.46875)
@test M7_1(0b001_1_000) == M7_1(0.5)
@test M7_1(0b001_1_001) == M7_1(0.5625)
@test M7_1(0b001_1_010) == M7_1(0.625)
@test M7_1(0b001_1_011) == M7_1(0.6875)
@test M7_1(0b001_1_100) == M7_1(0.75)
@test M7_1(0b001_1_101) == M7_1(0.8125)
@test M7_1(0b001_1_110) == M7_1(0.875)
@test M7_1(0b001_1_111) == M7_1(0.9375)
@test M7_1(0b010_0_000) == M7_1(1.0)
@test M7_1(0b010_0_001) == M7_1(1.125)
@test M7_1(0b010_0_010) == M7_1(1.25)
@test M7_1(0b010_0_011) == M7_1(1.375)
@test M7_1(0b010_0_100) == M7_1(1.5)
@test M7_1(0b010_0_101) == M7_1(1.625)
@test M7_1(0b010_0_110) == M7_1(1.75)
@test M7_1(0b010_0_111) == M7_1(1.875)
@test M7_1(0b010_1_000) == M7_1(2.0)
@test M7_1(0b010_1_001) == M7_1(2.25)
@test M7_1(0b010_1_010) == M7_1(2.5)
@test M7_1(0b010_1_011) == M7_1(2.75)
@test M7_1(0b010_1_100) == M7_1(3.0)
@test M7_1(0b010_1_101) == M7_1(3.25)
@test M7_1(0b010_1_110) == M7_1(3.5)
@test M7_1(0b010_1_111) == M7_1(3.75)
@test M7_1(0b0110_0_00) == M7_1(4.0)
@test M7_1(0b0110_0_01) == M7_1(5.0)
@test M7_1(0b0110_0_10) == M7_1(6.0)
@test M7_1(0b0110_0_11) == M7_1(7.0)
@test M7_1(0b0110_1_00) == M7_1(8.0)
@test M7_1(0b0110_1_01) == M7_1(10.0)
@test M7_1(0b0110_1_10) == M7_1(12.0)
@test M7_1(0b0110_1_11) == M7_1(14.0)
@test M7_1(0b01110_0_0) == M7_1(16.0)
@test M7_1(0b01110_0_1) == M7_1(24.0)
@test M7_1(0b01110_1_0) == M7_1(32.0)
@test M7_1(0b01110_1_1) == M7_1(48.0)
@test M7_1(0b011110_0)  == M7_1(64.0)
@test M7_1(0b011110_1)  == M7_1(128.0)
@test M7_1(0b0111110)   == M7_1(256.0)
@test M7_1(0b0111111)   == M7_1(1024.0)

@test M7_2(0b0000000)   == M7_2(0.0)
@test M7_2(0b0000001)   == M7_2(0.00000095367431640625)
@test M7_2(0b000001_0)  == M7_2(0.0000152587890625)
@test M7_2(0b000001_1)  == M7_2(0.00006103515625)
@test M7_2(0b00001_00)  == M7_2(0.000244140625)
@test M7_2(0b00001_01)  == M7_2(0.00048828125)
@test M7_2(0b00001_10)  == M7_2(0.0009765625)
@test M7_2(0b00001_11)  == M7_2(0.001953125)
@test M7_2(0b0001_00_0) == M7_2(0.00390625)
@test M7_2(0b0001_00_1) == M7_2(0.005859375)
@test M7_2(0b0001_01_0) == M7_2(0.0078125)
@test M7_2(0b0001_01_1) == M7_2(0.01171875)
@test M7_2(0b0001_10_0) == M7_2(0.015625)
@test M7_2(0b0001_10_1) == M7_2(0.0234375)
@test M7_2(0b0001_11_0) == M7_2(0.03125)
@test M7_2(0b0001_11_1) == M7_2(0.046875)
@test M7_2(0b001_00_00) == M7_2(0.0625)
@test M7_2(0b001_00_01) == M7_2(0.078125)
@test M7_2(0b001_00_10) == M7_2(0.09375)
@test M7_2(0b001_00_11) == M7_2(0.109375)
@test M7_2(0b001_01_00) == M7_2(0.125)
@test M7_2(0b001_01_01) == M7_2(0.15625)
@test M7_2(0b001_01_10) == M7_2(0.1875)
@test M7_2(0b001_01_11) == M7_2(0.21875)
@test M7_2(0b001_10_00) == M7_2(0.25)
@test M7_2(0b001_10_01) == M7_2(0.3125)
@test M7_2(0b001_10_10) == M7_2(0.375)
@test M7_2(0b001_10_11) == M7_2(0.4375)
@test M7_2(0b001_11_00) == M7_2(0.5)
@test M7_2(0b001_11_01) == M7_2(0.625)
@test M7_2(0b001_11_10) == M7_2(0.75)
@test M7_2(0b001_11_11) == M7_2(0.875)
@test M7_2(0b010_00_00) == M7_2(1.0)
@test M7_2(0b010_00_01) == M7_2(1.25)
@test M7_2(0b010_00_10) == M7_2(1.5)
@test M7_2(0b010_00_11) == M7_2(1.75)
@test M7_2(0b010_01_00) == M7_2(2.0)
@test M7_2(0b010_01_01) == M7_2(2.5)
@test M7_2(0b010_01_10) == M7_2(3.0)
@test M7_2(0b010_01_11) == M7_2(3.5)
@test M7_2(0b010_10_00) == M7_2(4.0)
@test M7_2(0b010_10_01) == M7_2(5.0)
@test M7_2(0b010_10_10) == M7_2(6.0)
@test M7_2(0b010_10_11) == M7_2(7.0)
@test M7_2(0b010_11_00) == M7_2(8.0)
@test M7_2(0b010_11_01) == M7_2(10.0)
@test M7_2(0b010_11_10) == M7_2(12.0)
@test M7_2(0b010_11_11) == M7_2(14.0)
@test M7_2(0b0110_00_0) == M7_2(16.0)
@test M7_2(0b0110_00_1) == M7_2(24.0)
@test M7_2(0b0110_01_0) == M7_2(32.0)
@test M7_2(0b0110_01_1) == M7_2(48.0)
@test M7_2(0b0110_10_0) == M7_2(64.0)
@test M7_2(0b0110_10_1) == M7_2(96.0)
@test M7_2(0b0110_11_0) == M7_2(128.0)
@test M7_2(0b0110_11_1) == M7_2(192.0)
@test M7_2(0b01110_00)  == M7_2(256.0)
@test M7_2(0b01110_01)  == M7_2(512.0)
@test M7_2(0b01110_10)  == M7_2(1024.0)
@test M7_2(0b01110_11)  == M7_2(2048.0)
@test M7_2(0b011110_0)  == M7_2(4096.0)
@test M7_2(0b011110_1)  == M7_2(16384.0)
@test M7_2(0b0111110)   == M7_2(65536.0)
@test M7_2(0b0111111)   == M7_2(1048576.0)

#test cross-conversions.
for tile in M7_1
  @test tile == M7_1(Float64(tile))
end

for tile in M7_2
  @test tile == M7_2(Float64(tile))
end

#and a general purpose function for testing an operation against a matrix
function testop(op, iterable; exceptions = [])
  fails = 0
  totalsize = size(iterable) * size(iterable)
  for STile1 in iterable, STile2 in iterable
    fval1 = Float64(STile1)
    fval2 = Float64(STile2)

    texp = op(fval1, fval2)

    if isnan(texp) || (length(exceptions) > 0) && ([fval1, fval2] == exceptions) #catch NaN values.
      ok::Bool = false
      try
        op(STile1, STile2)
      catch e
        ok = isa(e, SigmoidNumbers.NaNError)
      end
      if !ok
        println("$iterable: $(fval1) $op $(fval2) failed to trap as NaN")
        fails += 1
      end
      continue
    end

    exp = iterable(texp)
    fexp = Float64(exp)

    res = op(STile1, STile2)
    fres = Float64(res)

    if (exp != res)
      println("$iterable: $(fval1) $op $(fval2) failed as $(fres); should be $(fexp)")
      fails += 1
    end

  end
  println("$iterable: $op $fails / $(totalsize) = $(100 * fails/totalsize)% failure!")
end

testop(+, M7, exceptions = [Inf, Inf])
testop(-, M7, exceptions = [Inf, Inf])
testop(*, M7)
testop(/, M7)

testop(+, M7_1, exceptions = [Inf, Inf])
testop(-, M7_1, exceptions = [Inf, Inf])
testop(*, M7_1)
testop(/, M7_1)

testop(+, M7_2, exceptions = [Inf, Inf])
testop(-, M7_2, exceptions = [Inf, Inf])
testop(*, M7_2)
testop(/, M7_2)

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
# Test variational information

using Base.Test
using Clustering

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
a3 = [3, 3, 3, 2, 2, 2, 1, 1, 1, 1]

(ARI, RI, MI, HI) = randindex(a1, a1)

@test isapprox(ARI, 1.0, atol=1.0e-12)
@test isapprox(RI , 1.0, atol=1.0e-12)
@test isapprox(MI , 0.0, atol=1.0e-12)
@test isapprox(HI , 1.0, atol=1.0e-12)

(ARI, RI, MI, HI) = randindex(a1, a3)

@test isapprox(ARI, 1.0, atol=1.0e-12)
@test isapprox(RI , 1.0, atol=1.0e-12)
@test isapprox(MI , 0.0, atol=1.0e-12)
@test isapprox(HI , 1.0, atol=1.0e-12)

(ARI, RI, MI, HI) = randindex(a1, a2)

@test isapprox(ARI, 0.403669, atol=1.0e-5)
@test isapprox(RI , 0.711111, atol=1.0e-5)
@test isapprox(MI , 0.288888, atol=1.0e-5)
@test isapprox(HI , 0.422222, atol=1.0e-5)

@test randindex(a1, a2) == randindex(a2, a1)


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function Control(t,x,rate_vector,data_dictionary)
# ---------------------------------------------------------------------- #
# Control.jl was generated using the Kwatee code generation system.
# Username: jeffreyvarner
# Type: GRN-JULIA
# Version: 1.0
# Generation timestamp: 08-08-2016 20:01:41
#
# Arguments:
# t  - current time
# x  - state vector
# rate_vector - vector of reaction rates
# data_dictionary  - Data dictionary instance (holds model parameters)
# ---------------------------------------------------------------------- #

# Set a default value for the allosteric control variables -
EPSILON = 1.0e-3;
#EPSILON = 1.0e-1;
number_of_reactions = length(rate_vector);
control_vector = ones(number_of_reactions);
control_parameter_array = data_dictionary["CONTROL_PARAMETER_ARRAY"];

# Alias the species vector -
gene_Smad = x[1];
gene_Ecadherin = x[2];
gene_Bcatenin = x[3];
gene_PLCy = x[4];
gene_PKC = x[5];
gene_RAF = x[6];
gene_RAS = x[7];
gene_AP1_SP1 = x[8];
gene_GSK3 = x[9];
gene_tubulin = x[10];
gene_SNAIL_SLUG = x[11];
gene_Calcineurin = x[12];
gene_APC_AXIN = x[13];
gene_PI3K = x[14];
mRNA_Smad = x[15];
mRNA_Ecadherin = x[16];
mRNA_Bcatenin = x[17];
mRNA_PLCy = x[18];
mRNA_PKC = x[19];
mRNA_RAF = x[20];
mRNA_RAS = x[21];
mRNA_AP1_SP1 = x[22];
mRNA_Calcineurin = x[23];
mRNA_SNAIL_SLUG = x[24];
mRNA_tubulin = x[25];
mRNA_APC_AXIN = x[26];
mRNA_GSK3 = x[27];
mRNA_PI3K = x[28];
PLCy = x[29];
Active_PLCy = x[30];
Calcineurin = x[31];
Active_Calcineurin = x[32];
NFATc = x[33];
Active_NFATc = x[34];
PKC = x[35];
Active_PKC = x[36];
RAF = x[37];
RAF_P = x[38];
Smad = x[39];
pSmad = x[40];
PI3K = x[41];
PI3K_P = x[42];
GSK3 = x[43];
GSK3_P = x[44];
RAS = x[45];
RAS_GTP = x[46];
MAPK = x[47];
pMAPK = x[48];
AP1_SP1 = x[49];
AP1_SP1_P = x[50];
Ecadherin = x[51];
Bcatenin = x[52];
Ecadherin_Bcatenin = x[53];
APC_AXIN = x[54];
APC_AXIN_Bcatenin = x[55];
TCF4 = x[56];
Active_TCF4 = x[57];
LEF1 = x[58];
Active_LEF1 = x[59];
PPase = x[60];
YREG1 = x[61];
SNAIL_SLUG = x[62];
vimentin = x[63];
tubulin = x[64];
RNAP = x[65];
RIBOSOME = x[66];
TGFB3_Ab = x[67];
ss_virus = x[68];
ecad_virus = x[69];
U0126 = x[70];
vivit = x[71];
gene_TCF4 = x[72];
gene_TGFB3 = x[73];
gene_YREG1 = x[74];
gene_LEF1 = x[75];
gene_vimentin = x[76];
gene_R = x[77];
gene_VEGFR = x[78];
gene_NFATc = x[79];
gene_MAPK = x[80];
mRNA_TCF4 = x[81];
mRNA_TGFB3 = x[82];
mRNA_YREG1 = x[83];
mRNA_LEF1 = x[84];
mRNA_vimentin = x[85];
mRNA_R = x[86];
mRNA_VEGFR = x[87];
mRNA_NFATc = x[88];
mRNA_MAPK = x[89];
R = x[90];
TGFB12 = x[91];
TGFB3 = x[92];
VEGF = x[93];
TGFB12_R = x[94];
TGFB3_R = x[95];
VEGFR = x[96];
VEGF_VEGFR = x[97];

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# TGFB3_inhibition target: TGFB3_binding actor: TGFB3_Ab type: inhibition
if (TGFB3_Ab<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[1,1]*(TGFB3_Ab)^control_parameter_array[1,2])/(1+control_parameter_array[1,1]*(TGFB3_Ab)^control_parameter_array[1,2]));
end

control_vector[3] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# NFATc_inhibition target: Calcineurin_activation_NFATc actor: vivit type: inhibition
if (vivit<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[2,1]*(vivit)^control_parameter_array[2,2])/(1+control_parameter_array[2,1]*(vivit)^control_parameter_array[2,2]));
end

control_vector[9] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# GSK3_P_inhibition_Binding_APC_AXIN_forw target: Binding_APC_AXIN_forw actor: GSK3_P type: inhibition
if (GSK3_P<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[3,1]*(GSK3_P)^control_parameter_array[3,2])/(1+control_parameter_array[3,1]*(GSK3_P)^control_parameter_array[3,2]));
end

control_vector[23] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# Active_TCF4_inhibition_Ecad target: TCF4_Bcatenin_complex_formation actor: Ecadherin_Bcatenin type: inhibition
if (Ecadherin_Bcatenin<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[4,1]*(Ecadherin_Bcatenin)^control_parameter_array[4,2])/(1+control_parameter_array[4,1]*(Ecadherin_Bcatenin)^control_parameter_array[4,2]));
end

# Active_TCF4_inhibition_APC target: TCF4_Bcatenin_complex_formation actor: APC_AXIN_Bcatenin type: inhibition
if (APC_AXIN_Bcatenin<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[5,1]*(APC_AXIN_Bcatenin)^control_parameter_array[5,2])/(1+control_parameter_array[5,1]*(APC_AXIN_Bcatenin)^control_parameter_array[5,2]));
end

control_vector[25] = mean(transfer_function_vector); # was maximum
#@show transfer_function_vector
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];
transfer_function_vector_a = Float64[];

# activation_Ecadherin_by_Active_NFATc target: induction_gene_Ecadherin actor: Active_NFATc type: induction
push!(transfer_function_vector_a,(control_parameter_array[6,1]*(Active_NFATc)^control_parameter_array[6,2])/(1+control_parameter_array[6,1]*(Active_NFATc)^control_parameter_array[6,2]));
# repression_Ecadherin_by_pSmad target: induction_gene_Ecadherin actor: pSmad type: repression
# if (pSmad<EPSILON);
# 	push!(transfer_function_vector,1.0);
# else
# 	push!(transfer_function_vector,1.0 - (control_parameter_array[7,1]*(pSmad)^control_parameter_array[7,2])/(1+control_parameter_array[7,1]*(pSmad)^control_parameter_array[7,2]));
# end

# repression_Ecadherin_by_SNAIL_SLUG target: induction_gene_Ecadherin actor: SNAIL_SLUG type: repression
if (SNAIL_SLUG<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[8,1]*(Active_LEF1*SNAIL_SLUG)^control_parameter_array[8,2])/(1+control_parameter_array[8,1]*(Active_LEF1*SNAIL_SLUG)^control_parameter_array[8,2]));
end

# repression_Ecadherin_by_Active_LEF1 target: induction_gene_Ecadherin actor: Active_LEF1 type: repression
if (pSmad<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[9,1]*(Active_LEF1*pSmad)^control_parameter_array[9,2])/(1+control_parameter_array[9,1]*(Active_LEF1*pSmad)^control_parameter_array[9,2]));
end

# activation_Ecadherin_by_ecad_virus target: induction_gene_Ecadherin actor: ecad_virus type: induction
push!(transfer_function_vector_a,(control_parameter_array[10,1]*(ecad_virus)^control_parameter_array[10,2])/(1+control_parameter_array[10,1]*(ecad_virus)^control_parameter_array[10,2]));
#control_vector[125] = maximum([minimum(transfer_function_vector),maximum(transfer_function_vector_a)]);
#control_vector[125] = mean([minimum(transfer_function_vector),maximum(transfer_function_vector_a)]); # more stable?
# function product(A)
# prod = 1
# for x in A
# prod = prod*x
# end
# return prod;
# end
#control_vector[125] = mean([product(transfer_function_vector),maximum(transfer_function_vector_a)]); # more stable?
#control_vector[125] = mean([minimum(transfer_function_vector[1:2]),maximum(transfer_function_vector_a),transfer_function_vector[3]]); # more stable?
#push!(transfer_function_vector,0.5)

#control_vector[125] = maximum([minimum(transfer_function_vector),maximum(transfer_function_vector_a),transfer_function_vector[3]]);
#control_vector[125] = maximum([mean([minimum(transfer_function_vector),transfer_function_vector[3]]),maximum(transfer_function_vector_a)]);

#control_vector[125] = mean([minimum(transfer_function_vector),maximum(transfer_function_vector_a),transfer_function_vector[3]]);
control_vector[125] = mean([minimum(transfer_function_vector),maximum(transfer_function_vector_a)]);

#control_vector[125] = mean([maximum([minimum(transfer_function_vector),transfer_function_vector[3]]),maximum(transfer_function_vector_a)]);

#control_vector[125] = maximum([minimum(transfer_function_vector),maximum(transfer_function_vector_a)]);

#control_vector[125] = mean([minimum(transfer_function_vector),maximum(transfer_function_vector_a),transfer_function_vector[3]]);

transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_TGFB3_by_TCF4 target: induction_gene_TGFB3 actor: Active_TCF4 type: induction
push!(transfer_function_vector,(control_parameter_array[11,1]*(Active_TCF4)^control_parameter_array[11,2])/(1+control_parameter_array[11,1]*(Active_TCF4)^control_parameter_array[11,2]));
# repression_TGFB3_by_YREG1 target: induction_gene_TGFB3 actor: YREG1 type: repression
if (YREG1<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[12,1]*(YREG1)^control_parameter_array[12,2])/(1+control_parameter_array[12,1]*(YREG1)^control_parameter_array[12,2]));
end

control_vector[128] = minimum(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# repression_YREG1_by_SNAIL_SLUG target: induction_gene_YREG1 actor: SNAIL_SLUG type: repression
if (SNAIL_SLUG<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[13,1]*(SNAIL_SLUG)^control_parameter_array[13,2])/(1+control_parameter_array[13,1]*(SNAIL_SLUG)^control_parameter_array[13,2]));
end

control_vector[129] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_LEF1_by_pSmad target: induction_gene_LEF1 actor: pSmad type: induction
push!(transfer_function_vector,(control_parameter_array[14,1]*(pSmad)^control_parameter_array[14,2])/(1+control_parameter_array[14,1]*(pSmad)^control_parameter_array[14,2]));
# repression_LEF1_by_YREG1 target: induction_gene_LEF1 actor: YREG1 type: repression
if (YREG1<EPSILON);
	push!(transfer_function_vector,1.0);
else
	push!(transfer_function_vector,1.0 - (control_parameter_array[15,1]*(YREG1)^control_parameter_array[15,2])/(1+control_parameter_array[15,1]*(YREG1)^control_parameter_array[15,2]));
end

#control_vector[130] = minimum([mean(transfer_function_vector),transfer_function_vector[2]]); #[2]
#control_vector[130] = minimum([mean(transfer_function_vector),transfer_function_vector[2]*transfer_function_vector[1]]); #[2]
#control_vector[130] = mean([transfer_function_vector[2]*transfer_function_vector[1],transfer_function_vector[1]])
control_vector[130] = mean(transfer_function_vector)
#control_vector[130]=transfer_function_vector[2]*transfer_function_vector[1]

#control_vector[130] = maximum([transfer_function_vector[2],transfer_function_vector[1]])

#control_vector[130] = minimum([transfer_function_vector[2]*transfer_function_vector[1],transfer_function_vector[1]]); #[2]

#control_vector[130] = mean(transfer_function_vector);
#@show mean(transfer_function_vector),transfer_function_vector[2]
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_vimentin_by_Active_LEF1 target: induction_gene_vimentin actor: Active_LEF1 type: induction
push!(transfer_function_vector,(control_parameter_array[16,1]*(Active_LEF1*pSmad)^control_parameter_array[16,2])/(1+control_parameter_array[16,1]*(Active_LEF1*pSmad)^control_parameter_array[16,2]));
# activation_vimentin_by_AP1_SP1_P target: induction_gene_vimentin actor: AP1_SP1_P type: induction
push!(transfer_function_vector,(control_parameter_array[17,1]*(AP1_SP1_P)^control_parameter_array[17,2])/(1+control_parameter_array[17,1]*(AP1_SP1_P)^control_parameter_array[17,2]));
control_vector[132] = mean(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# ----------------------------------------------------------------------------------- #
transfer_function_vector = Float64[];

# activation_SNAIL_SLUG_by_AP1_SP1_P target: induction_gene_SNAIL_SLUG actor: AP1_SP1_P type: induction
push!(transfer_function_vector,(control_parameter_array[18,1]*(AP1_SP1_P)^control_parameter_array[18,2])/(1+control_parameter_array[18,1]*(AP1_SP1_P)^control_parameter_array[18,2]));
# activation_SNAIL_SLUG_by_ss_virus target: induction_gene_SNAIL_SLUG actor: ss_virus type: induction
push!(transfer_function_vector,(control_parameter_array[19,1]*(ss_virus)^control_parameter_array[19,2])/(1+control_parameter_array[19,1]*(ss_virus)^control_parameter_array[19,2]));
control_vector[143] = maximum(transfer_function_vector);
transfer_function_vector = 0;
# ----------------------------------------------------------------------------------- #

# Modify the rate_vector with the control variables -
rate_vector = rate_vector.*control_vector;

# Return the modified rate vector -
return rate_vector;
end


text:
using DifferentialEquations
using StochasticDelayDiffEq
using DiffEqCallbacks
using Random
using SparseArrays
using Tables
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_dde!(du,u,h,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    infection = β*c*I/N*S
    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
    Nd = Sd+Id+Rd
    recovery = β*c*Id/Nd*Sd
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
    end
    nothing
end;

# Define a sparse matrix by making a dense matrix and setting some values as not zero
A = zeros(3,2)
A[1,1] = 1
A[2,1] = 1
A[2,2] = 1
A[3,2] = 1
A = SparseArrays.sparse(A);

# Make `g` write the sparse matrix values
function sir_delayed_noise!(du,u,h,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    infection = β*c*I/N*S
    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
    Nd = Sd+Id+Rd
    recovery = β*c*Id/Nd*Sd
    du[1,1] = -sqrt(infection)
    du[2,1] = sqrt(infection)
    du[2,2] = -sqrt(recovery)
    du[3,2] = sqrt(recovery)
end;

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[2]
end;

function affect!(integrator)
  integrator.u[2] = 0.0
end;

cb = ContinuousCallback(condition,affect!);

function affect_initial_recovery!(integrator)
    integrator.u[2] -= u0[2]
    integrator.u[3] += u0[2]

    reset_aggregated_jumps!(integrator)
end
cb_initial_recovery = DiscreteCallback((u,t,integrator) -> t == p[3], affect_initial_recovery!)

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0]; # S,I,R

function sir_history(p, t)
    [1000.0, 0.0, 0.0]
end;

p = [0.05,10.0,4.0]; # β,c,τ

Random.seed!(1234);

prob_sdde = SDDEProblem(sir_dde!,sir_delayed_noise!,u0,sir_history,tspan,p;noise_rate_prototype=A);

sol_sdde = solve(prob_sdde,LambaEM(),callback=CallbackSet(cb,cb_initial_recovery));

df_sdde = DataFrame(Tables.table(sol_sdde(t)'))
rename!(df_sdde,["S","I","R"])
df_sdde[!,:t] = t;

@df df_sdde plot(:t,
    [:S,:I,:R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

#@benchmark solve(prob_sdde,LambaEM(),callback=CallbackSet(cb,cb_initial_recovery))


text:
import DataStructures: Deque

mutable struct OrderBook
    price_levels::PriceLevels
    order_id::UInt128 # Does it need to be this big?
    orders::Dict{UInt128, LimitOrder}

    function OrderBook(max_price::UInt128, min_price::UInt128)::OrderBook
        new(PriceLevels(max_price, min_price), UInt128(0), Dict{UInt128, LimitOrder}())
    end
end

function limitorder!(ob::OrderBook, side::Bool, price::UInt128,
                     size::UInt128, trader_id::UInt128)::Tuple{Deque{LimitOrder}, UInt128}
    # Init trades
    trades = Deque{LimitOrder}()

    # Matching
    if side == BUY
        ask = best_ask(ob)
        while (price >= ask) & (size > 0)
            price_level = get_level(ob.price_levels, ask)
            trades, size = match!(price_level, size, trades)
            if price_level.size == 0
                ask = update_ask!(ob.price_levels)
            end
        end
    # Matches a sell order to buy orders
    else
        bid = best_bid(ob)
        while (price <= bid) & (size > 0)
            price_level = get_level(ob.price_levels, bid)
            trades, size = match!(price_level, size, trades)
            if price_level.size == 0
                bid = update_bid!(ob.price_levels)
            end
        end
    end

    if size > 0
        # After Matching, just input the limit order
        ob.order_id = ob.order_id + 1
        order = LimitOrder(price, size, trader_id, side, ob.order_id)
        ob.orders[ob.order_id] = order

        # Input order
        insert_order!(ob.price_levels, order)

        return trades, ob.order_id
    else
        return trades, UInt128(0)
    end
end

function marketorder!(ob::OrderBook, side::Bool, size::UInt128,
                     trader_id::UInt128)

    return 0
end

function cancel!(ob::OrderBook, order_id::UInt128)

    @assert order_id in keys(ob.orders) "order_id not in orderbook"

    order = pop!(ob.orders, order_id)
    price_level = get_level(ob.price_levels, order.price)
    remove!(price_level, order_id)

    if (order.price == best_bid(ob)) & (size(price_level) == 0)
        update_bid!(ob.price_levels)
    elseif (order.price == best_ask(ob)) & (size(price_level) == 0)
        update_ask!(ob.price_levels)
    end
end

function update!(ob::OrderBook, order_id::UInt128, size::UInt128)
    order = ob.orders[order_id]
    diff = size - order.size
    order.size = size
    pl = get_level(ob.price_levels, order.price)
    pl.size = pl.size + diff
end

function best_bid(ob::OrderBook)::UInt128
    return ob.price_levels.bid
end

function best_bid_size(ob::OrderBook)::UInt128
    bid = best_bid(ob)
    return size(get_level(ob.price_levels, bid))
end

function best_ask(ob::OrderBook)::UInt128
    return ob.price_levels.ask
end

function best_ask_size(ob::OrderBook)::UInt128
    ask = best_ask(ob)
    return size(get_level(ob.price_levels, ask))
end


text:
"""
Copyright 2020 Oghogho Odemwingie <odemwingieog@gmail.com>. All rights reserved.
Use of this source code is governed by an MIT
license that can be found in the LICENSE file.

Package blowfish implements Bruce Schneier's Blowfish encryption algorithm.

Blowfish is a legacy cipher and its short block size makes it vulnerable to
birthday bound attacks (see https://sweet32.info). It should only be used
where compatibility with legacy systems, not security, is the goal.
"""

module Blowfish

export Blocksize, Cipher, KeySizeError, NewCipher, NewSaltedCipher, Encrypt!, Decrypt!, ExpandKey!

# The code is a port of Bruce Schneier's C implementation.
# See https://www.schneier.com/blowfish.html.

const Blocksize = 8

"""
	A KeySizeError is an Exception which indicates that the length of a given cipher key does not conform to the bounds key -> 1 < key < 56
"""
struct KeySizeError <: Exception
	size::Int
end
Base.show(kErr::KeySizeError) = "Blowfish.jl: invalid key size ($(kErr.size))"

"""
	A Cipher is an instance of Blowfish encryption using a particular key.
"""
mutable struct Cipher
	p::Array{UInt32, 1}
	s0::Array{UInt32, 1}
	s1::Array{UInt32, 1}
	s2::Array{UInt32, 1}
	s3::Array{UInt32, 1}

	Cipher() = new(collect(p), collect(s0), collect(s1), collect(s2), collect(s3))
end

include("Block.jl")
include("Const.jl")

"""
	NewCipher(key::Array{UInt8, 1})

Generates a new instance of a blowfish cipher from a given cipher key

# Examples
```julia-repl
julia> NewCipher(UInt8['a', 'e', 'i', 'o', 'u'])
Cipher(UInt32[0x7a48b0ec, 0xb7c5649b, 0xfee9c8a4, 0x985e8a2a, 0xd6a5938e, 0x145c1599, 0x51b0c8db, 0x91b92643, 0x53c0affc, 0x3642f567, 0x8100d42b, 0x49179092, 0x7d3bec7b, 0x3ba790e3, 0xc3e03bcd, 0x85640a20, 0xdeb3c3dc, 0xa49501f5], UInt32[0xde18af37, 0xd5438974, 0x461bbb53, 0x392f99c3, 0xfb69ec53, 0x049fe21a, 0xd4489b31, 0xc73d87f2, 0xfdb21b2e, 0xb78761e0  …  0x458175a4, 0xc13c6dd7, 0x1fa5ee57, 0x38925136, 0x29f87d84, 0x36c2ce5c, 0x76ec4eb0, 0x7b4a5c7c, 0x8fb372ec, 0xc775630f], UInt32[0x703686bc, 0xd1acb9d6, 0x2832c1f9, 0xecb4b1f3, 0xecb60b8d, 0xc22e7889, 0x79fb10f0, 0xb7230280, 0xbd6772e9, 0xb320259b  …  0x18610527, 0xb6e95540, 0xdd5a53b3, 0x6c144d00, 0x8254f230, 0xcc9cd03b, 0xda916333, 0x64818319, 0x4d123b0e, 0x17c4af84], UInt32[0xe1919f40, 0x50b5ebc1, 0x8673f868, 0x2e7e8c6d, 0xc4e0bf50, 0x4769fa9f, 0x7e55db03, 0x21e58059, 0xc4d0c578, 0x8239570e  …  0xaffbcc32, 0x480e8087, 0x0a41bc51, 0xb4571997, 0xc01ecce2, 0x32d89616, 0xc671c03e, 0xaf2f65a2, 0x2f073414, 0xa780d378], UInt32[0x5f19d24e, 0x335b6476, 0x9fb61c67, 0x17e3c8f5, 0x98aab626, 0x7c4c102f, 0x03bbf3a2, 0xade7d94a, 0x6f0bc1d2, 0x9a6da0e8  …  0xc85d33cf, 0x73a34514, 0x30600f9e, 0x5afc4493, 0x64fbb2ff, 0x1d0dacce, 0x9468581e, 0xd3ea966a, 0x10a864b3, 0xc4d48364])
```
"""
function NewCipher(key::AbstractArray{UInt8, 1}) :: Cipher
	(length(key) < 1 || length(key) > 56) && throw(KeySizeError(length(key)))
	res = Cipher()
	ExpandKey!(key, res)
	return res
end

"""
	NewSaltedCipher(key::Array{UInt8, 1})

Generates a Cipher that folds a salt into its key schedule. For most purposes, NewCipher, instead of NewSaltedCipher, is sufficient and desirable. For bcrypt compatibility, the key can be over 56 bytes.

# Examples
```julia-repl
julia> NewSaltedCipher(UInt8['s', 'a', 'l', 't'],UInt8['a', 'e', 'i', 'o', 'u'])
Cipher(UInt32[0xab5a9c73, 0xa49c2d2a, 0xc87e8ec8, 0xaa432b79, 0x1d0b7db7, 0x21cbb7fe, 0x2f36fc0f, 0x2871116d, 0xe7f0000b, 0x10264298, 0x545b2686, 0x3eb8dcde, 0xf6e314ab, 0x7d4f6dd5, 0x2e89833f, 0xe16df037, 0x61e28e24, 0xc5409a31], UInt32[0xc53c53e7, 0xfb819947, 0x23c33744, 0x795a4bf1, 0xdd4183b4, 0xae68d5c4, 0x9eb3ed52, 0x7b3e0314, 0x16e20c08, 0x779e49a0  …  0x1bceff46, 0x50a2f67b, 0x30c5d692, 0xb6aa33c1, 0x201dd700, 0x640abd32, 0x123498f4, 0x84383581, 0x6e6b9f74, 0x8ae4ae80], UInt32[0xe6f14e74, 0x2072fdce, 0x3318cee1, 0xb31f5a86, 0x33ad0de5, 0xab1423b6, 0x661392ad, 0x437999f9, 0x6b2e4675, 0xd72bb28c  …  0xba487480, 0x2f82dfcb, 0xff1c012d, 0x56705613, 0xe21b7f92, 0xf4a3ee42, 0x1f9c5542, 0x01e8f0ac, 0x000e41eb, 0x8ec4511b], UInt32[0xebbd88cb, 0xde0bc3b4, 0x23c31968, 0xa25fd135, 0xec8689f8, 0x70fcb881, 0xc74964a3, 0xe59851c2, 0x3cdb369f, 0xc2c4f08a  …  0x464ba53e, 0xef69b586, 0x196b4813, 0x1ec8abc6, 0x8ab71719, 0x34a15e09, 0x7675589d, 0x9d7f2a78, 0x92540a14, 0x1452cb64], UInt32[0x62053d78, 0xcc13a352, 0xb4b0d494, 0x92d66e7f, 0x2db90c93, 0xf5243558, 0xe0831a00, 0xca5f77b7, 0x99f6388a, 0x45906ae3  …  0xdd979889, 0x961e8abd, 0x5c821e51, 0x8a2fbfcc, 0x292c2a3c, 0xdbc24193, 0xdfdcd502, 0x22868810, 0x2d0a2a98, 0xa5fc3e68])
```
"""
function NewSaltedCipher(key::AbstractArray{UInt8, 1}, salt::AbstractArray{UInt8, 1}) :: Cipher
	isempty(salt) && return NewCipher(key)
	length(key) < 1 && throw(KeySizeError(length(key)))
	res = Cipher()
	expand_key_with_salt!(key, salt, res)
	return res
end

"""
	Encrypt!(c::Cipher, dst::Array{UInt8, 1}, src::Array{UInt8, 1})

Encrypt encrypts the 8-byte buffer src using the key k and stores the result in dst.
Note that for amounts of data larger than a block, it is not safe to just call Encrypt on successive blocks.
"""
function Encrypt!(c::Cipher, dst::AbstractArray{UInt8, 1}, src::AbstractArray{UInt8, 1})
	l = UInt32(src[1])<<24 | UInt32(src[2])<<16 | UInt32(src[3])<<8 | UInt32(src[4])
	r = UInt32(src[5])<<24 | UInt32(src[6])<<16 | UInt32(src[7])<<8 | UInt32(src[8])
	l, r = encrypt_block(l, r, c)
	dst[1], dst[2], dst[3], dst[4] = UInt8(l>>24 & 0xff), UInt8(l>>16 & 0xff), UInt8(l>>8 & 0xff), UInt8(l & 0xff)
	dst[5], dst[6], dst[7], dst[8] = UInt8(r>>24 & 0xff), UInt8(r>>16 & 0xff), UInt8(r>>8 & 0xff), UInt8(r & 0xff)
end

"""
	Decrypt!(c::Cipher, dst::Array{UInt8, 1}, src::Array{UInt8, 1})

Decrypt decrypts the 8-byte buffer src using the key k and stores the result in dst.
"""
function Decrypt!(c::Cipher, dst::AbstractArray{UInt8, 1}, src::AbstractArray{UInt8, 1})
	l = UInt32(src[1])<<24 | UInt32(src[2])<<16 | UInt32(src[3])<<8 | UInt32(src[4])
	r = UInt32(src[5])<<24 | UInt32(src[6])<<16 | UInt32(src[7])<<8 | UInt32(src[8])
	l, r = decrypt_block(l, r, c)
	dst[1], dst[2], dst[3], dst[4] = UInt8(l>>24 & 0xff), UInt8(l>>16 & 0xff), UInt8(l>>8 & 0xff), UInt8(l & 0xff)
	dst[5], dst[6], dst[7], dst[8] = UInt8(r>>24 & 0xff), UInt8(r>>16 & 0xff), UInt8(r>>8 & 0xff), UInt8(r & 0xff)
end

end # module


text:
function binpacking_vectorofvariables_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    w1 = 2
    w2 = 2

    c1 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x2, x3]),
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [w1, w2]),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 4
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 0
    end
end

function binpacking_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING, Int},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    w1 = 2
    w2 = 2

    c1 = MOI.add_constraint(
        model,
        MOIU.vectorize([x1, x2, x3]),
        CP.BinPacking{CP.NO_CAPACITY_BINPACKING}(1, 2, [w1, w2]),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, c1)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 4
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 0
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 0
    end
end

# @testset "BinPacking: ScalarAffineFunction with variable sizes" begin
#     model = OPTIMIZER
#     MOI.empty!(model)

#     @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
#     @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Int}, MOI.EqualTo{Int})
#     @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Int}, CP.BinPacking)

#     x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
#     x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
#     x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
#     x4, _ = MOI.add_constrained_variable(model, MOI.Integer())
#     x5, _ = MOI.add_constrained_variable(model, MOI.Integer())

#     c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1], [x4]), 0), MOI.EqualTo(2))
#     c1 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1], [x5]), 0), MOI.EqualTo(2))

# TODO: currently, throws an assertion. https://github.com/dourouc05/ConstraintProgrammingExtensions.jl/issues/4
#     @test_throws(
#         MOI.AddConstraintNotAllowed{MOI.VectorAffineFunction{Int64}, CP.BinPacking},
#         MOI.add_constraint(model,
#                            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 2, 3, 4, 5], MOI.ScalarAffineTerm.([1, 1, 1, 1, 1], [x1, x2, x3, x4, x5])), [0, 0, 0, 0, 0]),
#                            CP.BinPacking(1, 2))
#     )
# end

const binpackingtests = Dict(
    "binpacking_vectorofvariables" => binpacking_vectorofvariables_test,
    "binpacking_scalaraffinefunction" =>
        binpacking_scalaraffinefunction_test,
)

MOIT.@moitestset binpacking


text:
@testset "$(rpad("QuadratureRule",80))" begin

    o = 2
    b = [0.5, 0.5]
    c = [0.0, 1.0]

    quad = QuadratureRule(o,c,b)

    @test order(quad) == o
    @test nodes(quad) == c
    @test weights(quad) == b
    @test nnodes(quad) == length(b)
    @test eachindex(quad) == 1:2
    @test eltype(quad) == Float64

    @test quad(x -> 1)   == 1
    @test quad(x -> x)   == 1//2
    @test quad(x -> x^2) == 1//2

    quad2 = QuadratureRule{Float64,2}(o, c, b)
    quadϵ = QuadratureRule{Float64,2}(o, c .+ 8eps() * rand(2), b .+ 8eps() * rand(2))

    @test hash(quad) == hash(quad2)
    @test hash(quad) != hash(quadϵ)

    @test  isequal(quad, quad2)
    @test !isequal(quad, quadϵ)

    @test isapprox(quad, quad2)
    @test isapprox(quad, quadϵ)

    @test quad ≈   quad2
    @test quad ==  quad2
    @test quad === quad2

    @test quad ≈   quadϵ
    @test quad !=  quadϵ
    @test quad !== quadϵ

end


text:
# Some global variables that will stay constant
α = 0.013
α_q = (1-(1-α)^3)
b_param = 0.0124
d_param = 0.00822
β = 0.98
γ = 1.0
σ = 2.0

# The default wage distribution: a discretized log normal
log_wage_mean, wage_grid_size, max_wage = 20, 200, 170
w_vec = linspace(1e-3, max_wage, wage_grid_size + 1)
logw_dist = Normal(log(log_wage_mean), 1)
cdf_logw = cdf.(logw_dist, log.(w_vec))
pdf_logw = cdf_logw[2:end] - cdf_logw[1:end-1]
p_vec = pdf_logw ./ sum(pdf_logw)
w_vec = (w_vec[1:end-1] + w_vec[2:end]) / 2

"""
Compute the reservation wage, job finding rate and value functions of the
workers given c and τ.

"""
function compute_optimal_quantities(c::AbstractFloat, τ::AbstractFloat)
    mcm = McCallModel(α_q,
                      β,
                      γ,
                      c-τ,                # post-tax compensation
                      σ,
                      collect(w_vec-τ),  # post-tax wages
                      p_vec)

    w_bar, V, U = compute_reservation_wage(mcm, return_values=true)
    λ = γ * sum(p_vec[w_vec - τ .> w_bar])

    return w_bar, λ, V, U
end

"""
Compute the steady state unemployment rate given c and tau using optimal
quantities from the McCall model and computing corresponding steady state
quantities

"""
function compute_steady_state_quantities(c::AbstractFloat, τ::AbstractFloat)
    w_bar, λ_param, V, U = compute_optimal_quantities(c, τ)

    # Compute steady state employment and unemployment rates
    lm = LakeModel(λ=λ_param, α=α_q, b=b_param, d=d_param)
    x = rate_steady_state(lm)
    u_rate, e_rate = x

    # Compute steady state welfare
    w = sum(V .* p_vec .* (w_vec - τ .> w_bar)) / sum(p_vec .* (w_vec - τ .> w_bar))
    welfare = e_rate .* w + u_rate .* U

    return u_rate, e_rate, welfare
end

"""
Find tax level that will induce a balanced budget.

"""
function find_balanced_budget_tax(c::Real)
    function steady_state_budget(t::Real)
      u_rate, e_rate, w = compute_steady_state_quantities(c, t)
      return t - u_rate * c
    end

    τ = brent(steady_state_budget, 0.0, 0.9 * c)

    return τ
end

# Levels of unemployment insurance we wish to study
Nc = 60
c_vec = linspace(5.0, 140.0, Nc)

tax_vec = Vector{Float64}(Nc)
unempl_vec = Vector{Float64}(Nc)
empl_vec = Vector{Float64}(Nc)
welfare_vec = Vector{Float64}(Nc)

for i = 1:Nc
    t = find_balanced_budget_tax(c_vec[i])
    u_rate, e_rate, welfare = compute_steady_state_quantities(c_vec[i], t)
    tax_vec[i] = t
    unempl_vec[i] = u_rate
    empl_vec[i] = e_rate
    welfare_vec[i] = welfare
end

fig, axes = subplots(2, 2, figsize=(15, 10))

plots = [unempl_vec, empl_vec, tax_vec, welfare_vec]
titles = ["Unemployment", "Employment", "Tax", "Welfare"]

for (ax, plot, title) in zip(axes, plots, titles)
    ax[:plot](c_vec, plot, "b-", lw=2, alpha=0.7)
    ax[:set](title=title)
end

fig[:tight_layout]()


text:
using CSyntax
using GLTF

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))
include(joinpath(@__DIR__, "camera.jl"))

# init window
width, height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# we use glTF file format because it's OpenGL friendly and we don't need to care about
# how to parse and extract mesh data from other mesh file format like .obj format.
# load glTF file
suzanne = GLTF.load(joinpath(@__DIR__, "monkey2.gltf"))
suzanne_data = [read(joinpath(@__DIR__, b.uri)) for b in suzanne.buffers]

# load suzanne position metadata
search_name(x, keyword) = x[findfirst(x->occursin(keyword, x.name), x)]
pos = search_name(suzanne.accessors, "positions")
pos_bv = suzanne.bufferViews[pos.bufferView]
# load suzanne index metadata
indices = search_name(suzanne.accessors, "indices")
idx_bv = suzanne.bufferViews[indices.bufferView]
# load suzanne texture coordinate metadata
texcoords = search_name(suzanne.accessors, "texcoords")
tex_bv = suzanne.bufferViews[texcoords.bufferView]
# load suzanne normal metadata
normals = search_name(suzanne.accessors, "normals")
normal_bv = suzanne.bufferViews[normals.bufferView]

# create VBO and EBO
pos_vbo = GLuint(0)
@c glGenBuffers(1, &pos_vbo)
glBindBuffer(pos_bv.target, pos_vbo)
glBufferData(pos_bv.target, pos_bv.byteLength, C_NULL, GL_STATIC_DRAW)
pos_data = suzanne_data[pos_bv.buffer]
@c glBufferSubData(pos_bv.target, 0, pos_bv.byteLength, &pos_data[pos_bv.byteOffset])

normal_vbo = GLuint(0)
@c glGenBuffers(1, &normal_vbo)
glBindBuffer(normal_bv.target, normal_vbo)
glBufferData(normal_bv.target, normal_bv.byteLength, C_NULL, GL_STATIC_DRAW)
normal_data = suzanne_data[normal_bv.buffer]
@c glBufferSubData(normal_bv.target, 0, normal_bv.byteLength, &normal_data[normal_bv.byteOffset])

tex_vbo = GLuint(0)
@c glGenBuffers(1, &tex_vbo)
glBindBuffer(tex_bv.target, tex_vbo)
glBufferData(tex_bv.target, tex_bv.byteLength, C_NULL, GL_STATIC_DRAW)
tex_data = suzanne_data[tex_bv.buffer]
@c glBufferSubData(tex_bv.target, 0, tex_bv.byteLength, &tex_data[tex_bv.byteOffset])

idx_ebo = GLuint(0)
@c glGenBuffers(1, &idx_ebo)
glBindBuffer(idx_bv.target, idx_ebo)
glBufferData(idx_bv.target, idx_bv.byteLength, C_NULL, GL_STATIC_DRAW)
idx_data = suzanne_data[idx_bv.buffer]
@c glBufferSubData(idx_bv.target, 0, idx_bv.byteLength, &idx_data[idx_bv.byteOffset])

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
# bind position vbo
glBindBuffer(pos_bv.target, pos_vbo)
glVertexAttribPointer(0, 3, pos.componentType, pos.normalized, pos_bv.byteStride, Ptr{Cvoid}(pos.byteOffset))
glEnableVertexAttribArray(0)
# bind normal vbo
glBindBuffer(normal_bv.target, normal_vbo)
glVertexAttribPointer(1, 3, normals.componentType, normals.normalized, normal_bv.byteStride, Ptr{Cvoid}(normals.byteOffset))
glEnableVertexAttribArray(1)
# bind texture coordinate vbo
glBindBuffer(tex_bv.target, tex_vbo)
glVertexAttribPointer(2, 2, texcoords.componentType, texcoords.normalized, tex_bv.byteStride, Ptr{Cvoid}(texcoords.byteOffset))
glEnableVertexAttribArray(2)

# camera
camera = PerspectiveCamera()
setposition!(camera, [0.0, 0.0, 5.0])

# create shader program
vert_shader = createshader(joinpath(@__DIR__, "mesh.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "mesh.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)
view_loc = glGetUniformLocation(shader_prog, "view")
proj_loc = glGetUniformLocation(shader_prog, "proj")
glUseProgram(shader_prog)
glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))

# # enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CCW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

let
updatefps = FPSCounter()
# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glBindBuffer(idx_bv.target, idx_ebo)
    glDrawElements(GL_TRIANGLES, indices.count, indices.componentType, Ptr{Cvoid}(0))
    # check and call events
    GLFW.PollEvents()
    # move camera
    updatecamera!(window, camera)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
module Boards

using PyCall
@pyimport gpu_lock2 as gg

export
    getboard

function getboard()
    board_id = gg.obtain_lock_id()

    if board_id < 0
        error("could not obtain lock")
    end

    board_id
end

end # module


text:
using RecurrenceAnalysis, SparseArrays
using Test, Statistics, LinearAlgebra, DelimitedFiles

cells = [(1,1),(2,2),(3,2),(6,2),(7,2),(8,2),(6,3),(7,3),(10,3),(2,4),(4,4),
    (1,5),(2,5),(3,5),(5,5),(9,5),(10,5),(2,6),(3,6),(6,6),(7,6),(10,6),
    (3,7),(8,7),(4,8),(1,9),(5,9),(7,9),(9,9),(9,10),(10,10),(3,11),(7,11)]
i = [x[1] for x in cells]
j = [x[2] for x in cells]
rmat = CrossRecurrenceMatrix(sparse(i,j,trues(length(i))))
### Graph of the recurrence matrix:
## · : non-recurrent (white) points
## x : recurrent (black) points
## o : recurrent points within Theiler window = 2
#
#    1 2 3 4 5 6 7 8 9 A B
#  1 o · · · x · · · x · ·
#  2 · o · x x x · · · · ·
#  3 · o · · x x x · · · x
#  4 · · · o · · · x · · ·
#  5 · · · · o · · · x · ·
#  6 · x x · · o · · · · ·
#  7 · x x · · o · · x · x
#  8 · x · · · · o · · · ·
#  9 · · · · x · · · o o ·
#  A · · x · x x · · · o ·
#
##

@testset "Recurrence structures" begin
    @testset "Default parameters" begin
        histograms = recurrencestructures(rmat)
        true_histograms = Dict("diagonal" => [11,7,1,0,1],
                               "vertical" => [15,6,2],
                               "recurrencetimes" => [3,1,4,4])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=0, lmin=1, border=1)
        @test rqa_params[:RR] == 33/110
        @test rqa_params[:DET] == 1.0
        @test rqa_params[:L] == 33/20
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.996 atol=0.001
        @test rqa_params[:TREND] ≈ -33.8 atol=0.1
        @test rqa_params[:LAM] == 1.0
        @test rqa_params[:TT] == 33/23
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 33/12
        @test rqa_params[:RTE] ≈ 1.286 atol=0.001
        @test rqa_params[:NMPRT] == 4
    end
    @testset "With Theiler window" begin
        histograms = recurrencestructures(rmat, theiler=2)
        true_histograms = Dict("diagonal" => [9,4,0,0,1],
                               "vertical" => [10,3,2],
                               "recurrencetimes" => [1,0,3,0,0,0,2])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=2, lmin=1, border=1)
        @test rqa_params[:RR] == 22/110
        @test rqa_params[:DET] == 1.0
        @test rqa_params[:L] == 22/14
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.830 atol=0.001
        @test rqa_params[:TREND] ≈ -12.2 atol=0.1
        @test rqa_params[:LAM] == 1.0
        @test rqa_params[:TT] == 22/15
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 24/6
        @test rqa_params[:RTE] ≈ 1.011 atol=0.001
        @test rqa_params[:NMPRT] == 3
    end
    @testset "With minimum line" begin
        histograms = recurrencestructures(rmat, lmin=2)
        true_histograms = Dict("diagonal" => [11,7,1,0,1],
                               "vertical" => [15,6,2],
                               "recurrencetimes" => [3,1,4,4])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=0, lmin=2, border=1)
        @test rqa_params[:RR] == 33/110
        @test rqa_params[:DET] == 22/33
        @test rqa_params[:L] == 22/9
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.684 atol=0.001
        @test rqa_params[:TREND] ≈ -33.8 atol=0.1
        @test rqa_params[:LAM] == 18/33
        @test rqa_params[:TT] == 18/8
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 33/12
        @test rqa_params[:RTE] ≈ 1.286 atol=0.001
        @test rqa_params[:NMPRT] == 4
    end
    @testset "Theiler and minimum line" begin
        histograms = recurrencestructures(rmat, theiler=2, lmin=2)
        true_histograms = Dict("diagonal" => [9,4,0,0,1],
                               "vertical" => [10,3,2],
                               "recurrencetimes" => [1,0,3,0,0,0,2])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=2, lmin=2, border=1)
        @test rqa_params[:RR] == 22/110
        @test rqa_params[:DET] == 13/22
        @test rqa_params[:L] == 13/5
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.500 atol=0.001
        @test rqa_params[:TREND] ≈ -12.2 atol=0.1
        @test rqa_params[:LAM] == 12/22
        @test rqa_params[:TT] == 12/5
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 24/6
        @test rqa_params[:RTE] ≈ 1.011 atol=0.001
        @test rqa_params[:NMPRT] == 3
    end
end

## Recurrence network
@testset "Recurrence networks" begin
    # 7 edges, 12 linked triples, 1 triangle (1-2-4)
    adjmat7 =  [0 1 0 1 0 0
                1 0 1 1 1 0
                0 1 0 0 0 0
                1 1 0 0 0 1
                0 1 0 0 0 1
                0 0 0 1 1 0]
    rna_dict = rna(RecurrenceMatrix(adjmat7))
    @test rna_dict[:density] == 7/15
    @test rna_dict[:transitivity] == 0.25 # (3/12)
    # taken from Donner et al. https://doi.org/10.1088/1367-2630/12/3/033025
    adjmat10 = [0 0 0 1 0 0 0 1 0 0
                0 0 0 0 1 0 0 0 1 0
                0 0 0 0 0 1 0 0 0 1
                1 0 0 0 0 0 1 0 0 0
                0 1 0 0 0 0 0 1 0 0
                0 0 1 0 0 0 0 0 1 0
                0 0 0 1 0 0 0 0 0 1
                1 0 0 0 1 0 0 0 0 0
                0 1 0 0 0 1 0 0 0 0
                0 0 1 0 0 0 1 0 0 0]
    dismat10 = [0 3 4 1 2 5 2 1 4 3
                3 0 3 4 1 2 5 2 1 4
                4 3 0 3 4 1 2 5 2 1
                1 4 3 0 3 4 1 2 5 2
                2 1 4 3 0 3 4 1 2 5
                5 2 1 4 3 0 3 4 1 2
                2 5 2 1 4 3 0 3 4 1
                1 2 5 2 1 4 3 0 3 4
                4 1 2 5 2 1 4 3 0 3
                3 4 1 2 5 2 1 4 3 0]
    rna_dict = rna(RecurrenceMatrix(adjmat10))
    @test rna_dict[:density] == 2/9
    @test rna_dict[:transitivity] == 0
    @test rna_dict[:averagepath] ≈ sum(dismat10) / (10*9)
    @test rna_dict[:diameter] == maximum(dismat10)
    adjmat9 =  [0 0 0 0 1 0 0 1 1
                0 0 0 0 0 1 0 0 1
                0 0 0 1 0 0 1 0 0
                0 0 1 0 1 0 1 1 0
                1 0 0 1 0 0 1 1 1
                0 1 0 0 0 0 0 0 0
                0 0 1 1 1 0 0 0 0
                1 0 0 1 1 0 0 0 1
                1 1 0 0 1 0 0 1 0]
    dismat9 =  [0 1 2 2 1 2 3 1 1
                1 0 2 3 1 1 2 2 1
                2 2 0 4 1 1 1 3 2
                2 3 4 0 3 4 5 1 3
                1 1 1 3 0 1 2 2 1
                2 1 1 4 1 0 1 3 2
                3 2 1 5 2 1 0 4 3
                1 2 3 1 2 3 4 0 2
                1 1 2 3 1 2 3 2 0]
    rna_dict = rna(RecurrenceMatrix(adjmat9))
    triples = [3, 1, 1, 6, 10, 0, 3, 6, 6]
    triangles = [3, 0, 1, 3, 5, 0, 2, 4, 3]
    @test rna_dict[:density] == 7/18
    @test rna_dict[:transitivity] ≈ sum(triangles) / sum(triples)
    @test rna_dict[:averagepath] ≈ sum(dismat9) / (9*8)
    @test rna_dict[:diameter] == maximum(dismat9)
    adjmat11 = [0 0 0 1 0 0 0 1 0 0
                0 0 0 0 1 0 0 0 1 0
                0 0 0 0 0 1 0 0 0 0
                1 0 0 0 0 0 1 0 0 0
                0 1 0 0 0 0 0 1 0 0
                0 0 1 0 0 0 0 0 1 0
                0 0 0 1 0 0 0 0 0 0
                1 0 0 0 1 0 0 0 0 0
                0 1 0 0 0 1 0 0 0 0
                0 0 0 0 0 0 0 0 0 0]
    dismat11 =  [0 3 6 1 2 5 2 1 4 Inf
                 3 0 3 4 1 2 5 2 1 Inf
                 6 3 0 7 4 1 8 5 2 Inf
                 1 4 7 0 3 6 1 2 5 Inf
                 2 1 4 3 0 3 4 1 2 Inf
                 5 2 1 6 3 0 7 4 1 Inf
                 2 5 8 1 4 7 0 3 6 Inf
                 1 2 5 2 1 4 3 0 3 Inf
                 4 1 2 5 2 1 6 3 0 Inf
                 Inf Inf Inf Inf Inf Inf Inf Inf Inf 0]
    rna_dict = rna(RecurrenceMatrix(adjmat11))
    @test !isinf(rna_dict[:averagepath])
    @test rna_dict[:averagepath] ≈ sum(dismat11[1:9,1:9]) / (9*10)
end


text:
#=

   The Euro coin problem.

   From Think Bayes, page 33ff
   """
   A statistical statement appeared in "The Guardian" on Friday January 4, 2002:
      When spun on edge 250 times, a Belgian one-euro coin
      came up heads 140 times and tails 110. 'It looks very
      suspicious to me,' said Barry Blight, a statistics lecturer
      at the London School of Economics. 'If the coin were
      unbiased, the chance of getting a result as extreme as
      that would be less than 7%.'

   But do these data give evidence that the coin is biased rather than fair?
   """

  Model 1
  mean(prob): 0.5639764325445531
  Distributions of variable probLt0_5 (num:0)
  0.00000 =>    9896  (0.989600)
  1.00000 =>     104  (0.010400)

  Distributions of variable probGt0_5 (num:0)
  1.00000 =>    9896  (0.989600)
  0.00000 =>     104  (0.010400)

  Model 2
  Summary Statistics
  parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

        prob     0.5496    0.0296     0.0003    0.0026     52.7368    1.0154        6.8901
       heads   139.9668    1.8488     0.0185    0.0332   2384.8100    1.0002      311.5770
   probLt0_5     0.0006    0.0245     0.0002    0.0006   1544.1256    1.0005      201.7410
   probGt0_5     0.9994    0.0245     0.0002    0.0006   1544.1256    1.0005      201.7410

   Cf ~/blog/euro_coin_problem.blog
      ~/psi/euro_coin_problem.blog
      ~/webppl/euro_coin_problem.jl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function euro_coin_problem1()
    n = 250

    # Probability of throwing head
    prob ~ Beta(2,2) # Uniform(0,1)
    coin ~ filldist(flip(prob),n)

    sum250 ~ Dirac(sum([coin[i] == true ? 1 : 0 for i in 1:n]))
    true ~ Dirac(sum250 == 140)

    probLt0_5 ~ Dirac(prob < 0.5)
    probGt0_5 ~ Dirac(prob > 0.5)

end

# Simpler and faster
@model function euro_coin_problem2()
    prob ~ Beta(2,2)
    heads ~ Binomial(250,prob)

    true ~ Dirac(heads==140)

    probLt0_5 ~ Dirac(prob < 0.5)
    probGt0_5 ~ Dirac(prob > 0.5)

end

println("Model 1")
model = euro_coin_problem1()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# display(chns)
# display(plot(chns))

println("mean(prob): ", mean(chns[:prob]))
show_var_dist_pct(chns, :probLt0_5)
show_var_dist_pct(chns, :probGt0_5)

println("\nModel 2")
model = euro_coin_problem2()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :probLt0_5)
show_var_dist_pct(chns, :probGt0_5)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

function relaxation(distance::AbstractArray, auxdistances::AbstractArray, cutoff::Real)
  # patterns enabled in the training image
  enabled = !isinf(distance); npatterns = sum(enabled)

  # candidates with good overlap
  dbsize = all(distance[enabled] .== 0) ? npatterns : ceil(Int, cutoff*npatterns)
  overlapdb = selectperm(distance[:], 1:dbsize)

  # candidates in accordance with auxiliary data
  naux = length(auxdistances)
  softdb = fill(Int[], naux)

  patterndb = []
  softdistance = [copy(auxdistances[i]) for i=1:naux]
  frac = .1 * (dbsize / npatterns)
  while true
    softdbsize = ceil(Int, frac*npatterns)

    patterndb = overlapdb
    for n=1:length(softdistance)
      softdistance[n][softdb[n]] = Inf
      softdb[n] = [softdb[n]; selectperm(softdistance[n][:], 1:softdbsize-length(softdb[n]))]

      patterndb = quick_intersect(patterndb, softdb[n], length(distance))

      isempty(patterndb) && break
    end

    !isempty(patterndb) && break
    frac = min(frac + .1, 1)
  end

  patterndb
end

function quick_intersect(A::AbstractVector{Int}, B::AbstractVector{Int},
                         nbits::Integer)
  bitsA = falses(nbits)
  bitsB = falses(nbits)
  bitsA[A] = true
  bitsB[B] = true

  find(bitsA & bitsB)
end