text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
Either `Rational` or `Euclidean`, used for [`Spline`](@ref)s and [`SplineSurface`](@ref)s.
"""
abstract type CurveType end
abstract type Rational <: CurveType end
abstract type Euclidean <: CurveType end
export CurveType

"""
    Spline{P<:CurveType,S<:Number,N,M}

`M` is the curve order, i.e., the highest power of the parameterizing variable, u.
`P` determines the [`CurveType`](@ref).

All Spline types must implement:
```julia
point(curve,u)
```
and have field `controlpolygon`
"""
abstract type Spline{P<:CurveType,S<:Number,N,M} end

"""
    SplineSurface{P,S,N,M} <: ParametricSurface{S,N}

Curve order, `M`, is the same in the u and v direction and fixed over all spans.
`P` determines the [`CurveType`](@ref).
"""
abstract type SplineSurface{P<:CurveType,S,N,M} <: ParametricSurface{S,N} end
export Spline, SplineSurface

function toeuclidean(point)
    return SVector{length(point) - 1}([point[i] / point[end] for i in 1:(lastindex(point) - 1)])
end

function euclideancontrolpoints(curve::Spline{Rational,S,N,M}) where {S,N,M}
    return toeuclidean.(curve.controlpolygon)
end

function euclideancontrolpoints(curve::Spline{Euclidean,S,N,M}) where {S,N,M}
    return curve.controlpolygon
end

#Bernstein polynomial functions
#B(i,n,u) is broken into Bc,Bu because these two functions will be needed when constructing the matrics to compute moving lines or moving planes
Bc(i, n) = factorial(n) ÷ (factorial(i) * factorial(n - i)) #coefficient of the i,n Bernstein basis polynomial
Bu(i, n, u) = u^i * (1 - u)^(n - i) #polynomial part of Bernstein basis polynomial
B(i, n, u) = Bc(i, n) * Bu(i, n, u)


text:
using SDFResults
using Test
using Unitful
using PICDataStructures
using RecursiveArrayTools: recursive_bottom_eltype

@testset "SDFResults.jl" begin
    dir = "gauss"
    sim = read_simulation(dir)
    @test sim isa EPOCHSimulation
    @test isconcretetype(typeof(sim))
    file = sim[1]
    @test file isa SDFFile
    @test size(sim) == (1,)
    @test length(sim) == 1

    Ex, Ey = file[:ex, :ey]
    Ez = file[:ez]
    @test Ex isa ScalarField{3}
    @test Ey isa ScalarField{3}
    @test Ez isa ScalarField{3}
    @test unit(eltype(Ex)) == u"V/m"
    @test unit(recursive_bottom_eltype(getdomain(Ex))) == u"m"
    @test nameof(Ex) == "Electric Field/Ex"
    @test nameof(Ey) == "Electric Field/Ey"
    @test nameof(Ez) == "Electric Field/Ez"
    @test axisnames(getdomain(Ex)) == ["x (m)", "y (m)", "z (m)"]

    vars = (:grid, Symbol("py/electron"), Symbol("pz/electron"))
    @task all(vars .∈ keys(file))
    (x,y,z), py, pz = read(file, vars...)
    @test all(unit.(x) .== u"m")
    @test all(unit.(py) .== u"kg*m/s")

    # test for different code paths in expensive grid detection
    @testset begin
        pz = file["pz/electron"]
        px, py = file["px/electron", "py/electron"]
        @test getdomain(px) == getdomain(py) == getdomain(pz)
        @test length(sim.particle_cache) == 1
    end

    t = get_time(file)
    @test (t |> u"fs") ≈ 10u"fs" atol = 0.1u"fs"

    nx, ny, nz = get_parameter.((file,), (:nx, :ny, :nz))
    @test nx == ny == nz == 10

    λ = get_parameter(file, :laser, :lambda)
    @test (λ |> u"nm") ≈ 800u"nm"

    @test ndims(file) == 3
    @test cell_volume(file) ≠ 0
    @test cell_length(file, :x) != 0
    @test cell_length(file, :y) == cell_length(file, :z)

    # TODO: check that the value is correct
    @test timestep(sim) ≈ 0.36590830829382u"fs"

    @testset "show" begin
        io = IOBuffer()

        show(io, MIME"text/plain"(), sim)
        @test startswith(String(take!(io)), "Epoch3d simulation with 1 files from " *
        "10.1 fs to 10.1 fs.")

        show(io, MIME"text/plain"(), sim[1])
        @test String(take!(io)) == "SDFFile with 44 entries at t = 10.1 fs"
    end
end

ext: jl


text:
"""
    Carries training information for a given tree node
"""
mutable struct TrainNodeGPU{T<:AbstractFloat}
    gain::T
    𝑖::Union{Nothing, AbstractVector{UInt32}}
    ∑::AbstractVector{T}
    h::AbstractArray{T,3}
    hL::AbstractArray{T,3}
    hR::AbstractArray{T,3}
    gains::AbstractMatrix{T}
end

function TrainNodeGPU(nvars, nbins, K, T)
    node = TrainNodeGPU{T}(
            zero(T),
            nothing,
            CUDA.zeros(T, 2*K+1),
            CUDA.zeros(T, (2*K+1) , nbins, nvars),
            CUDA.zeros(T, (2*K+1) , nbins, nvars),
            CUDA.zeros(T, (2*K+1) , nbins, nvars),
            CUDA.zeros(T, nbins, nvars))
    return node
end

struct TreeNodeGPU{T<:AbstractFloat, S, B<:Bool}
    left::S
    right::S
    feat::S
    cond::T
    gain::T
    pred::Vector{T}
    split::B
end

TreeNodeGPU(left::S, right::S, feat::S, cond::T, gain::T, K) where {T<:AbstractFloat, S} = TreeNodeGPU(left, right, feat, cond, gain, zeros(T,K), true)
TreeNodeGPU(pred::Vector{T}) where {T} = TreeNodeGPU(UInt32(0), UInt32(0), UInt32(0), zero(T), zero(T), pred, false)

"""
    single tree is made of a a vector of nodes
"""
struct TreeGPU{T<:AbstractFloat}
    feat::CuVector{Int}
    cond_bin::CuVector{UInt8}
    cond_float::CuVector{T}
    gain::CuVector{T}
    pred::CuMatrix{T}
    split::CuVector{Bool}
end

TreeGPU(x::CuVector{T}) where {T <: AbstractFloat} = TreeGPU(CUDA.zeros(Int, 1), CUDA.zeros(UInt8, 1), CUDA.zeros(T, 1), CUDA.zeros(T, 1), reshape(x, :, 1), CUDA.zeros(Bool, 1))
TreeGPU(depth, K, ::T) where {T <: AbstractFloat} = TreeGPU(CUDA.zeros(Int, 2^depth-1), CUDA.zeros(UInt8, 2^depth-1), CUDA.zeros(T, 2^depth-1), CUDA.zeros(T, 2^depth-1), CUDA.zeros(T, K, 2^depth-1), CUDA.zeros(Bool, 2^depth-1))

# gradient-boosted tree is formed by a vector of trees
struct GBTreeGPU{T<:AbstractFloat}
    trees::Vector{TreeGPU{T}}
    params::EvoTypes
    metric::Metric
    K::Int
    levels
end


text:
"""
Calculates Bézout coefficients (see `gcdx`)
"""
function bezout(a::R, b::R) where {R}
    rev = a < b
    x, y = rev ? (a,b) : (b,a)

    s0, s1 = oneunit(R), zero(R)
    t0, t1 = zero(R), oneunit(R)

    while y != zero(R)
        q = div(x, y)
        x, y = y, x - y * q
        s0, s1 = s1, s0 - q * s1
        t0, t1 = t1, t0 - q * t1
    end

    s, t = rev ? (s0, t0) : (t0, s0)
    g = x

    if g == a
        s = one(R)
        t = zero(R)
    elseif g == -a
        s = -one(R)
        t = zero(R)
    end

    return s, t, g
end

max_stars_repo_name: kevseryolcu/PolyhedralOmega

max_issues_repo_name: kevseryolcu/PolyhedralOmega

max_forks_repo_name: kevseryolcu/PolyhedralOmega


text:
#module ValidatedNumericsTests

using ValidatedNumerics
using FactCheck
using Compat

FactCheck.roughly(a::Interval) = b -> (dist(a, b) < 2*max(eps(a), eps(b)))
roughly = FactCheck.roughly

# Interval tests:

include("interval_tests/intervals.jl")
include("multidim_tests/multidim.jl")
include("decoration_tests/decoration_tests.jl")

# Display tests:
include("display_tests/display.jl")

# Root-finding tests:

include("root_finding_tests/root_finding.jl")

# ITF1788 tests

include("ITF1788_tests/ITF1788_tests.jl")

FactCheck.exitstatus()

#end


text:
export GLMLearner, update!, grad!, predict!, loss, predict, GLMNetLearner, linpred, linpred!

include("glmmodel.jl")

abstract AbstractGLMLearner <: Learner

predict!(obj::AbstractGLMLearner, pr::DenseVector{Float64}, x::DSMat{Float64}; offset=emptyvector(Float64)) =
    predict!(obj.m, pr, obj.coefs, x, offset=offset)

predict(obj::AbstractGLMLearner, x::DSMat{Float64}; offset=emptyvector(Float64)) =
    predict!(obj, Array(Float64, size(x, 1)), x, offset=offset)

linpred!(obj::AbstractGLMLearner, pr::DenseVector{Float64}, x::DSMat{Float64}; offset=emptyvector(Float64)) =
    linpred!(pr, obj.coefs, x, offset=offset)

linpred(obj::AbstractGLMLearner, x::DSMat{Float64}; offset=emptyvector(Float64)) =
    linpred!(obj, Array(Float64, size(x, 1)), x, offset=offset)

loss(obj::AbstractGLMLearner, x::DSMat{Float64}, y::Vector{Float64}; offset=emptyvector(Float64)) =
    loss(obj.m, predict(obj, x, offset=offset), y)

################################
## GLM without regularization ##
################################
type GLMLearner <: AbstractGLMLearner
    m::GLMModel
    p::Int
    coefs::Vector{Float64}
    gr::Vector{Float64}
    optimizer::AbstractSGD
    initialized::Bool
    function GLMLearner(m::GLMModel, optimizer::AbstractSGD)
        obj = new()
        obj.m = m
        obj.optimizer = optimizer
        obj.initialized = false
        obj
    end
end

function Base.show(io::IO, obj::GLMLearner)
    print(io, "Model: ")
    show(io, obj.m)
    print(io, " Optimizer: ")
    show(io, obj.optimizer)
end

#Allocate obj.coefs and obj.gr on first call to update!
function init!(obj::GLMLearner, p)
    obj.initialized && error("already initialized")
    obj.p = p
    obj.coefs = zeros(p)
    obj.gr = Array(Float64, p)
    obj.initialized = true
    obj
end

function update!(obj::GLMLearner, x::DSMat{Float64}, y::Vector{Float64}; offset=emptyvector(Float64), scratch=emptyvector(Float64))
    obj.initialized || init!(obj, size(x, 2))
    if length(scratch) == length(y)
        grad_scratch!(obj.m, obj.gr, which_weights(obj.optimizer, obj.coefs), x, y, scratch, offset=offset)
    else
        grad!(obj.m, obj.gr, which_weights(obj.optimizer, obj.coefs), x, y, offset=offset)
    end
    update!(obj.optimizer, obj.coefs, obj.gr)
    obj
end

#############################
## GLM with regularization ##
#############################

type GLMNetLearner <: AbstractGLMLearner
    m::GLMModel
    lambda1::Float64
    lambda2::Float64
    p::Int
    coefs::Vector{Float64}
    gr::Vector{Float64}
    gr2::Vector{Float64}
    mu::Vector{Float64}
    nu::Vector{Float64}
    optimizer::AbstractSGD
    optimizer2::AbstractSGD
    initialized::Bool
    function GLMNetLearner(m::GLMModel, optimizer::AbstractSGD, lambda1::Float64, lambda2::Float64)
        lambda1 < 0.0 || lambda2 < 0.0 && error("lambda1 and lambda2 should be non-negative")
        obj = new()
        obj.m = m
        obj.lambda1 = lambda1
        obj.lambda2 = lambda2
        obj.optimizer = optimizer
        obj.initialized = false
        obj
    end
end

GLMNetLearner(m::GLMModel, optimizer::AbstractSGD, lambda1 = 0.0, lambda2 = 0.0) = GLMNetLearner(m, optimizer, lambda1, lambda2)

function init!(obj::GLMNetLearner, p)
    obj.initialized && error("already initialized")
    obj.p = p
    obj.coefs = zeros(p)
    obj.gr = Array(Float64, p)
    if obj.lambda1 > 0.0
        obj.mu = zeros(p)
        obj.nu = zeros(p)
        obj.gr2 = Array(Float64, p)
        obj.optimizer2 = deepcopy(obj.optimizer)
    end
    obj.initialized = true
    obj
end

function update!(obj::GLMNetLearner, x::DSMat{Float64}, y::Vector{Float64})
    obj.initialized || init!(obj, size(x,2))

    #calculate the gradient like usual
    grad!(obj.m, obj.gr,  which_weights(obj.optimizer, obj.coefs), x, y)
    if obj.lambda2 > 0.0
        #for l2 regularization, add l2 penalty
        fma!(obj.gr, obj.coefs, obj.lambda2) #gr = gr + lambda2 * coef
    end

    if obj.lambda1 > 0.0
        #for l1 regularization, calculate obj.gr - lambda1, stored in obj.gr2
        map!(Subtract(), obj.gr2, obj.lambda1, obj.gr)
        # then store obj.gr + lambda1 in obj.gr
        add!(obj.gr, obj.lambda1)
        update!(obj.optimizer, obj.mu, obj.gr) #update weights for both grads in mu, nu
        update!(obj.optimizer2, obj.nu, obj.gr2)
        map1!(MaxFun(), obj.mu, 0.0) #trunc temp weights at 0.0
        map1!(MaxFun(), obj.nu, 0.0)
        map!(Subtract(), obj.coefs, obj.mu, obj.nu) #calc final weights as mu - nu
    else
        #if no l1 regularization (so either only l2, or no reg at all)
        #just do standard update
        update!(obj.optimizer, obj.coefs, obj.gr)
    end
    obj
end


text:
#=
#rigi according to 2016 paper
y=679520.05
x=212273.44
#using BenchmarkTools
@btime latLong($x,$y)

@code_warntype latLong(x,y)
lat,long=latLong(x,y)
latA,longA=approxLatLong(x,y)
lat_seconds=((lat-47)*60-3)*60
long_seconds=((long-8)*60-29)*60
latA_seconds=((latA-47)*60-3)*60
longA_seconds=((longA-8)*60-29)*60
-3*60

latLong(x,y)
a,b=latLong(x,y)
a,b=latLong(200_000,600_000)

a,b=latLong(-100_000,100_000)
compound(a),compound(b)

a,b=latLong(191775.030,602030.680)
compound(a),compound(b)
#we have CH1903+ ellipsoid coordinates
X,Y,Z=geoCentricCartesianCoordinates(a,b,897.361)
#X,Y,Z=geoCentricCartesianCoordinates(a,b,500)

#convert to ETRS89
Xnew=X+674.374
Ynew=Y+15.056
Znew=Z+405.346

aNew,bNew,hNew=geoCentricCartesianToEllipsoidCoordinates(Xnew,Ynew,Znew)
compound(aNew),compound(bNew)

@code_warntype geoCentricCartesianToEllipsoidCoordinates(Xnew,Ynew,Znew)
@btime geoCentricCartesianToEllipsoidCoordinates($Xnew,$Ynew,$Znew)

xx=191775.030
yy=602030.680
@btime LV03toETRS89($xx,$yy)
@code_warntype LV03toETRS89(xx,yy)

=#

#=
#Bern, Sidlerstrasse 5
#approxCHCoord(lat,long)
#http://geodesy.geo.admin.ch/reframe/lv03towgs84?easting=600000&northing=200000&altitude=550.0&format=json

lat,long = 46.95108288705888,7.438632502714563
chCoord(lat,long)

y,x=600000,200000
lat_res,long_res=approxLatLong(x,y)
lat-lat_res
long-long_res

approxLatLong(200_000,600_000)

=#

#=
    @btime approxCHCoord($lat,$long)
    @code_warntype approxCHCoord(lat,long)

    x,y=approxCHCoord(lat,long)
    approxLatLong(x,y)

    @btime approxLatLong($x,$y)
    @code_warntype approxLatLong(x,y)
=#

#=
#Rigi
#47°03'28.95659233"N 8°29'11.11127154"E
lat=47+(3*60+28.95659233)/3600
long=8+(29*60+11.11127154)/3600
#y,x=679602.034,212422.148
x,y=chCoord(lat,long)
lat_res,long_res=approxLatLong(x,y)
lat-lat_res
long-long_res
x-chCoord(lat_res,long_res)[1]
y-chCoord(lat_res,long_res)[2]
=#
#=
    @code_warntype chCoord(lat,long)
    @btime chCoord($lat,$long) #198.879ns (199.165,199.165)
=#


text:
module MInterpolations

immutable LinInterp
    grid::Array
    vals::Array
end

function (points::LinInterp)(x::Real)
    grid = points.grid; vals = points.vals
    lower_index = searchsortedlast(points.grid, x)

    # 補外に対応
    if lower_index == 0
        lower_index = 1
    elseif lower_index == length(grid)
        lower_index = length(grid) - 1
    end

    upper_index = lower_index + 1
    grid_step = grid[upper_index] - grid[lower_index]
    val_step = vals[upper_index] - vals[lower_index]
    iterp_val = vals[lower_index] + (val_step / grid_step) * (x - grid[lower_index])

    return iterp_val
end

end # module

max_stars_repo_name: myuuuuun/MInterpolations

max_issues_repo_name: myuuuuun/MInterpolations

max_forks_repo_name: myuuuuun/MInterpolations


text:
# Create example data
x = LinRange(-2, 2, 40)
y = x.^3 .+ x.^2 .+ x .+ 6
# Plot x and y, with dashed stems ended in a star
stem(x, y, "--p")
# Move the baseline to 5
stem(x, y, baseline = 5)


text:
# Choose the underlying optimisation solver.
# opt_solver = CPLEX.Optimizer
opt_solver = Gurobi.Optimizer

function build_gl_solver(gl_solver_sym::Symbol, λ::Float64, η::Float64, n_iter::Int, solve_all_budgets_at_once::Bool)
  if gl_solver_sym == :budgeted
    return OSSBSubgradientBudgeted(CombinatorialBandits.OSSBSubgradientBudgetedFixedParameterChoice(λ, η), n_iter, solve_all_budgets_at_once)
  elseif gl_solver_sym == :exact
    return OSSBExact()
  elseif gl_solver_sym == :exact_naive
    return OSSBExactNaive(opt_solver)
  elseif gl_solver_sym == :exact_smart
    return OSSBExactSmart()
  else
    error("Unknown option: gl_solver_sym = $(gl_solver_sym)")
  end
end

@testset "Graves-Lai" begin
  @testset "m-sets in dimension $d" for d in [2, 3, 4]
    @testset "Solver: $gl_solver_sym" for gl_solver_sym in [:exact, :exact_smart, :exact_naive, :budgeted]
      @testset "m = $m" for m in 1:floor(Int, d/2)
        # Instance parameters
        θ = Dict{Int, Float64}()
        for i in 1:m
          θ[i] = 2
        end
        for i in (m + 1):d
          θ[i] = 1
        end
        distr = Distribution[Normal(p, 1/2) for p in values(θ)]

        # Compute the theoretical value.
        # TODO: Move this to a new GL solver (the original paper should have the corresponding solution?).
        rewards = sort(collect(values(θ)), rev=true)
        @show rewards
        th_val = sum(1.0 / (rewards[m] - rewards[i]) for i in (m + 1):d)

        # Determine the parameters for this problem.
        λ = 1000.0
        η = 0.5
        n_iter = 15

        if m == 1 # d ∈ [2, 3, 4]
          λ = 1.0
        elseif d == 4 && m == 2
          λ = 0.2
          η = 0.2
          n_iter = 200
        end

        # Start the approximation algorithm.
        s = MSetLPSolver(get_solver())
        i = MSet(distr, m, s)
        gl_solver = build_gl_solver(gl_solver_sym, λ, η, n_iter, true)
        gl = CombinatorialBandits.optimise_graves_lai(i, gl_solver, θ, with_trace=false)
        if gl_solver isa OSSBSubgradientBudgeted
          @test gl.objective ≈ th_val atol=1.0e-0
        else
          @test gl.objective ≈ th_val atol=1.0e-4
        end
      end
    end
  end

  bipartite_target_values = Dict{Int, Float64}(2 => 1.0, 3 => 3.0, 4 => 6.0, 5 => 27.8697)
  # Starting at 4, too naïve solvers take really a lot of time.

  @testset "Bipartite matchings in dimension $n" for n in [2, 3]
    @testset "Solver: $gl_solver_sym" for gl_solver_sym in [:exact, :exact_smart, :exact_naive, :budgeted]
      # Instance parameters
      graph = complete_bipartite_graph(n, n)
      θ = Dict((i, j) => (i == j) ? 2.0 : 1.0 for i in 1:n, j in 1:n)
      distr = Distribution[Normal(θ[i, j], 1/2) for i in 1:n, j in 1:n]

      # Determine the parameters for this problem.
      if n == 2
        λ = 1.0
        η = 0.05
        n_iter = 200
      elseif n == 3
        λ = 0.1
        η = 0.2
        n_iter = 400
      end

      # Start the approximation algorithm.
      s = PerfectBipartiteMatchingLPSolver(get_solver())
      i = UncorrelatedPerfectBipartiteMatching(distr, s)
      gl_solver = build_gl_solver(gl_solver_sym, λ, η, n_iter, true)
      gl = CombinatorialBandits.optimise_graves_lai(i, gl_solver, θ, with_trace=false)
      if gl_solver isa OSSBSubgradientBudgeted
        @test gl.objective ≈ bipartite_target_values[n] atol=1.0e-0
      else
        @test gl.objective ≈ bipartite_target_values[n] atol=1.0e-4
      end
    end
  end
end


text:
include("main.jl")

struct LocalPotential           # This is Dirac potential
    V::Float64                  # Additive potential
    loc::Location               # Location of the unit cell
end

struct AtomsSystem
    μ::Float64                          # Chemical potential
    T::Float64                          # Temperature
    potential::Vector{LocalPotential}   # Local potential
end

function diag_V(potential::Vector{LocalPotential})
    U_ = map(x -> x.V, potential) |> Diagonal
    return U_
end


text:
import DataStructures
import XMLDict

include("MortalityTable.jl")

function open_and_read(path)
    s = open(path) do file
        read(file,String)
    end
end

function getXML(open_file)

    return xml = XMLDict.xml_dict(open_file)

end

# get potentially missing value out of dict
function getVal(dict,key)
    try
        return val = parse(Float64,dict[key])
    catch y
        if isa(y,KeyError)
            return val = missing
        else
            throw(y)
        end
    end
end

struct XTbMLTable
    select::DataStructures.DefaultOrderedDict
    ultimate::DataStructures.DefaultOrderedDict
    d::TableMetaData
end

function parseXTbMLTable(x,path)
    md = x["XTbML"]["ContentClassification"]
    name = get(md,"TableName",nothing)
    id = get(md,"TableIdentity",nothing)
    provider = get(md,"ProviderName",nothing)
    reference = get(md,"TableReference",nothing)
    description = get(md,"TableDescription",nothing)
    comments = get(md,"Comments",nothing)
    source_path = path
    d = TableMetaData(name=name,id=id,provider=provider,reference=reference,
        description=description,comments=comments,source_path=source_path)
    tbl = XTbMLTable(DataStructures.DefaultOrderedDict(missing),DataStructures.DefaultOrderedDict(missing),d)

    if isa(x["XTbML"]["Table"],Vector)
        # for a select and ultimate table, will have multiple tables
        # parsed into a vector of tables
        for ai in x["XTbML"]["Table"][1]["Values"]["Axis"]
            issue_age = parse(Int,ai[:t])
            tbl.select[issue_age] = DataStructures.DefaultOrderedDict(missing)
            for aj in ai["Axis"]["Y"]
                duration = parse(Int,aj[:t])
                rate = getVal(aj,"")
                if !ismissing(rate)
                    tbl.select[issue_age][duration]= rate
                end
            end
        end

        for ai in x["XTbML"]["Table"][2]["Values"]["Axis"]["Y"]
            age = parse(Int,ai[:t])
            rate = getVal(ai,"")
            if !ismissing(rate)
                tbl.ultimate[age]= rate
            end
        end
    elseif isa(x["XTbML"]["Table"],DataStructures.OrderedDict)
        # a table without select period will just have one set of values, which
        # are loaded into the OrderedDict

        for ai in x["XTbML"]["Table"]["Values"]["Axis"]["Y"]
            age = parse(Int,ai[:t])
            rate = getVal(ai,"")
            if !ismissing(rate)
                tbl.ultimate[age]= rate
            end
        end
    else
        error("don't know how to handle table: "  *  name)
    end

    return tbl
end

"""
    qx(table::XTbMLTable, issueAge, duration)
Given a mortality table, an issue age, and a duration, returns the appropriate select or ultimate qx.
"""
function qx(table::XTbMLTable, issueAge::Int, duration::Int)
    if length(table.select) > 0
        s = table.select[issueAge]
        if ismissing(s)
            q = missing
        else
            q = s[duration]
        end
    else
        q = missing
    end

    if ismissing(q)
        q = table.ultimate[issueAge + duration - 1]
    end
    return q
end

"""
    qx(table::XTbMLTable, age)
Given a mortality table and an age returns the appropriate ultimate qx.
"""
function qx(table::XTbMLTable, age)
    return table.ultimate[age]
end

function XTbML_Table_To_MortalityTable(tbl::XTbMLTable)

    return MortalityTable(
    OffsetArray([qx(tbl,issue_age,dur) for issue_age=0:mort_max_issue_age,dur=1:mort_max_dur],-1,0),
    OffsetArray([qx(tbl,age) for age=0:mort_max_issue_age],-1),
    tbl.d)
end

function readXTbML(path)
    x = open_and_read(path) |> getXML
    XTbML_Table_To_MortalityTable(parseXTbMLTable(x,path))
end

# Load Available Tables ###

"""
    tables(dir=nothing)

    Loads the XtbML (the SOA XML data format for mortality tables) stored in the
    given path. If no path is specified, will load the packages in the
    MortalityTables package directory. To see where your system keeps packages,
    run `DEPOT_PATH` from a Julia REPL.
"""
function tables(dir=nothing)
    if isnothing(dir)
        table_dir = joinpath(dirname(pathof(MortalityTables)), "tables", "SOA")
    else
        table_dir = dir
    end

    tables = Dict()
    for (root, dirs, files) in walkdir(table_dir)
        for file in files

            if basename(file)[end-3:end] == ".xml"
                tbl = readXTbML(joinpath(root,file))
                tables[tbl.d.name] = tbl
            end
        end
    end
    return tables
end


text:
begin
# Binary Search

# Julia code to recursive binary search x in arr[].
# If x is present then return its location,
# otherwise return -1
import Base:length,show,AbstractArray
arr = [2, 3, 4, 10, 40]
x =
n = length(arr)
l=1
function binarysearch(arr::AbstractArray{},l,r,x)
    if r >= l
		value=r-l//2
        mid = round(Int,l + (r - l) // 2)
        # If element is present at the middle itself
        if arr[mid] == x
            return mid
        # If element is smaller than mid, then it
        # can only be present in left subarray
		elseif arr[mid] > x
            return binarysearch(arr, l, mid-1, x)
        # Else the element can only be present
        # in right subarray
        else
            return binarysearch(arr, mid + 1, r, x)
		end
    else
        # Element is not present in the array
        return -1
	end
end
global result=-1

result=binarysearch(arr,l, n, x)

println("Result after function execution:",result)
if result == -1
    println("Element:",x," is not present in array",arr)
else
    println("Element:",x," is present at index:", result," in array",arr)
end
end

max_stars_repo_name: Learning-Julia/Algorithms

max_issues_repo_name: Learning-Julia/Algorithms

max_forks_repo_name: Learning-Julia/Algorithms


text:
#=
Methods for slicing and dicing TS objects
=#

@doc """
Get the first `n` observations of a TS object
""" ->
head{V,T}(x::TS{V,T}, n::Int=5) = x[1:n,:]

@doc """
Get the last `n` observations of a TS object
""" ->
tail{V,T}(x::TS{V,T}, n::Int=5) = x[end-n+1:end,:]

@doc """
Get the indexes of all rows in an Array containing NaN values
""" ->
function nanrows(x::Array{Float64}; fun::Function=any)
	@assert fun == any || fun == all "Argument `fun` must be either `any` or `all`"
	cutrows = falses(size(x,1))
	@inbounds for i = 1:size(x,1)
		if fun(isnan(x[i,:]))
			cutrows[i] = true
		end
	end
	return cutrows
end

@doc """
Get the indexes of all columns in an Array containing NaN values
""" ->
function nancols(x::Array{Float64}; fun::Function=any)
	@assert fun == any || fun == all "Argument `fun` must be either `any` or `all`"
	cutcols = falses(size(x,2))
	@inbounds for j = 1:size(x,2)
		if fun(isnan(x[:,j]))
			cutcols[j] = true
		end
	end
	return cutcols
end

nanrows{V,T}(x::TS{V,T}) = ts(nanrows(x.values), x.index, x.fields)
nancols(x::TS) = nancols(x.values)

@doc """
Drop missing (NaN) values from an Array
""" ->
function dropnan(x::Array{Float64}; dim::Int=1, fun::Function=any)
	@assert dim == 1 || dim == 2 || dim == 3 "Argument `dim` must be 1 (rows), 2 (cols), or 3 (both)."
	if dim == 1  # rows only
		return x[!nanrows(x, fun=fun),:]
	elseif dim == 2  # columns only
		return x[:,!nancols(x, fun=fun)]
	elseif dim == 3  # rows and columns
		c = !nancols(x, fun=fun)
		return x[!nanrows(x[:,c], fun=fun),:]
	end
end

@doc """
Drop missing (NaN) values from a TS object
""" ->
function dropnan{V,T}(x::TS{V,T}; dim::Int=1, fun::Function=any)
	@assert dim == 1 || dim == 2 || dim == 3 "Argument `dim` must be 1 (rows), 2 (columns), or 3 (both)."
	if dim == 1
		return x[!nanrows(x.values, fun=fun)]
	elseif dim == 2
		return x[:,!nancols(x.values, fun=fun)]
	elseif dim == 3
		return x[!nanrows(x.values, fun=fun), !nancols(x.values, fun=fun)]
	end
end

function ffill!{Float64}(x::AbstractArray{Float64,1})
	i = findfirst(!isnan(x))
	@inbounds for i = i+1:size(x,1)
		isnan(x[i]) ? x[i] = x[i-1] : nothing
	end
	return x
end

function ffill!{Float64}(x::AbstractArray{Float64,2})
	@inbounds for j = 1:size(x,2)
		x[:,j] = ffill!(x[:,j])
	end
end

function bfill!{Float64}(x::AbstractArray{Float64,1})
	i = findlast(!isnan(x))
	@inbounds for i = i-1:-1:1
		isnan(x[i]) ? x[i] = x[i+1] : nothing
	end
	return x
end

function bfill!{Float64}(x::AbstractArray{Float64,2})
	@inbounds for j = 1:size(x,2)
		x[:,j] = bfill!(x[:,j])
	end
end

function interpolate(x1::Int, x2::Int, y1::Float64, y2::Float64)
	m = (y2-y1)/(x2-x1)
	b = y1 - m*x1
	x = collect(x1:1.0:x2)
	y = m*x + b
	return y
end

function linterp!{Float64}(x::AbstractArray{Float64,1})
	@assert size(x,1) > 3 "Must have 3 or more elements to interpolate."
	isval = !isnan(x)
	if all(isval)
		return x
	end
	@assert sum(isval) > 2 "Must have at least 2 non-missing values to interpolate."
	idx = find(isval)
	a = idx[1]
	b = idx[2]
	@inbounds for i = 2:size(idx,1)
		x[a:b] = interpolate(a, b, x[a], x[b])
		a = b
		b = idx[i+1]
	end
	return x
end

function linterp!{Float64}(x::AbstractArray{Float64,2})
	@inbounds for j = 1:size(x,2)
		x[:,j] = linterp!(x[:,j])
	end
end

# TODO: make this more efficient (learn how to assign and not mutate)
@doc """
Fill missing (NaN) values from a TS object
""" ->
function fillnan{V,T}(x::TS{V,T}, method::Symbol=:ffill)
	c = nancols(x.values)
	if !any(c)
		return x
	end
	v = copy(x.values)
	if method == :ffill
		ffill!(v)
	elseif method == :bfill
		bfill!(v)
	elseif method == :linear
		@inbounds for j = find(c)
			linterp!(v)
		end
	#TODO: elseif method == :spline
	else
		error("Invalid method argument: Must be `:ffill`, `:bfill`, `:linear`, or `:spline`.")
	end
	return ts(v, x.index, x.fields)
end

@doc """
Replace missing (NaN) values from a TS object with filled values.
""" ->
function fillnan!{V,T}(x::TS{V,T}, method::Symbol=:ffill)
    c = nancols(x.values)
    if !any(c)
        return x
    end
    v = x.values
    if method == :ffill
        ffill!(v)
    elseif method == :bfill
        bfill!(v)
    elseif method == :linear
        linterp!(v)
    # elseif method == :spline
    end
    return ts(v, x.index, x.fields)
end


text:
function rearrangePOP2Grid(
    data     :: Array{T},
    grid     :: Symbol;
    has_time :: Bool = true,
) where T

    dim = size(data)
    dim_len = length(dim)

    if has_time
        sdim_len = dim_len - 1
    else
        sdim_len = dim_len
    end

    if ! (sdim_len in (2, 3) )
        throw(ErrorException("I can only handle 2D and 3D data"))
    end

    Nx = dim[1]
    Ny = dim[2]

    # new horizontal / vertical spatial dimensions
    new_hsdim = Array{Any}(undef,0)
    new_vsdim = Array{Any}(undef,0)
    vsdim_colon = Array{Any}(undef,0)

    push!(new_hsdim, Nx)

    if grid == :V || grid == :UV
        push!(new_hsdim, Ny+1)
    else
        push!(new_hsdim, Ny)
    end

    if sdim_len == 3
        Nz = dim[3]
        push!(vsdim_colon, Colon())
        if grid == :W
            push!(new_vsdim, Nz+1)
        else
            push!(new_vsdim, Nz)
        end
    end

    dummy_dim = Array{Any}(undef,0)
    dummy_colon = Array{Any}(undef,0)

    if has_time
        Nt = dim[end]
        push!(dummy_dim, Nt)
        push!(dummy_colon, Colon())
    end

    vsdim_zero_circ = [0 for i=1:sdim_len-2]
    dummy_zero_circ = [0 for i=1:length(dummy_dim)]

    new_data = zeros(Float64, new_hsdim..., new_vsdim..., dummy_dim...)

    if grid == :T
        new_data[1:Nx, 1:Ny, vsdim_colon..., dummy_colon...] = data
    elseif grid == :UV
        new_data[1:Nx, 2:Ny+1, vsdim_colon..., dummy_colon...] = circshift(data, (1, 0, vsdim_zero_circ..., dummy_zero_circ...))
    elseif grid == :U
        new_data[1:Nx, 1:Ny, vsdim_colon..., dummy_colon...] = circshift(data, (1, 0, vsdim_zero_circ..., dummy_zero_circ...))
    elseif grid == :V
        new_data[1:Nx, 2:Ny+1, vsdim_colon..., dummy_colon...] = data
    elseif grid == :W
        new_data[1:Nx, 1:Ny, 1:Nz] = data
    else
        throw(ErrorException("Unknown grid: ", string(grid)))
    end

    permute_dims = [1, 2]
    if sdim_len == 3
        pushfirst!(permute_dims, 3)
    end

    if has_time
        push!(permute_dims, maximum(permute_dims)+1)  # it can be 3 or 4
    end

    return permutedims(new_data, permute_dims)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module Statistics

	using Flux

	"""
		error(model::Chain, x::Vector, y::Vector)

	Calculates the percentage of errored predictions.
	"""
	function error(model::Chain, x::Vector, y::Vector)
		errored = 0
		for (xi, yi) in zip(x, y)
			if yi != (argmax(model(xi))-1)
				errored += 1
			end
		end
		errored / length(y)
	end

	export error

end


text:
using PenPlots
using Test
import Documenter.doctest
import GeometryBasics: Point

import Colors.RGB

function Base.isapprox(path1::Path, path2::Path)
    all(isapprox.(path1.points, path2.points))
end

function Base.isapprox(multipath1::MultiPath, multipath2::MultiPath)
    all(isapprox.(multipath1, multipath2))
end

paths = [
    Path([Point(4, 5), Point(5, 13), Point(2, 10)]),
    Path([Point(1, 2), Point(9, 1), Point(4, 17)]),
]

@testset "Doctests" begin
    doctest(PenPlots)
end

@testset "Geometric Algebra Tests" begin
    @test 4 * paths == [[p * 4 for p in path] for path in paths]

    @test Point(5, 8) * paths == [[p * Point(5, 8) for p in path] for path in paths]

    @test isapprox(degree_rotation(90), frac_rotation(0.25))

    @test isapprox(degree_rotation(90), radian_rotation(pi / 2))

    @test isapprox(
        degree_rotation(90) * paths,
        [[Point(-p[2], p[1]) for p in path] for path in paths],
    )

    @test Point(3, 6) + paths == [[Point(3, 6) + p for p in path] for path in paths]
end

@testset "Plot Construction Tests" begin
    paths1 = [Path([Point(4, 5), Point(5, 6), Point(9, 10)])]

    paths2 = [Path([Point(1, 2), Point(7, 1), Point(9, 1)])]
    paths3 = [Path([Point(2, 3), Point(1, 1), Point(8, 8)])]

    red = RGB(1, 0, 0)

    plot = PenPlot(
        # Explicit layer
        Layer("mylayer", paths1, red),

        # Implicit layer
        paths2,

        # Named implicit layer
        "mylayer3" => paths3,
    )

    to_svg(plot)
end

@testset "Extent tests" begin
    @test extent(paths[1]) == Extent(Point(2, 5), Point(5, 13))

    @test extent(paths) == Extent(Point(1, 1), Point(9, 17))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ cb23b3e4-30ba-11eb-048b-33940fa83826
using ModelingToolkit

# ╔═╡ 23f048f2-30bb-11eb-2d54-9fb1c6aed0a5
using Latexify

# ╔═╡ 0d2468ee-30be-11eb-0d75-ed21f3712b85
using Plots

# ╔═╡ 44e05d48-30ba-11eb-1820-c92042f7fc1f
md"ModelingToolkit.jl is a modeling language for high-performance symbolic-numeric computation in scientific computing and scientific machine learning. It allows for users to give a high-level description of a model for symbolic preprocessing to analyze and enhance the model. ModelingToolkit can automatically generate fast functions for model components like Jacobians and Hessians, along with automatically sparsifying and parallelizing the computations. Automatic transformations, such as index reduction, can be applied to the model to make it easier for numerical solvers to handle.
"

# ╔═╡ f6b027e8-30ba-11eb-3279-cfb2ba45c50b
@variables x y

# ╔═╡ 0a1f16ba-30bb-11eb-2116-218fecf8c8e5
z = x^3 + y

# ╔═╡ 10315162-30bb-11eb-1462-d3fb963a29ec
latexify(z)

# ╔═╡ 0887b66c-30bc-11eb-0353-93ba2f2fcba1
z1 = t + t^2

# ╔═╡ 2f5af0ba-30bc-11eb-21d7-43f31111d22e
@derivatives D'~t

# ╔═╡ 1ccf0eae-30bc-11eb-1d8a-cd7f727bb69b
D(z1)

# ╔═╡ 40072080-30bc-11eb-2b27-fb501cc055c8
dz1=expand_derivatives(D(z1))

# ╔═╡ 55a84814-30bc-11eb-25df-cff76501368c
latexify(z1)

# ╔═╡ 4d6c6c00-30bc-11eb-1da6-216840f18165
latexify(dz1)

# ╔═╡ c710d79e-30bc-11eb-06df-8f322f48cda0
fz1=build_function(z1,t,expression=false)

# ╔═╡ 881a68a6-30bd-11eb-35c2-4dc061e74ce2
t1=[1:20...]

# ╔═╡ 1b06d722-30bd-11eb-0bbc-27e45fde5a46
ft1=fz1.([1:20...])

# ╔═╡ f8f7d7c2-30bd-11eb-13ee-e9e682866368
plot(t1,ft1)

# ╔═╡ 9e94b908-30bb-11eb-19f0-33840780b3ee
build_function(z,[x,y])

# ╔═╡ a41d0c4a-30bb-11eb-1434-4982edcff617

# ╔═╡ 440889ec-30bb-11eb-3c5c-a519f291111a
A = [x^2+y 0 2x
     0     0 2y
     y^2+x 0 0]

# ╔═╡ 4b7b92b2-30bb-11eb-342f-0fb5b3f27ec5
latexify(A)

# ╔═╡ Cell order:
# ╠═44e05d48-30ba-11eb-1820-c92042f7fc1f
# ╠═cb23b3e4-30ba-11eb-048b-33940fa83826
# ╠═23f048f2-30bb-11eb-2d54-9fb1c6aed0a5
# ╠═0d2468ee-30be-11eb-0d75-ed21f3712b85
# ╠═f6b027e8-30ba-11eb-3279-cfb2ba45c50b
# ╠═0a1f16ba-30bb-11eb-2116-218fecf8c8e5
# ╠═10315162-30bb-11eb-1462-d3fb963a29ec
# ╠═0887b66c-30bc-11eb-0353-93ba2f2fcba1
# ╠═2f5af0ba-30bc-11eb-21d7-43f31111d22e
# ╠═1ccf0eae-30bc-11eb-1d8a-cd7f727bb69b
# ╠═40072080-30bc-11eb-2b27-fb501cc055c8
# ╠═55a84814-30bc-11eb-25df-cff76501368c
# ╠═4d6c6c00-30bc-11eb-1da6-216840f18165
# ╠═c710d79e-30bc-11eb-06df-8f322f48cda0
# ╠═881a68a6-30bd-11eb-35c2-4dc061e74ce2
# ╠═1b06d722-30bd-11eb-0bbc-27e45fde5a46
# ╠═f8f7d7c2-30bd-11eb-13ee-e9e682866368
# ╠═9e94b908-30bb-11eb-19f0-33840780b3ee
# ╠═a41d0c4a-30bb-11eb-1434-4982edcff617
# ╠═440889ec-30bb-11eb-3c5c-a519f291111a
# ╠═4b7b92b2-30bb-11eb-342f-0fb5b3f27ec5


text:
using SketchyCGAL
using Test

@testset "SketchyCGAL.jl" begin

end


text:
using Random, LinearAlgebra, StatsBase, Distributions

##FUNCTION data_construction
"""Generate synthetic data
INPUTS
- ℓ             Loss function
- n             Number of samples
- k             Number of relevant features
- p             Number of features/covariates
- SNR           Squareroot of the Signal-to-Noise Ratio
- rho           Correlation coefficient between features
OUTPUTS
- X             Matrix of observed covariates
- Y             Vector of responses
- indices_true  Set of relevant features
- w_true        True regressor"""
function data_construction(ℓ::LossFunction, n::Int, k::Int, p::Int, SNR, ρ)
  indices_true, w_true = generate_estimator(p, k)
  #Auto-regresive AR(1) filter preserving noise power
  X = randn(n, p)
  for i in 2:p
      X[:, i] = X[:, i] + ρ*X[:, i-1]
  end
  # rho = ρ
  # R"""
  # X=matrix(NA,$n,$p)
  # for(i in 1:$n){
  #     X[i,]=arima.sim(n=$p,list(ar=$rho))
  # }
  # """
  # @rget X
  # X = √12 * rand(n, p) - (√12)/2

  Y = generate_output(ℓ, noisy_signal(X, indices_true, w_true, SNR))

  return X, Y, indices_true, w_true
end

function hard_data_construction(ℓ::LossFunction, n::Int, k::Int, p::Int, SNR, ρ)
  indices_true, w_true = generate_estimator(p, k)

  indices_true = collect(1:k)
  w_true = ones(k)./√k #sample(-1:2:1, k)

  #Covariance matrix
  θ = (1/√k + 1/k)/2
  Σ = SparseMatrixCSC{Float64}(I, p, p)
  for i in 1:k
    Σ[i,k+1] = θ
    Σ[k+1,i] = θ
  end

  # perm = shuffle(1:p)
  # indices_true = [findfirst( perm .== i) for i in indices_true]
  # Σ[:,:] = Σ[perm, perm]

  d = MvNormal(full(Σ))

  X = rand(d, n)'

  Y = generate_output(ℓ, noisy_signal(X, indices_true, w_true, SNR))

  # E = randn(n, p)  #Noise
  # for j in 1:p
  #   E[:,j] *= norm(X[:,j])/norm(E[:,j])/SNR*2
  # end

  return X, Y, indices_true, w_true
end

# FUNCTIONS Generate support and regressor w
function generate_estimator(p::Int, k::Int)
  weights = ones(p) ./ p
  indices = sort(sample(1:p, StatsBase.Weights(weights), k, replace=false))
  w = sample(-1:2:1, k)
  return indices, w
end

function noisy_signal(X, indices, w, SNR)
  n = size(X,1)
  S = X[:, indices]*w #Signal
  E = randn(n)  #Noise
  E *= (norm(S)/norm(E))/SNR #Scaling according to SNR
  return S+E
end

#FUNCTIONS Generate Output Y
function generate_output(ℓ::Regression, S)
  return S
end

function generate_output(ℓ::L1SVM, S)
  P = 1. ./( 1 .+ exp.(-S) )
  return 1.0*(P .> .5) - 1.0*(P .< .5)
end

function generate_output(ℓ::LogReg, S)
    n = size(S,1)
    P = 1. ./( 1 .+ exp.(-S) )
    return (-ones(n)).^(rand(n) .> P)
end


text:
abstract type Elastic end
abstract type Acoustic end
"""
2-D/ 3-D elastic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdElastic <: Elastic end
"""
2-D/ 3-D acoustic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdAcoustic <: Acoustic end
"""
2-D/ 3-D viscoacoustic forward modeling using staggered-grid velocity-stress finite-difference formulation.
"""
struct FdtdAcousticVisco end
"""
2-D Linearized forward modeling using a finite-difference simulation of the acoustic wave-equation.
"""
struct FdtdAcousticBorn <: Acoustic end
"""
2-D/ 3-D Linearized forward modeling using a analytical solutions for homogeneous acoustic media.
"""
struct AcousticBorn <: Acoustic end
"""
2-D/ 3-D Linearized forward modeling using a analytical solutions for homogeneous elastic media.
"""
struct ElasticBorn <: Elastic end

# define supersource, source and receiver structs
struct Srcs
    n::Int
end
Srcs() = Srcs(0)
struct SSrcs
    n::Int
end
SSrcs() = SSrcs(0)
struct Recs
    n::Int
end
Recs() = Recs(0)


text:
"
  Name     : c4_18.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 1/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
"

Pkg.add("Gadfly")
using Gadfly
draw(SVG("output.svg", 6inch, 3inch), plot([sin, cos], 0, 25))


text:
@testset "Water tests." begin

    water = Water()

    Ts = collect([273.153:10.0:140.3+273.15]...)
    Ps = collect([[101325.0]..., 5.0E6:19.0E6:100.0E6...])

    rhos_CP = [
        999.843288668736
        999.7022064047571
        998.206531114165
        995.648547727687
        992.2152054029182
        988.0336893152614
        983.1942808217368
        977.7629140511095
        971.7885282590977
        965.3075721956923
        1002.321301004123
        1002.0310410719416
        1000.4389567920352
        997.8207863375475
        994.3540212626696
        990.1600120842329
        985.325295236465
        979.9133791196485
        973.9716718278099
        967.5357377570548
        960.6319955267068
        953.2794593328697
        945.4908642455383
        937.2733731926065
        928.628980933184
    ]

    Cps_CP = [
        4219.434432031926
        4195.15402061601
        4184.048832139037
        4179.819080630115
        4179.415075226967
        4181.343152236771
        4184.954590522909
        4190.068861555105
        4196.75552379441
        4205.208431029631
        4195.6772230196875
        4176.951566752204
        4169.021633960576
        4166.726174432052
        4167.529584370848
        4170.198442209381
        4174.226023012843
        4179.512509552163
        4186.1735172851095
        4194.427743361589
        4204.535630211352
        4216.770290859532
        4231.408130310603
        4248.731120668077
        4269.035936461646
    ]

    ks_CP = [
        0.5556573479629365
        0.57878368629924
        0.598017660579779
        0.6143967504241387
        0.6284896176744261
        0.6406244513552712
        0.6510031477223245
        0.6597606493966989
        0.6669962636640698
        0.6727901194178599
        0.5593479468610951
        0.5819555952701735
        0.6008790707937224
        0.6170777712382046
        0.6310751274302673
        0.6431717005931662
        0.653551885098157
        0.6623392583494603
        0.6696256425276464
        0.6754862275416963
        0.6799871759111179
        0.6831891713047554
        0.6851488181656393
        0.6859189511453694
        0.6855484362606141
    ]

    mus_CP=[
    0.0017915688995613212
    0.0013057868784517098
    0.0010015225439034675
    0.0007971708792954054
    0.0006526919175727183
    0.0005464887396638285
    0.00046601393786597313
    0.0004035315809411959
    0.0003540373865383654
    0.00031416450875558944
    0.0017808395544338116
    0.0013014716443436355
    0.0010000705264838394
    0.0007971003389448608
    0.0006533214417919576
    0.0005474817473247409
    0.00046719513565785306
    0.00040480615876922705
    0.00035535285387594236
    0.0003154916084978094
    0.0002828958170615312
    0.0002559024852849985
    0.00023329557520111158
    0.00021416827628534935
    0.0001978326774638032
    ]

    rhos = Float64[]
    Cps = Float64[]
    ks = Float64[]
    mus = Float64[]

    ips = [
        (101325.0, 273.153)
        (101325.0, 283.153)
        (101325.0, 293.153)
        (101325.0, 303.153)
        (101325.0, 313.153)
        (101325.0, 323.153)
        (101325.0, 333.153)
        (101325.0, 343.153)
        (101325.0, 353.153)
        (101325.0, 363.153)
        (5.0e6, 273.153)
        (5.0e6, 283.153)
        (5.0e6, 293.153)
        (5.0e6, 303.153)
        (5.0e6, 313.153)
        (5.0e6, 323.153)
        (5.0e6, 333.153)
        (5.0e6, 343.153)
        (5.0e6, 353.153)
        (5.0e6, 363.153)
        (5.0e6, 373.153)
        (5.0e6, 383.153)
        (5.0e6, 393.153)
        (5.0e6, 403.153)
        (5.0e6, 413.153)
    ]

    for (P,T) ∈ ips
            rho = DrillMudsThermalProps.rho(P, T, water)
            cp = DrillMudsThermalProps.Cp(P, T, water)
            k = DrillMudsThermalProps.therm_cond(P, T, water)
            mu = DrillMudsThermalProps.visc(P,T,water)

            append!(rhos, rho)
            append!(Cps, cp)
            append!(ks, k)
            append!(mus,mu)
    end

    dev_rho = abs.(rhos_CP .- rhos) ./ rhos_CP

    @testset "Density" begin
        for (dev,ip) in zip(dev_rho,ips)
            @test (dev ≤ 0.01)
        end
    end

    dev_Cp=abs.(Cps_CP.-Cps)./Cps
    @testset "Cp" begin
        for (dev,ip) in zip(dev_Cp,ips)
            @test (dev ≤ 0.05)
        end
    end

    dev_mu=abs.(mus.-mus_CP)./mus_CP
    @testset "Viscosity" begin
        for (dev,ip) in zip(dev_mu,ips)
            @test (dev ≤ 0.07)
        end
    end

    dev_ks=abs.(ks-ks_CP)./ks_CP
    @testset "therm. cond." begin
        for (dev,ip) in zip(dev_ks,ips)
            @test (dev ≤ 0.05)
        end
    end

end


text:
#quaternion initialization and beta
q0 = 1.0
q1 = 0.0
q2 = 0.0
q3 = 0.0
beta = 0.1

function invsqrt(x)
	return 1/sqrt(x)
end

function AHRS_update_IMU( gx,  gy,  gz,  ax,  ay,  az, q0, q1, q2, q3)

	gx = gx * 0.0174533
	gy = gy * 0.0174533
	gz = gz * 0.0174533

	qDot1 = 0.5 * (-q1 * gx - q2 * gy - q3 * gz)
	qDot2 = 0.5 * (q0 * gx + q2 * gz - q3 * gy)
	qDot3 = 0.5 * (q0 * gy - q1 * gz + q3 * gx)
	qDot4 = 0.5 * (q0 * gz + q1 * gy - q2 * gx)

	if !((ax == 0.0) && (ay == 0.0) && (az == 0.0))

		norm = invsqrt(ax * ax + ay * ay + az * az)
		ax = ax * norm
		ay = ay * norm
		az = az * norm

		two_q0 = 2.0 * q0
		two_q1 = 2.0 * q1
		two_q2 = 2.0 * q2
		two_q3 = 2.0 * q3
		four_q0 = 4.0 * q0
		four_q1 = 4.0 * q1
		four_q2 = 4.0 * q2
		eight_q1 = 8.0 * q1
		eight_q2 = 8.0 * q2
		q0q0 = q0 * q0
		q1q1 = q1 * q1
		q2q2 = q2 * q2
		q3q3 = q3 * q3

		s0 = four_q0 * q2q2 + two_q2 * ax + four_q0 * q1q1 - two_q1 * ay
		s1 = four_q1 * q3q3 - two_q3 * ax + 4.0 * q0q0 * q1 - two_q0 * ay - four_q1 + eight_q1 * q1q1 + eight_q1 * q2q2 + four_q1 * az
		s2 = 4.0 * q0q0 * q2 + two_q0 * ax + four_q2 * q3q3 - two_q3 * ay - four_q2 + eight_q2 * q1q1 + eight_q2 * q2q2 + four_q2 * az
		s3 = 4.0 * q1q1 * q3 - two_q1 * ax + 4.0 * q2q2 * q3 - two_q2 * ay
		norm = invsqrt(s0 * s0 + s1 * s1 + s2 * s2 + s3 * s3)
		print(s0," ", s1," ", s2," ", s3, " ", norm, " \n")
		s0 = s0 * norm
		s1 = s1 * norm
		s2 = s2 * norm
		s3 = s3 * norm

		qDot1 = qDot1 - beta * s0
		qDot2 = qDot2 - beta * s1
		qDot3 = qDot3 - beta * s2
		qDot4 = qDot4 - beta * s3
	end

	#print(norm ,"\n")
	#print(s0," ", s1," ", s2," ", s3, " \n")
	#print(qDot1," ", qDot2," ", qDot3," ", qDot4, " \n")
	q0 = q0 +  qDot1 * (1.0 / 512.0)
	q1 = q1 +  qDot2 * (1.0 / 512.0)
	q2 = q2 +  qDot3 * (1.0 / 512.0)
	q3 = q3 +  qDot4 * (1.0 / 512.0)

	norm = invsqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3)
	q0 = q0 * norm
	q1 = q1 * norm
	q2 = q2 * norm
	q3 = q3 * norm
	return q0, q1, q2, q3

end

function AHRS_update( gx,  gy,  gz,  ax,  ay,  az,  mx,  my,  mz, q0, q1, q2, q3)
	# Usa IMU para o caso se a medição do magnetometro ser invalida
	if (mx == 0.0) && (my == 0.0) && (mz == 0.0)
		q0, q1, q2, q3 = AHRS_update_IMU(gx, gy, gz, ax, ay, az, q0, q1, q2, q3)
		return q0, q1, q2, q3
	end
	# De graus/sec pra rad/sec
	gx = gx * 0.0174533;
	gy = gy * 0.0174533;
	gz = gz * 0.0174533;

	# Taxa de variação do quaternion pelo giroscopio
	qDot1 = 0.5 * (-q1 * gx - q2 * gy - q3 * gz)
	qDot2 = 0.5 * (q0 * gx + q2 * gz - q3 * gy)
	qDot3 = 0.5 * (q0 * gy - q1 * gz + q3 * gx)
	qDot4 = 0.5 * (q0 * gz + q1 * gy - q2 * gx)

	if !((ax == 0.0) && (ay == 0.0) && (az == 0.0))

		# Normalizando dados do acelerometro
		norm = invsqrt(ax * ax + ay * ay + az * az)
		ax = ax * norm
		ay = ay * norm
		az = az * norm

		# Normaliza magnetometro
		norm = invsqrt(mx * mx + my * my + mz * mz)
		mx = mx * norm
		my = my * norm
		mz = mz * norm

		# Pra nao repetir calculos
		two_q0mx = 2.0 * q0 * mx
		two_q0my = 2.0 * q0 * my
		two_q0mz = 2.0 * q0 * mz
		two_q1mx = 2.0 * q1 * mx
		two_q0 = 2.0 * q0
		two_q1 = 2.0 * q1
		two_q2 = 2.0 * q2
		two_q3 = 2.0 * q3
		two_q0q2 = 2.0 * q0 * q2
		two_q2q3 = 2.0 * q2 * q3
		q0q0 = q0 * q0
		q0q1 = q0 * q1
		q0q2 = q0 * q2
		q0q3 = q0 * q3
		q1q1 = q1 * q1
		q1q2 = q1 * q2
		q1q3 = q1 * q3
		q2q2 = q2 * q2
		q2q3 = q2 * q3
		q3q3 = q3 * q3

		# compensação da direção do campo magnetico
		hx = mx * q0q0 - two_q0my * q3 + two_q0mz * q2 + mx * q1q1 + two_q1 * my * q2 + two_q1 * mz * q3 - mx * q2q2 - mx * q3q3
		hy = two_q0mx * q3 + my * q0q0 - two_q0mz * q1 + two_q1mx * q2 - my * q1q1 + my * q2q2 + two_q2 * mz * q3 - my * q3q3
		two_bx = sqrt(hx * hx + hy * hy)
		two_bz = -two_q0mx * q2 + two_q0my * q1 + mz * q0q0 + two_q1mx * q3 - mz * q1q1 + two_q2 * my * q3 - mz * q2q2 + mz * q3q3
		four_bx = 2.0 * two_bx
		four_bz = 2.0 * two_bz

		# Gradiente descendente
		s0 = -two_q2 * (2.0 * q1q3 - two_q0q2 - ax) + two_q1 * (2.0 * q0q1 + two_q2q3 - ay) - two_bz * q2 * (two_bx * (0.5 - q2q2 - q3q3) + two_bz * (q1q3 - q0q2) - mx) + (-two_bx * q3 + two_bz * q1) * (two_bx * (q1q2 - q0q3) + two_bz * (q0q1 + q2q3) - my) + two_bx * q2 * (two_bx * (q0q2 + q1q3) + two_bz * (0.5 - q1q1 - q2q2) - mz)
		s1 = two_q3 * (2.0 * q1q3 - two_q0q2 - ax) + two_q0 * (2.0 * q0q1 + two_q2q3 - ay) - 4.0 * q1 * (1 - 2.0 * q1q1 - 2.0 * q2q2 - az) + two_bz * q3 * (two_bx * (0.5 - q2q2 - q3q3) + two_bz * (q1q3 - q0q2) - mx) + (two_bx * q2 + two_bz * q0) * (two_bx * (q1q2 - q0q3) + two_bz * (q0q1 + q2q3) - my) + (two_bx * q3 - four_bz * q1) * (two_bx * (q0q2 + q1q3) + two_bz * (0.5 - q1q1 - q2q2) - mz)
		s2 = -two_q0 * (2.0 * q1q3 - two_q0q2 - ax) + two_q3 * (2.0 * q0q1 + two_q2q3 - ay) - 4.0 * q2 * (1 - 2.0 * q1q1 - 2.0 * q2q2 - az) + (-four_bx * q2 - two_bz * q0) * (two_bx * (0.5 - q2q2 - q3q3) + two_bz * (q1q3 - q0q2) - mx) + (two_bx * q1 + two_bz * q3) * (two_bx * (q1q2 - q0q3) + two_bz * (q0q1 + q2q3) - my) + (two_bx * q0 - four_bz * q2) * (two_bx * (q0q2 + q1q3) + two_bz * (0.5 - q1q1 - q2q2) - mz)
		s3 = two_q1 * (2.0 * q1q3 - two_q0q2 - ax) + two_q2 * (2.0 * q0q1 + two_q2q3 - ay) + (-four_bx * q3 + two_bz * q1) * (two_bx * (0.5 - q2q2 - q3q3) + two_bz * (q1q3 - q0q2) - mx) + (-two_bx * q0 + two_bz * q2) * (two_bx * (q1q2 - q0q3) + two_bz * (q0q1 + q2q3) - my) + two_bx * q1 * (two_bx * (q0q2 + q1q3) + two_bz * (0.5 - q1q1 - q2q2) - mz)

		#Normalizando
		norm = invsqrt(s0 * s0 + s1 * s1 + s2 * s2 + s3 * s3)
		s0 = s0 * norm
		s1 = s1 * norm
		s2 = s2 * norm
		s3 = s3 * norm

		# passo do feedback
		qDot1 = qDot1 - beta * s0
		qDot2 = qDot2 - beta * s1
		qDot3 = qDot3 - beta * s2
		qDot4 = qDot4 - beta * s3
	end

	# aplicando no quaterinon
	q0 = q0 + qDot1 * (1.0 / 512.0)
	q1 = q1 + qDot2 * (1.0 / 512.0)
	q2 = q2 + qDot3 * (1.0 / 512.0)
	q3 = q3 + qDot4 * (1.0 / 512.0)

	# Normalizando
	norm = invsqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3)
	q0 = q0 * norm
	q1 = q1 * norm
	q2 = q2 * norm
	q3 = q3 * norm
	return q0, q1, q2, q3
end

gx = 30.0
gy = 30.0
gz = 30.0
ax = 1.0
ay = 1.0
az = 1.0
mx = 30.0
my = 30.0
mz = 30.0
x = 1
y = 0
print(q0," ", q1," ", q2," ", q3, " \n")

for i = 0:9
    q0, q1, q2, q3 = AHRS_update(gx, gy, gz, ax, ay, az, mx, my, mz, q0, q1, q2, q3)
    roll = atan2(q0*q1 + q2*q3, 0.5 - q1*q1 - q2*q2);
	pitch = asin(-2.0 * (q1*q3 - q0*q2));
	yaw = atan2(q1*q2 + q0*q3, 0.5 - q2*q2 - q3*q3);
	print("  Angles ", roll," ", pitch," ", yaw," ", " \n")
    #print("Quaternion ", q0," ", q1," ", q2," ", q3, " \n")
end


text:
const ARGS_TOKENS = LittleDict{Char, Vector{Pair{TokenFinder, Symbol}}}(
    '\\' => [
        forward_match("\\\"") => :SKIP,
    ],
    '"' => [
        forward_match("\"", ['"'], false)     => :SINGLE_QUOTE,
        forward_match("\"\"\"", ['"'], false) => :TRIPLE_QUOTE
    ]
)


text:
# Just dumping partially written ideas here

# function reverse_propagate_delay!(solver,cache,schedule,delay_vec)
#     buffer = zeros(nv(schedule))
#     for v in reverse(topological_sort_by_dfs(get_graph(schedule)))
#         Δt_min = get_path_spec(schedule,v).min_duration
#         buffer[v] = (cache.tF[v] - (cache.t0[v] + Δt_min))
#         for v2 in outneighbors(schedule,v)
#             delay_vec[v] = max(delay_vec[v], delay_vec[v2] - buffer[v2])
#             delay_vec[v] = max(0, delay_vec[v] - (cache.t0[v2] - cache.tF[v]))
#         end
#     end
#     delay_vec
# end
# function get_delay_vec(solver,cache,schedule,v0)
#     delay = zeros(nv(schedule))
#     for v in reverse(topological_sort_by_dfs(get_graph(schedule)))
#         Δt_min = get_path_spec(schedule,v).min_duration
#         delay[v] = (cache.tF[v] - (cache.t0[v] + Δt_min + minimum(cache.local_slack[v])))
#     end
#     # zero out all vtxs that are not predecessors of v
#     delay_vec = zeros(Int,nv(schedule))
#     for v in map(e->e.dst,collect(edges(bfs_tree(get_graph(schedule),v0;dir=:in))))
#         if cache.max_deadline[v] > 0
#             delay_vec[v] = min(delay[v],typemax(Int))
#         end
#     end
#     delay_vec
# end
# function backtrack_deadlines(solver,cache,schedule,v)
#     frontier = Set{Int}([v])
#     delay_cut = Set{Int}()
#     while length(frontier) > 0
#         v = pop!(frontier)
#         Δt_min = get_path_spec(schedule,v).min_duration
#         buffer = (cache.tF[v] - (cache.t0[v] + Δt_min))
#         if get_path_spec(schedule,v).fixed
#             continue
#         elseif cache.max_deadline[v] > 0 # still has room for some delay
#             push!(delay_cut,v)
#         elseif indegree(schedule,v) == 0
#             # if v is a root_node, the deadlines cannot be tightened anymore
#             @log_info(-1,solver.l3_verbosity,"ISPS: deadlines cannot be tightened any more.")
#             return Set{Int}()
#         else
#             for v2 in inneighbors(schedule,v)
#                 push!(frontier,v2)
#             end
#         end
#     end
#     return delay_cut
# end
# function tighten_deadline!(solver,env,route_plan,v,dt=1)
#     cache = get_cache(env)
#     schedule = get_schedule(env)
#     active_set = cache.active_set
#     closed_set = cache.closed_set
#     node_queue = cache.node_queue
#     graph = get_graph(schedule)
#     # adjust allowable_slack
#     cache.max_deadline[v] = max(0, cache.max_deadline[v]-1)
#     # remove from closed set
#     push!(active_set, v)
#     setdiff!(closed_set,v)
#     # Trim schedule at v.t0
#     agent_id = get_path_spec(schedule,v).agent_id
#     if agent_id != -1
#         cbs_env = typeof(env.env)(
#             graph = env.env.graph,
#             agent_idx = agent_id,
#             cost_model = get_cost_model(env.env),
#             heuristic = get_heuristic_model(env.env)
#         )
#         new_path = trim_path(cbs_env,get_paths(route_plan)[agent_id],cache.t0[v])
#         set_solution_path!(route_plan,new_path,agent_id)
#         set_path_cost!(route_plan,new_path.cost,agent_id)
#     end
#     # Update node queue
#     for v2 in map(e->e.dst,collect(edges(bfs_tree(graph,v;dir=:out))))
#         setdiff!(closed_set,v2)
#         setdiff!(active_set,v2)
#     end
#     empty!(node_queue)
#     for v2 in active_set
#         node_queue[v2] = isps_queue_cost(schedule,cache,v2)
#     end
#     env,route_plan
# end


text:
using Constructions
using Documenter

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

makedocs(;
    modules=[Constructions],
    authors="Andreas Tell <atell@soundtheory.com> and contributors",
    repo="https://github.com/ATell-SoundTheory/Constructions.jl/blob/{commit}{path}#{line}",
    sitename="Constructions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ATell-SoundTheory.github.io/Constructions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ATell-SoundTheory/Constructions.jl",
    devbranch="main",
)


text:
include("protobuf/OSMPBF.jl")

struct BBox
    bottom_lat::Float64
    left_lon::Float64
    top_lat::Float64
    right_lon::Float64
end

struct LatLon
    lat::Float64
    lon::Float64
end

struct Node
    latlon::LatLon
    tags::Union{Dict{Symbol,String}, Nothing}
end

struct Way
    refs::Vector{Int64}
    tags::Union{Dict{Symbol,String}, Nothing}
end

struct Relation
    refs::Vector{Int64}
    types::Vector{Symbol}
    roles::Vector{String}
    tags::Union{Dict{Symbol,String}, Nothing}
end

struct Map
    nodes::Dict{Int64,Node}
    ways::Dict{Int64,Way}
    relations::Dict{Int64,Relation}
    meta::Dict{Symbol,Any}

    Map() = new(Dict(), Dict(), Dict(), Dict())
end


text:
module GameZero
using Colors
using Random

export Actor, Game, game, draw, schduler, schedule_once, schedule_interval, schedule_unique, unschedule,
        collide, angle, distance, play_music, play_sound, line, clear, rungame
export Keys, MouseButtons, KeyMods
export Line, Rect, Circle

using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

include("keyboard.jl")
include("timer.jl")
include("window.jl")
include("event.jl")
include("resources.jl")
include("screen.jl")
include("actor.jl")

#Magic variables
const HEIGHTSYMBOL = :HEIGHT
const WIDTHSYMBOL = :WIDTH
const BACKSYMBOL = :BACKGROUND

mutable struct Game
    screen::Screen
    location::String
    game_module::Module
    keyboard::Keyboard
    render_function::Function
    update_function::Function
    onkey_function::Function
    onmousedown_function::Function
    onmouseup_function::Function
    onmousemove_function::Function
    Game() = new()
end

# const EVENT_HANDLER_FN = Dict(
#     MOUSEBUTTONDOWN => "on_mouse_down",
#     MOUSEBUTTONUP => "on_mouse_up",
#     MOUSEMOTION => "on_mouse_move",
#     KEYDOWN => "on_key_down",
#     KEYUP => "on_key_up",
#     MUSIC_END => "on_music_end"
# )

const timer = WallTimer()

const game = Ref{Game}()
const playing = Ref{Bool}(false)
const paused = Ref{Bool}(false)

function __init__()

end

function initscreen(gm::Module, name::String)
    h = getifdefined(gm, HEIGHTSYMBOL, 400)
    w = getifdefined(gm, WIDTHSYMBOL, 600)
    color = getifdefined(gm, BACKSYMBOL, ARGB(colorant"white"))
    s = Screen(name, w, h, color)
    clear(s)
    return s
end

getifdefined(m, s, v) = isdefined(m, s) ? getfield(m, s) : v

mainloop(g::Ref{Game}) = mainloop(g[])

function mainloop(g::Game)
    start!(timer)
    while (true)
      #Don't run if game is paused by system (resizing, lost focus, etc)
      while window_paused[] != 0
          _ = pollEvent!()
          sleep(0.5)
      end

      # Handle Events
        errorMsg = ""
        try
            hadEvents = true
            while hadEvents
                e, hadEvents = pollEvent!()
                t = getEventType(e)
                handleEvents!(g, e, t)
            end
        catch e
            rethrow()
        end

      # Render
      #if (debug && debugText) renderFPS(renderer,last_10_frame_times) end
        SDL2.RenderClear(g.screen.renderer)
        Base.invokelatest(g.render_function, g)
        SDL2.RenderPresent(g.screen.renderer)

        dt = elapsed(timer)
      # Don't let the game proceed at fewer than this frames per second. If an
      # update takes too long, allow the game to actually slow, rather than
      # having too big of frames.
        min_fps = 20.0
        max_fps = 60.0
        dt = min(dt/1e9, 1.0 / min_fps)
        start!(timer)
        Base.invokelatest(g.update_function, g, dt)
        tick!(scheduler[])
        if (playing[] == false)
            throw(QuitException())
        end
    end
end

function handleEvents!(g::Game, e, t)
    global playing, paused
    if (t == SDL2.KEYDOWN || t == SDL2.KEYUP)
        handleKeyPress(g::Game, e, t)
    elseif (t == SDL2.MOUSEBUTTONUP || t == SDL2.MOUSEBUTTONDOWN)
        handleMouseClick(g::Game, e, t)
    #TODO elseif (t == SDL2.MOUSEWHEEL); handleMouseScroll(e)
    elseif (t == SDL2.MOUSEMOTION)
        handleMousePan(g::Game, e, t)
    elseif (t == SDL2.QUIT)
        paused[] = playing[] = false
    end
end

function handleKeyPress(g::Game, e, t)
    keySym = getKeySym(e)
    keyMod = getKeyMod(e)
    @debug "Keyboard" keySym, keyMod
    if (t == SDL2.KEYDOWN)
        push!(g.keyboard, keySym)
        Base.invokelatest(g.onkey_function, g, Keys.Key(keySym), keyMod)
    elseif (t == SDL2.KEYUP)
        delete!(g.keyboard, keySym)
    end
    #keyRepeat = (getKeyRepeat(e) != 0)
end

function handleMouseClick(g::Game, e, t)
    button = getMouseButtonClick(e)
    x = getMouseClickX(e)
    y = getMouseClickY(e)
    @debug "Mouse Button" button, x, y
    if (t == SDL2.MOUSEBUTTONUP)
        Base.invokelatest(g.onmouseup_function, g, (x, y), MouseButtons.MouseButton(button))
    elseif (t == SDL2.MOUSEBUTTONDOWN)
        Base.invokelatest(g.onmousedown_function, g, (x, y), MouseButtons.MouseButton(button))
    end
end

function handleMousePan(g::Game, e, t)
    x = getMouseMoveX(e)
    y = getMouseMoveY(e)
    @debug "Mouse Move" x, y
    Base.invokelatest(g.onmousemove_function, g, (x, y))
end

getKeySym(e) = bitcat(UInt32, e[24:-1:21])
getKeyRepeat(e) = bitcat(UInt8, e[14:-1:14])
getKeyMod(e) = bitcat(UInt16, e[26:-1:25])

getMouseButtonClick(e) = bitcat(UInt8, e[17:-1:17])
getMouseClickX(e) =  bitcat(Int32, e[24:-1:21])
getMouseClickY(e) = bitcat(Int32, e[28:-1:25])

getMouseMoveX(e) = bitcat(Int32, e[24:-1:21])
getMouseMoveY(e) = bitcat(Int32, e[28:-1:25])

function rungame(jlf::String)
    global playing, paused
    g = initgame(jlf::String)
    try
        playing[] = paused[] = true
        mainloop(g)
    catch e
        if !isa(e, QuitException) && !isa(e, InterruptException)
            @error e exception = (e, catch_backtrace())
        end
    finally
        GameZero.quitSDL(game[])
    end
end

function initgame(jlf::String)
    if !isfile(jlf)
        ArgumentError("File not found: $jlf")
    end
    name = titlecase(replace(basename(jlf), ".jl"=>""))
    module_name = Symbol(name*"_"*randstring(5))
    game_module = Module(module_name)
    @debug "Initialised Anonymous Game Module" module_name
    initSDL()
    game[] = Game()
    scheduler[] = Scheduler()
    g = game[]
    g.game_module = game_module
    g.location = dirname(jlf)
    g.keyboard = Keyboard()

    Base.include_string(game_module, "using GameZero")
    Base.include_string(game_module, "import GameZero.draw")
    Base.include_string(game_module, "using Colors")
    Base.include(game_module, jlf)

    g.update_function = getfn(game_module, :update, 2)
    g.render_function = getfn(game_module, :draw, 1)
    g.onkey_function = getfn(game_module, :on_key_down, 3)
    g.onmouseup_function = getfn(game_module, :on_mouse_up, 3)
    g.onmousedown_function = getfn(game_module, :on_mouse_down, 3)
    g.onmousemove_function = getfn(game_module, :on_mouse_move, 2)
    g.screen = initscreen(game_module, "GameZero::"*name)
    clear(g.screen)
    return g
end

function getfn(m::Module, s::Symbol, maxargs=3)
    if isdefined(m, s)
        fn = getfield(m, s)

        ms = copy(methods(fn).ms)
        filter!(x->x.module == m, ms)
        if length(ms) > 1
            sort!(ms, by=x->x.nargs, rev=true)
        end
        m = ms[1]
        if (m.nargs - 1) > maxargs
            error("Found a $s function with $(m.nargs-1) arguments. A maximum of $maxargs arguments are allowed.")
        end
        @debug "Event method" fn m.nargs
        #TODO Validate types for arguments
        if m.nargs - 1 == maxargs #required to handle the zero-arg case
            return fn
        end
        return (x...) -> fn(x[1:(m.nargs-1)]...)
    else
        return (x...) -> nothing
    end
end

# Having a QuitException is useful for testing, since an exception will simply
# pause the interpreter. For release builds, the catch() block will call quitSDL().
struct QuitException <: Exception end

function getSDLError()
    x = SDL2.GetError()
    return unsafe_string(x)
end

function initSDL()
    SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLEBUFFERS, 4)
    SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLESAMPLES, 4)
    r = SDL2.Init(UInt32(SDL2.INIT_VIDEO | SDL2.INIT_AUDIO))
    if r != 0
        error("Uanble to initialise SDL: $(getSDLError())")
    end
    SDL2.TTF_Init()

    mix_init_flags = SDL2.MIX_INIT_FLAC|SDL2.MIX_INIT_MP3|SDL2.MIX_INIT_OGG
    inited = SDL2.Mix_Init(Int32(mix_init_flags))
    if inited & mix_init_flags != mix_init_flags
        @warn "Failed to initialise audio mixer properly. All sounds may not play correctly\n$(getSDLError())"
    end

    device = SDL2.Mix_OpenAudio(Int32(22050), UInt16(SDL2.MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
    if device != 0
        @warn "No audio device available, sounds and music will not play.\n$(getSDLError())"
        SDL2.Mix_CloseAudio()
    end
end

function quitSDL(g::Game)
    # Need to close the callback before quitting SDL to prevent it from hanging
    # https://github.com/n0name/2D_Engine/issues/3
    @debug "Quitting the game"
    clear!(scheduler[])
    SDL2.DelEventWatch(window_event_watcher_cfunc[], g.screen.window);
    SDL2.DestroyRenderer(g.screen.renderer)
    SDL2.DestroyWindow(g.screen.window)
    #Run all finalisers
    GC.gc();GC.gc();
    quitSDL()
end

function quitSDL()
    SDL2.Mix_HaltMusic()
    SDL2.Mix_HaltChannel(Int32(-1))
    SDL2.Mix_CloseAudio()
    SDL2.TTF_Quit()
    SDL2.Mix_Quit()
    SDL2.Quit()
end

function main()
    if length(ARGS) < 1
        throw(ArgumentError("No file to run"))
    end
    jlf = ARGS[1]
    rungame(jlf)
end

end # module


text:
using LinearAlgebra
using Distributions
using Random
using Test
using prml
using prml: pdf

include("ch01.jl")
include("ch02.jl")
include("ch03.jl")
include("ch05.jl")
include("ch06.jl")

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type AutocorLenAlgorithm

Abstract type for integrated autocorrelation length estimation algorithms.
"""
abstract type AutocorLenAlgorithm end
export AutocorLenAlgorithm

"""
    bat_integrated_autocorr_len(
        v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
        algorithm::AutocorLenAlgorithm = GeyerAutocorLen()
    )

Estimate the integrated autocorrelation length of variate series `v`,
separately for each degree of freedom.

Returns a NamedTuple of the shape

```julia
(result = integrated_autocorr_len, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_integrated_autocorr_len`, add methods to
    `bat_integrated_autocorr_len_impl` instead.
"""
function bat_integrated_autocorr_len end
export bat_integrated_autocorr_len

function bat_integrated_autocorr_len_impl end

function bat_integrated_autocorr_len(
    v::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}},
    algorithm = bat_default_withdebug(bat_integrated_autocorr_len, Val(:algorithm), v)
)
    r = bat_integrated_autocorr_len_impl(v, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_integrated_autocorr_len), ::Val{:algorithm}, x::AutocorLenAlgorithm)
    "Using integrated autocorrelation length estimator $x"
end


text:
struct CopositiveInnerVariableBridge{T, VB} <: AbstractVariableBridge
    variable_bridge::VB
    # TODO store GreaterThan constraints
end

function add_variable_bridge(::Type{CopositiveInnerVariableBridge{T, VB}},
                             model::MOI.ModelLike, set::CopositiveInner) where {T, VB}
    func, variable_bridge = add_variable_bridge(VB, model, set.psd_inner)
    F = MOI.ScalarAffineFunction{T}
    # If `func isa F` then `convert` won't create a copy so we explictly copy it.
    # It is less costly to copy it before, e.g. if `q is MOI.SingleVariable`.
    Q = F[convert(F, copy(q)) for q in func]
    k = 0
    for j in 1:side_dimension(set)
        for i in 1:(j-1)
            k += 1
            # TODO: these should be deleted when the bridge is deleted
            Nij = MOI.SingleVariable(MOI.add_variable(model))
            MOI.add_constraint(model, Nij, MOI.GreaterThan(zero(T)))
            MOIU.operate!(+, T, Q[k], Nij)
        end
        k += 1 # for diagonal entry (i, i)
    end
    return Q, CopositiveInnerVariableBridge{T, VB}(variable_bridge)
end

function MOIB.added_constraint_types(::Type{CopositiveInnerVariableBridge{T, VB}}) where {T, VB}
    added = MOIB.added_constraint_types(VB)
    push!(added, (MOI.SingleVariable, MOI.GreaterThan{T}))
    return added
end

function variable_bridge_type(::Type{CopositiveInner{S}}, T::Type) where S
    return CopositiveInnerVariableBridge{T, variable_bridge_type(S, T)}
end

function MOI.delete(model::MOI.ModelLike, bridge::CopositiveInnerVariableBridge)
    # TODO remove GreaterThan constraints
    MOI.delete(model, bridge.variable_bridge)
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{MomentMatrixAttribute, GramMatrixAttribute},
                 bridge::CopositiveInnerVariableBridge)
    return MOI.get(modle, attr, bridge.variable_bridge)
end


text:
module MyLib

export increment32, increment64

function increment(count)
    count += 1
    return count
end

Base.@ccallable function increment32(count::Cint)::Cint
    count = increment(count)
    println("Incremented count: $count (Cint)")
    return count
end

Base.@ccallable function increment64(count::Clong)::Clong
    count = increment(count)
    println("Incremented count: $count (Clong)")
    return count
end

end


text:
export get_Jacobian

"""
$(SIGNATURES)

Obtain the symbolic Jacobian matrix of `eom` (either a `HarmonicEquation` or a `DifferentialEquation`).

"""
function get_Jacobian(eom::HarmonicEquation)

    lhs = Num.(getfield.(HarmonicBalance.rearrange_standard(eom).equations, :lhs))

    old_vars = get_variables(eom)
    bracket_rules = [var => HarmonicBalance.declare_variable(var_name(var)) for var in old_vars]

    lhs = substitute_all(lhs,bracket_rules)
    vars = getindex.(bracket_rules, 2)
    J = Matrix{Num}(undef, length(vars), length(vars))
    for idx in CartesianIndices(J)
        J[idx] = expand_derivatives(d(lhs[idx[1]], vars[idx[2]]))
    end
    J

end

#   Obtain a Jacobian from a `DifferentialEquation` by first converting it into a `HarmonicEquation`.
function get_Jacobian(diff_eom::DifferentialEquation)
    @variables T
    harmonic_eq = get_harmonic_equations(diff_eom, slow_time=T, fast_time=first(get_independent_variables(diff_eom)))
    get_Jacobian(harmonic_eq)
end

"""Convert a set of equations in variables `vars` into matrix form.
If the equations are linear, this matrix does not depends on the variables."""
function equations_to_matrix(eqs::Vector{Num}, vars::Vector{Num})
    length(eqs) != length(vars) ? error("System under- or over-constrained!") : nothing
    M = Matrix{Num}(undef, length(eqs), length(eqs))
    for idx in CartesianIndices(M)
        M[idx] = expand_derivatives(d(eqs[idx[1]], vars[idx[2]]))
    end
    M
end

# get the Jacobian of a set of equations `eqs` with respect to the variables `vars`
function get_Jacobian(eqs::Vector{Num}, vars::Vector{Num})
    length(eqs) == length(vars) || error("Jacobians are only defined for square matrices!")
    M = Matrix{Num}(undef, length(vars), length(vars))

    for idx in CartesianIndices(M)
        M[idx] = expand_derivatives(d(eqs[idx[1]], vars[idx[2]]))
    end
    M
end

get_Jacobian(eqs::Vector{Equation}, vars::Vector{Num}) = get_Jacobian(Num.(getfield.(eqs, :lhs) .- getfield.(eqs, :rhs)), vars)

function get_Jacobian_steady(eom::HarmonicEquation; differential_order=0)
    vars_simp = Dict([var => HarmonicBalance.declare_variable(var_name(var)) for var in get_variables(eom)])
    T = get_independent_variables(eom)[1]
    J = get_Jacobian(eom.equations, d(get_variables(eom), T, differential_order))
    expand_derivatives.(HarmonicBalance.substitute_all(J, vars_simp))
end

# COMPILE THIS!
function get_implicit_Jacobian(eom::HarmonicEquation)
    M = get_Jacobian_steady(eom, differential_order=0)
    Mp = get_Jacobian_steady(eom, differential_order=1)
    function J(soln::Dict)
        -inv(ComplexF64.(substitute_all(Mp, soln))) * ComplexF64.(substitute_all(M, soln))
    end
    J
end

###
# STUFF BELOW IS FOR THE CORRECTED JACOBIAN METHOD
###

"""
    get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)

Obtain the symbolic linear response matrix of a `diff_eq` corresponding to a perturbation frequency `freq`.
This routine cannot accept a `HarmonicEquation` since there, some time-derivatives are already dropped.
`order` denotes the highest differential order to be considered.

"""
function get_response_matrix(diff_eq::DifferentialEquation, freq::Num; order=2)::Matrix
    @variables T, i
    time = get_independent_variables(diff_eq)[1]

    eom = HarmonicBalance.harmonic_ansatz(diff_eq, time)

    # replace the time-dependence of harmonic variables by slow time BUT do not drop any derivatives
    eom = HarmonicBalance.slow_flow(eom, fast_time=time, slow_time=T, degree=order+1)

    eom = HarmonicBalance.fourier_transform(eom, time)

    M = get_Jacobian(eom.equations, get_variables(eom))
    for n in 1:order
        M += (i * freq)^n * get_Jacobian(eom.equations, d(get_variables(eom), T, n))
    end
    M = substitute_all(M, [var => HarmonicBalance.declare_variable(var_name(var)) for var in get_variables(eom)])
    substitute_all(expand_derivatives.(M), i => im)
end


text:
module Distributed

export
    MultiArch, child_architecture, reconstruct_global_grid,
    HaloCommunication, HaloCommunicationBC,
    inject_halo_communication_boundary_conditions,
    DistributedFFTBasedPoissonSolver

using MPI

using Oceananigans.Utils
using Oceananigans.Grids

include("distributed_utils.jl")
include("multi_architectures.jl")
include("partition_assemble.jl")
include("distributed_grids.jl")
include("distributed_kernel_launching.jl")
include("halo_communication_bcs.jl")
include("halo_communication.jl")
include("distributed_apply_flux_bcs.jl")
include("distributed_fields.jl")
include("distributed_fft_based_poisson_solver.jl")

end # module


text:
using HDF5
using Test

@testset "Raw Chunk I/O" begin

fn = tempname()

# Direct chunk write is no longer dependent on HL library
# Test direct chunk writing Cartesian index
h5open(fn, "w") do f
    d = create_dataset(f, "dataset", datatype(Int), dataspace(4, 4), chunk=(2, 2))
    HDF5.API.h5d_extend(d, HDF5.API.hsize_t[3,3]) # should do nothing (deprecated call)
    HDF5.API.h5d_extend(d, HDF5.API.hsize_t[4,4]) # should do nothing (deprecated call)
    raw = HDF5.ChunkStorage(d)
    raw[1,1] = 0, collect(reinterpret(UInt8, [1,2,5,6]))
    raw[3,1] = 0, collect(reinterpret(UInt8, [3,4,7,8]))
    raw[1,3] = 0, collect(reinterpret(UInt8, [9,10,13,14]))
    raw[3,3] = 0, collect(reinterpret(UInt8, [11,12,15,16]))
end

# Test read back
@test h5open(fn, "r") do f
    vec(f["dataset"][:,:])
end == collect(1:16)

# Test reading direct chunks via linear indexing
h5open(fn, "r") do f
    d = f["dataset"]
    raw = HDF5.ChunkStorage{IndexLinear}(d)
    @test size(raw) == (4,)
    @test length(raw) == 4
    @test axes(raw) == (Base.OneTo(4),)
    @test prod(HDF5.get_num_chunks_per_dim(d)) == HDF5.get_num_chunks(d)
    if v"1.10.5" ≤ HDF5.API._libhdf5_build_ver
        @test HDF5.get_chunk_length(d) == HDF5.API.h5d_get_chunk_info(d,1)[:size]
    end
    @test reinterpret(Int, raw[1][2]) == [1,2,5,6]
    @test reinterpret(Int, raw[2][2]) == [3,4,7,8]
    @test reinterpret(Int, raw[3][2]) == [9,10,13,14]
    @test reinterpret(Int, raw[4][2]) == [11,12,15,16]
    # Test 0-based indexed API
    @test HDF5.get_chunk_offset(d, 0) == (0, 0)
    @test HDF5.get_chunk_offset(d, 1) == (2, 0)
    @test HDF5.get_chunk_offset(d, 2) == (0, 2)
    @test HDF5.get_chunk_offset(d, 3) == (2, 2)
    # Test reverse look up of index from coords
    @test HDF5.get_chunk_index(d, (0, 0)) == 0
    @test HDF5.get_chunk_index(d, (2, 0)) == 1
    @test HDF5.get_chunk_index(d, (0, 2)) == 2
    @test HDF5.get_chunk_index(d, (2, 2)) == 3
    # Test internal coordinates
    @test HDF5.get_chunk_index(d, (1, 1)) == 0
    @test HDF5.get_chunk_index(d, (3, 1)) == 1
    @test HDF5.get_chunk_index(d, (1, 3)) == 2
    @test HDF5.get_chunk_index(d, (3, 3)) == 3
end

# Test direct write chunk writing via linear indexing
h5open(fn, "w") do f
    d = create_dataset(f, "dataset", datatype(Int64), dataspace(4, 6), chunk=(2, 3))
    raw = HDF5.ChunkStorage{IndexLinear}(d)
    raw[1] = 0, collect(reinterpret(UInt8, Int64[1,2,5,6, 9,10]))
    raw[2] = 0, collect(reinterpret(UInt8, Int64[3,4,7,8,11,12]))
    raw[3] = 0, collect(reinterpret(UInt8, Int64[13,14,17,18,21,22]))
    raw[4] = 0, collect(reinterpret(UInt8, Int64[15,16,19,20,23,24]))
end

@test h5open(fn, "r") do f
    f["dataset"][:,:]
end == reshape(1:24, 4, 6)

h5open(fn, "r") do f
    d = f["dataset"]
    raw = HDF5.ChunkStorage(d)
    chunk = HDF5.get_chunk(d)
    extent = HDF5.get_extent_dims(d)[1]

    @test chunk == (2, 3)
    @test extent == (4, 6)
    @test size(raw) == (2, 2)
    @test length(raw) == 4
    @test axes(raw) == (1:2:4, 1:3:6)
    @test prod(HDF5.get_num_chunks_per_dim(d)) == HDF5.get_num_chunks(d)

    # Test 0-based indexed API
    @test HDF5.get_chunk_offset(d, 0) == (0, 0)
    @test HDF5.get_chunk_offset(d, 1) == (2, 0)
    @test HDF5.get_chunk_offset(d, 2) == (0, 3)
    @test HDF5.get_chunk_offset(d, 3) == (2, 3)
    # Test reverse look up of index from coords
    @test HDF5.get_chunk_index(d, (0, 0)) == 0
    @test HDF5.get_chunk_index(d, (2, 0)) == 1
    @test HDF5.get_chunk_index(d, (0, 3)) == 2
    @test HDF5.get_chunk_index(d, (2, 3)) == 3
    # Test internal coordinates
    @test HDF5.get_chunk_index(d, (1, 1)) == 0
    @test HDF5.get_chunk_index(d, (3, 1)) == 1
    @test HDF5.get_chunk_index(d, (1, 4)) == 2
    @test HDF5.get_chunk_index(d, (3, 4)) == 3

    if v"1.10.5" ≤ HDF5.API._libhdf5_build_ver
        chunk_length = HDF5.get_chunk_length(d)
        origin = HDF5.API.h5d_get_chunk_info(d, 0)
        @test chunk_length == origin[:size]
        chunk_info = HDF5.API.h5d_get_chunk_info_by_coord(d, HDF5.API.hsize_t[0, 1])
        @test chunk_info[:filter_mask] == 0
        @test chunk_info[:size] == chunk_length

        # Test HDF5.get_chunk_offset equivalence to h5d_get_chunk_info information
        @test all(reverse(HDF5.API.h5d_get_chunk_info(d, 3)[:offset]) .== HDF5.get_chunk_offset(d, 3))

        # Test HDF5.get_chunk_index equivalence to h5d_get_chunk_info_by_coord information
        offset = HDF5.API.hsize_t[2,3]
        chunk_info = HDF5.API.h5d_get_chunk_info_by_coord(d, reverse(offset))
        @test HDF5.get_chunk_index(d, offset) == (chunk_info[:addr] - origin[:addr]) ÷ chunk_info[:size]

        @test HDF5.API.h5d_get_chunk_storage_size(d, HDF5.API.hsize_t[0, 1]) == chunk_length
        @test HDF5.API.h5d_get_storage_size(d) == sizeof(Int64)*24
        @test HDF5.API.h5d_get_space_status(d) == HDF5.API.H5D_SPACE_STATUS_ALLOCATED
    end

    # Manually reconstruct matrix
    A = Matrix{Int}(undef, extent)
    for (r,c) in Iterators.product(axes(raw)...)
        A[r:r+chunk[1]-1, c:c+chunk[2]-1] .= reshape( reinterpret(Int64, raw[r,c][2]), chunk)
    end
    @test A == reshape(1:24, extent)

end

rm(fn)

end # testset "Raw Chunk I/O"


text:
@testset "KSP.jl" begin
    begin
        vec = PetscVec()
        setSize!(vec, n_local=(Int32)(3))
        vec[1:3] = (Float64)[1, 2, 3]
        assemble!(vec)

        sol = PetscVec()
        setSize!(sol, n_local=(Int32)(3))
        assemble!(sol)

        mat = PetscMat()
        setSize!(mat, m_local=(Int32)(3), n_local=(Int32)(3))
        setPreallocation!(mat, (Int32)[1,1,1], (Int32)[])

        mat[1,1] = (Float64)(1)
        mat[2,2] = (Float64)(2)
        mat[3,3] = (Float64)(3)

        assemble!(mat)

        ksp = PetscKSP()

        setOperators(ksp, mat)

        solve!(ksp, vec, sol)

        for val in sol[1:3]
            @test abs(sol[1]-1.0) < 1e-9
        end
    end
end


text:
# =============================================================================
# Compare with original single-stage MIP, 329 pattern example
# =============================================================================
nbedrooms_range = collect(1:2)
nbedrooms_frequency = Weights([0.75, 0.25])
nbathrooms_range = collect(1:2)
nbathrooms_frequency = Weights([0.75; 0.25])
prices_range_pp = [800.0, 1000.0]
area_ranges = [0.0, 700.0, 850.0]
budget = 1000.0
market_data = StudentHousing.MarketData(nbedrooms_range, nbedrooms_frequency,
    nbathrooms_range, nbathrooms_frequency, prices_range_pp, area_ranges)
nhouses = 50
srand(1)
d = StudentHousingData(market_data, nhouses = nhouses, budget = budget,
                  demand_distribution = Uniform(0.0, 100.0))

# =============================================================================
# Solve problem with column generation algorithm, pattern subproblems
# =============================================================================
tic()
m3, U_generated, λ_generated, pattern_choice = solve_ia_p_generation(d)
toc()
# elapsed time: 15.290560206 seconds

# When we solved the same problem with one subproblem for each house, the solving
# time was around 2 seconds.

# =============================================================================
# Look at solution
# =============================================================================
println(getobjectivevalue(m3))
# 25.81671609417941
println("Shortage = ", sum(sum(d.demands)) - getobjectivevalue(m3))
# Shortage = 16615.18328390582
# ... this is equal to the LP relaxation

# When we solved the same problem with one subproblem for each house, the LP
# bound was 16615.366.

# Not integer, so get integer solution with existing columns
r = recover_integer_soln_p(m3, d)
println("Integer shortage = ", sum(sum(d.demands)) - sum(sum(r)))
# Integer shortage = 16616.0
# ... which matches the MIP solution

# =============================================================================
# Solve problem with column generation algorithm, pattern subproblems
# without looking for best reduced cost, stop @ any positive reduced cost
# =============================================================================
# srand(32)
# d = StudentHousingData(market_data, nhouses = 2, budget = budget,
#                   demand_distribution = Uniform(0.0, 100.0))
tic()
m3, U_generated, λ_generated, pattern_choice = solve_ia_p_generation(d, true)
toc()
# too long on PC


text:
## Imports

using Base.Threads
using Flux
using Graphs
using InferOpt
using MultiAgentPathFinding
using PythonCall
using ProgressMeter
using UnicodePlots

## Test

rail_generators = pyimport("flatland.envs.rail_generators")
line_generators = pyimport("flatland.envs.line_generators")
rail_env = pyimport("flatland.envs.rail_env")

rail_generator = rail_generators.sparse_rail_generator(; max_num_cities=4)
line_generator = line_generators.sparse_line_generator()

pyenv = rail_env.RailEnv(;
    width=30,
    height=30,
    number_of_agents=20,
    rail_generator=rail_generator,
    line_generator=line_generator,
    random_seed=11,
)
pyenv.reset();

## Data generation

nb_instances = 10

instances = Vector{typeof(flatland_mapf(pyenv))}(undef, nb_instances);
@showprogress "Generating instances: " for k in 1:nb_instances
    pyenv.reset()
    instances[k] = flatland_mapf(pyenv)
end

solutions = Vector{Solution}(undef, nb_instances);
@threads for k in 1:nb_instances
    @info "Instance $k solved by thread $(threadid())"
    solutions[k] = large_neighborhood_search(instances[k]; N=10, steps=100, progress=false)
end

T = maximum(max_time(solution) for solution in solutions) * 2
X = [agents_embedding(mapf) for mapf in instances];
Y = [solution_to_vec(solution, mapf; T=T) for (solution, mapf) in zip(solutions, instances)];

function maximizer(θ; mapf)
    permutation = sortperm(θ; rev=true)
    solution = cooperative_astar(mapf, permutation)
    ŷ = solution_to_vec(solution, mapf; T=T)
    return ŷ
end

## Initialization

encoder = Chain(Dense(size(X[1], 1), 1), vec)
model = Perturbed(maximizer; ε=0.1, M=10)
squared_loss(ŷ, y) = sum(abs2, y - ŷ) / T;

λ = 30
opt = ADAGrad();
par = Flux.params(encoder)

losses = Float64[]

## Training

for epoch in 1:100
    l = 0.0
    @showprogress "Epoch $epoch/100 - " for k in 1:nb_instances
        gs = gradient(par) do
            l += (
                squared_loss(model(encoder(X[k]); mapf=instances[k]), Y[k]) +
                (λ / nb_instances) * sum(abs, encoder[1].weight)
            )
        end
        Flux.update!(opt, par, gs)
    end
    push!(losses, l)
end;

println(lineplot(losses))

## Evaluation

solutions_naive = Vector{Solution}(undef, nb_instances);
solutions_pred = Vector{Solution}(undef, nb_instances);
for k in 1:nb_instances
    solutions_naive[k] = cooperative_astar(instances[k], 1:nb_agents(instances[k]))
    θ = encoder(X[k])
    permutation = sortperm(θ; rev=true)
    solutions_pred[k] = cooperative_astar(instances[k], permutation)
end

avg_flowtime =
    sum(flowtime(solution, mapf) for (solution, mapf) in zip(solutions, instances)) /
    nb_instances

avg_flowtime_pred =
    sum(flowtime(solution, mapf) for (solution, mapf) in zip(solutions_pred, instances)) /
    nb_instances

avg_flowtime_naive =
    sum(flowtime(solution, mapf) for (solution, mapf) in zip(solutions_naive, instances)) /
    nb_instances


text:
"""
function logmean(aL,aR)
"""

function logmean(aL,aR)
    return logmean(aL,aR,log(aL),log(aR))
end

"""
function logmean(aL,aR,logL,logR)
    Compute logarithmic mean using pre-computed log values
"""

function logmean(aL,aR,logL,logR)

    "from: Entropy stable num. approx. for the isothermal and polytropic Euler"
    # assumes the specific value γ = 1.4
    da = aR-aL;
    aavg = .5*(aR+aL);
    f = da/aavg;
    v = f^2;
    if abs(f)<1e-4
        return aavg*(1 + v*(-.2-v*(.0512 - v*0.026038857142857)))
    else
        return -da/(logL-logR)
    end
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using GRAYTR
using Test

# Geometry package tests
include("test_geometry.jl")

# Utility function tests
include("test_utility.jl")

# Tests of rays.jl
include("test_rays.jl")

# Tests of bounding.jl
include("test_bounding.jl")

# Tests of shapes
include("test_shapes.jl")

# Tests of spectra
include("test_spectrum.jl")

# Tests of brdf.jl
include("test_brdf.jl")

# Tests of geometric primitives
include("test_primitives.jl")

# Tests of dummy aggregate
include("test_dumbaggregate.jl")

# Tests of samplers.jl
include("test_samplers.jl")

# Tests of filters.jl
include("test_filters.jl")

# Tests of cameras.jl
include("test_cameras.jl")

# Full integration test for image rendering.
include("test_full_computation.jl")

# Full radiation pressure test computation
include("test_pressure.jl")

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
function node_depth(args...;ntip::Int=Int[], nnode::Int=Int[],
    e1::Array{Int64,1}=Int[],e2::Array{Int64,1}=Int[], nedge::Int=Int[],method::Int=Int[],xx::Array{Int64,1}=Int[])

    for i in 1:ntip
        xx[i] = 1
    end

    if method == 1
        for i in 1:nedge
            xx[e1[i]] = xx[e1[i]] + xx[e2[i]]
        end
    elseif method == 2
        for i in 1:nedge
            if (xx[e1[i]] != 0)
                if xx[e1[i]] >=xx[e1[i]] + 1
                    continue
                end
                xx[ei[i]] = xx[e2[i]] + 1
            end
        end
    end
    return xx
end


text:
lines = readlines("inputs/05.txt")
function day_5(lines)
  grid = zeros(Int64, 1000,1000)
  pairs = []

  for line in lines
    x1 = parse(Int64, split(line,',')[1])
    y2 = parse(Int64, split(line,',')[3])
    mid = split(line,',')[2]
    y1 = parse(Int64, split(mid, ' ')[1])
    x2 = parse(Int64, split(mid, ' ')[3])
    push!(pairs, [x1, y1, x2, y2])
  end

  for p in pairs
    x1, y1, x2, y2 = p
    if x1 == x2
      for y in min(y1, y2):max(y1, y2)
        grid[x1, y] += 1
      end
    elseif y1 == y2
      for x in min(x1, x2):max(x1, x2)
        grid[x, y1] += 1
      end
    end
  end
  count(i -> i > 1, grid)
end

println(day_5(lines))


text:
# AUTO GENERATED FILE - DO NOT EDIT

export drawer

"""
    drawer(;kwargs...)
    drawer(children::Any;kwargs...)
    drawer(children_maker::Function;kwargs...)

A Drawer component.
Display overlay area at any side of the screen. For more information, see: https://mantine.dev/core/drawer/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Drawer children components
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `hideCloseButton` (Bool; optional): Hides close button, modal still can be closed with escape key and by clicking outside
- `noCloseOnClickOutside` (Bool; optional): Disable onClock trigger for outside events
- `noCloseOnEscape` (Bool; optional): Disable onClock trigger for escape key press
- `noFocusTrap` (Bool; optional): Disables focus trap
- `noOverlay` (Bool; optional): Removes overlay entirely
- `noScrollLock` (Bool; optional): Disables scroll lock
- `opened` (Bool; optional): If true drawer is mounted to the dom
- `overlayColor` (String; optional): Sets overlay color, defaults to theme.black in light theme and to theme.colors.dark[9] in dark theme
- `overlayOpacity` (Real; optional): Sets overlay opacity, defaults to 0.75 in light theme and to 0.85 in dark theme
- `padding` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Drawer body padding from theme or number for padding in px
- `position` (a value equal to: "right", "center", "left", "apart"; optional): Drawer body position
- `size` (a value equal to: "right", "center", "left", "apart" | String | Real; optional): Drawer body width (right | left position) or height (top | bottom position), cannot exceed 100vh for height and 100% for width
- `title` (String; optional): Drawer title, displayed in header before close button
- `zIndex` (Real; optional): Popper zIndex
"""
function drawer(; kwargs...)
        available_props = Symbol[:children, :id, :className, :hideCloseButton, :noCloseOnClickOutside, :noCloseOnEscape, :noFocusTrap, :noOverlay, :noScrollLock, :opened, :overlayColor, :overlayOpacity, :padding, :position, :size, :title, :zIndex]
        wild_props = Symbol[]
        return Component("drawer", "Drawer", "dash_mantine_components", available_props, wild_props; kwargs...)
end

drawer(children::Any; kwargs...) = drawer(;kwargs..., children = children)
drawer(children_maker::Function; kwargs...) = drawer(children_maker(); kwargs...)


text:
using Parameters
import Base: length

using SparseArrays

abstract type AbstractPhase end
export AbstractPhase

abstract type IdealPhase <: AbstractPhase end
export IdealPhase

struct EmptyPhase <: AbstractPhase end
export EmptyPhase

include("Calculators/Ratevec.jl")
include("Calculators/Thermovec.jl")
include("Reaction.jl")

@with_kw struct IdealGas{W<:Tuple,W2} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    spcdict::Dict{String,Int64}
    stoichmatrix::W2
    Nrp::Array{Float64,1}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::NASAvec
    otherreactions::Array{ElementaryReaction,1}
    diffusionlimited::Bool = false
end
IdealGas(species,reactions; name="",diffusionlimited=false) = IdealGas(species=species,reactions=reactions,name=name,
diffusionlimited=diffusionlimited,spcdict=Dict([sp.name=>sp.index for sp in species]))

function IdealGas(species,reactions; name="",diffusionlimited=false)
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    M,Nrp = getstoichmatrix(species,rxns)
    return IdealGas(species=species,reactions=rxns,name=name,
        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,veckinetics=vectuple, veckineticsinds=posinds, vecthermo=therm, otherreactions=otherrxns,diffusionlimited=diffusionlimited,)
end
export IdealGas

@with_kw struct IdealDiluteSolution{W<:Tuple,W2} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    solvent::Solvent
    stoichmatrix::W2
    Nrp::Array{Float64,1}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::NASAvec
    otherreactions::Array{ElementaryReaction,1}
    spcdict::Dict{String,Int64}
    diffusionlimited::Bool = true
end
IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited=true) = IdealDiluteSolution(species=species,reactions=reactions,
solvent=solvent,name=name,diffusionlimited=diffusionlimited,spcdict=Dict([sp.name=>sp.index for sp in species]))
function IdealDiluteSolution(species,reactions,solvent; name="",diffusionlimited=true)
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    M,Nrp = getstoichmatrix(species,rxns)
    return IdealDiluteSolution(species=species,reactions=rxns,solvent=solvent,name=name,
        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,veckinetics=vectuple,veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,diffusionlimited=diffusionlimited)
end
export IdealDiluteSolution

@with_kw struct IdealSurface{W<:Tuple,W2} <: IdealPhase
    name::String = ""
    species::Array{Species,1}
    reactions::Array{ElementaryReaction,1}
    sitedensity::Float64
    stoichmatrix::W2
    Nrp::Array{Float64,1}
    veckinetics::W
    veckineticsinds::Array{Int64,1}
    vecthermo::NASAvec
    otherreactions::Array{ElementaryReaction,1}
    spcdict::Dict{String,Int64}
    diffusionlimited::Bool = false
end
function IdealSurface(species,reactions,sitedensity; name="",diffusionlimited=false)
    @assert diffusionlimited==false "diffusionlimited=true not supported for IdealSurface"
    vectuple,vecinds,otherrxns,otherrxninds,posinds = getveckinetics(reactions)
    rxns = vcat(reactions[vecinds],reactions[otherrxninds])
    rxns = [ElementaryReaction(index=i,reactants=rxn.reactants,reactantinds=rxn.reactantinds,products=rxn.products,
        productinds=rxn.productinds,kinetics=rxn.kinetics,radicalchange=rxn.radicalchange,pairs=rxn.pairs) for (i,rxn) in enumerate(rxns)]
    therm = getvecthermo(species)
    M,Nrp = getstoichmatrix(species,rxns)
    return IdealSurface(species=species,reactions=rxns,sitedensity=sitedensity,name=name,
        spcdict=Dict([sp.name=>sp.index for sp in species]),stoichmatrix=M,Nrp=Nrp,veckinetics=vectuple,veckineticsinds=posinds,vecthermo=therm,otherreactions=otherrxns,diffusionlimited=diffusionlimited)
end
export IdealSurface

"""
construct the stochiometric matrix for the reactions and the reaction molecule # change
"""
function getstoichmatrix(spcs,rxns)
    M = spzeros(length(rxns),length(spcs))
    Nrp = zeros(length(rxns))
    for (i,rxn) in enumerate(rxns)
        Nrp[i] = Float64(length(rxn.productinds) - length(rxn.reactantinds))
        for ind in rxn.reactantinds
            M[i,ind] += 1.0
        end
        for ind in rxn.productinds
            M[i,ind] -= 1.0
        end
    end
    return M,Nrp
end

"""
Split the reactions into groups with the same kinetics
"""
function splitreactionsbykinetics(rxns)
    tps = []
    rxnlists = []
    for (i,rxn) in enumerate(rxns)
        typ = getkineticstype(rxn.kinetics)
        tind = findfirst(x->x==typ,tps)
        if tind == nothing
            push!(rxnlists,[])
            push!(tps,typ)
            push!(rxnlists[end],rxn)
        else
            push!(rxnlists[tind],rxn)
        end
    end
    return (tps,rxnlists)
end
export splitreactionsbykinetics

"""
create vectorized kinetics calculators for the reactions
"""
function getveckinetics(rxns)
    tps,rxnlists = splitreactionsbykinetics(rxns)
    posinds = Array{Int64,1}()
    fs = []
    otherrxns = Array{ElementaryReaction,1}()
    vecinds = Array{Int64,1}()
    otherrxninds = Array{Int64,1}()
    for (i,tp) in enumerate(tps)
        if typeof(tp)<:Tuple
            typ = split(tp[1],".")[end] #this split needs done because in pyrms the type names come as ReactionMechanismSimulator.Arrhenius instead of Arrhenius in RMS proper
        else
            typ = split(tp,".")[end]
        end
        rinds = [findfirst(isequal(rxn),rxns) for rxn in rxnlists[i]]
        fcn = Symbol(typ * "vec")
        if !(fcn in allowedfcnlist) || occursin("Troe",typ) #no vectorized kinetics or for the moment stuff with efficiencies
            append!(otherrxns,rxnlists[i])
            append!(otherrxninds,rinds)
        else
            append!(vecinds,rinds)
            fcn = eval(fcn)
            x = fcn([x.kinetics for x in rxnlists[i]])
            push!(fs,x)
            if posinds == Array{Int64,1}()
                push!(posinds,length(rinds))
            else
                push!(posinds,length(rinds)+posinds[end])
            end
        end
    end
    vectuple = tuple(fs...)
    return (vectuple,vecinds,otherrxns,otherrxninds,posinds)
end
export getveckinetics

"""
create vectorized thermo calculator for species
"""
function getvecthermo(spcs)
    thermo = [sp.thermo for sp in spcs]
    typeassert.(thermo,NASA)
    return NASAvec([sp.thermo for sp in spcs])
end
export getvecthermo

length(p::T) where {T<:AbstractPhase} = 1
export length

iterate(p::T) where {T<:AbstractPhase} = p
export iterate

Broadcast.broadcastable(p::T) where {T<:AbstractPhase} = Ref(p)
export broadcastable


text:
@testset "Geometries" begin
  @testset "Rectangle" begin
    r = Rectangle((1.,1.), (2.,3.))
    @test extrema(r) == ([1.,1.], [2.,3.])
    @test sides(r) == [1.,2.]
    @test GeoStatsBase.center(r) == [3/2, 4/2]
    @test GeoStatsBase.diagonal(r) == sqrt(1^2 + 2^2)
    @test volume(r) == 1*2
  end
end