text:
using Test, ParticleDA, MPI, Statistics, Random

MPI.Init()
my_rank = MPI.Comm_rank(MPI.COMM_WORLD)
my_size = MPI.Comm_size(MPI.COMM_WORLD)

Random.seed!(1234 + my_rank)

nprt_per_rank = 5
n = nprt_per_rank * my_size

local_particles = rand(nprt_per_rank)

for i = 1:my_size
    if i == my_rank + 1
        println("rank ",my_rank,": particle states: ", local_particles)
    end
    MPI.Barrier(MPI.COMM_WORLD)
end

@testset for size in ((1,1,1), (1,1,1,1))

    stats = Array{ParticleDA.SummaryStat{Float64},length(size)}(undef, size)
    ParticleDA.get_mean_and_var!(stats,reshape(local_particles,size...,nprt_per_rank),0)
    @test_throws ErrorException ParticleDA.get_mean_and_var!(stats,reshape(local_particles,size[begin:end-1]...,nprt_per_rank),0)

    global_particles = MPI.Gather(local_particles, 0, MPI.COMM_WORLD)

    if my_rank == 0
        gather_mean = Statistics.mean(global_particles)
        gather_var = Statistics.var(global_particles, corrected=true)

        println("Mean     ",stats[1].avg, " -- ", stats[1].avg ≈ gather_mean)
        println("Variance ",stats[1].var, " -- ", stats[1].var ≈ gather_var)

        @test stats[1].avg ≈ gather_mean
        @test stats[1].var ≈ gather_var
    end

end


text:
import DataStructures: SortedDict

export priceEuropean, AdaptiveFilonCharFuncPricer

#Adaptive Filon with variable transform.
struct AdaptiveFilonCharFuncPricer{T}
    τ::T
    kcos::Array{T,2}
    vControl::T
    pi::T
    function AdaptiveFilonCharFuncPricer(
        p::CharFunc,
        τ::T;
        qTol::T = sqrt(eps(T)),
    ) where {T}
        mcos = Dict{T,Tuple{T,T}}()
        iPure = oneim(p)
        @inline function integrand(z::T)::Tuple{T,T} where {T}
            if z == one(T)
                #x=+infty
                return zero(T), zero(T)
            end
            x = (1 + z) / (1 - z)
            u = x - iPure / 2
            phi = evaluateCharFunc(p, u, τ)
            denom = (1 / 4 + x^2)
            phi /= denom
            factor = 2 / (1 - z)^2
            phi *= factor
            return real(phi), imag(phi)
        end

        @inline function fcos(x::T)::T where {T}
            v = integrand(x)
            if (x != 1)
                mcos[x] = v
            end
            return v[1]
        end
        @inline function fsin(x::T)::T where {T}
            # v = get!(mcos,x , integrand(x))
            local v
            if haskey(mcos, x)
                v = mcos[x]
            else
                v = integrand(x)
                if (x != 1)
                    mcos[x] = v
                end
            end
            return v[2]
        end
        integrateSimpsonGG(fcos, -one(T), one(T), qTol, 16)
        integrateSimpsonGG(fsin, -one(T), one(T), qTol, 16)
        sortedKeys = sort!(collect(keys(mcos)))
        kcos = zeros(T, (3, length(mcos) * 2 - 1))
        #convert to original coordinate
        for (i, u) in enumerate(sortedKeys)
            v = mcos[u]
            x = (1 + u) / (1 - u)
            factor = 2 / (1 - u)^2
            kcos[:, i*2-1] = [x, v[1] / factor, v[2] / factor]
        end
        #add mid points
        for i = 1:length(mcos)-1
            a = @view kcos[:, i*2-1]
            b = @view kcos[:, (i+1)*2-1]
            x = (a[1] + b[1]) / 2
            u = (x - 1) / (x + 1) #(1-u)*x = 1+u => u(1+x) = x-1 => u=(x-1)/(x+1)
            v = integrand(u)
            factor = 2 / (1 - u)^2
            kcos[:, i*2] = [x, v[1] / factor, v[2] / factor]
        end
        #end interval = kcos[1,end]
        return new{T}(τ, kcos, getControlVariance(p, τ), const_pi(p))
    end

end

function priceEuropean(
    p::AdaptiveFilonCharFuncPricer{T},
    isCall::Bool,
    strike::T,
    forward::T,
    τ::T,
    discountDf::T,
) where {T}
    if τ != p.τ
        throw(DomainError(τ, string("maturity is different from pricer maturity ", p.τ)))
    end
    #intervals are not equals anymore! => cubic hermite spline f, fp.
    #for now, double the points, note: need to transform f and f' back as well.
    freq = log(strike / forward) #-k
    sqrtEps = sqrt(eps(T))
    integral = zero(T)
    diff = zero(T)
    diff1 = zero(T)
    kcos = p.kcos
    nx = size(kcos)[2] #cols

    @inbounds for i = 1:2:nx-2
        fp = FilonParams(freq, kcos[1, i], kcos[1, i+2])
        integral2 = filonCosSin3Points(
            fp,
            kcos[1, i],
            kcos[1, i+1],
            kcos[1, i+2],
            kcos[2, i],
            kcos[2, i+1],
            kcos[2, i+2],
            kcos[3, i],
            kcos[3, i+1],
            kcos[3, i+2],
            freq,
        )
        integral += integral2
    end
    callPrice = -discountDf * sqrt(strike * forward) / p.pi * integral
    if p.vControl != 0
        #use control variate
        variance = p.vControl * p.τ
        sqrtVar = sqrt(variance)
        d1 = log(forward / strike) / sqrtVar + sqrtVar / 2
        d2 = d1 - sqrtVar
        nd1 = normcdf(d1)
        nd2 = normcdf(d2)
        price = discountDf * (forward * nd1 - strike * nd2)
        callPrice += price
    else
        callPrice += discountDf * forward
    end
    if isCall
        return callPrice
    end
    putPrice = callPrice - (forward - strike) * discountDf
    return putPrice
end

struct FilonParams{T}
    α::T
    β::T
    γ::T
    cost::T
    sint::T
end

function FilonParams(t::T, x0::T, x2::T) where {T}
    theta = t * (x2 - x0) / 2
    return computeFilonParams(theta)
end

function belowFilonThreshold(theta::T)::Bool where {T}
    return abs(theta) <= 15 * eps(T)^0.125
end

function belowFilonThreshold(theta::Float64)::Bool
    return theta <= 1.0 / 6.0
end

@inline function computeFilonParams(theta::T) where {T}
    theta2 = theta^2
    local α, β, γ
    sint, cost = sincos(theta)

    if belowFilonThreshold(theta)
        α = theta * (theta2 * (2 / 45 + theta2 * (-2 / 315 + 2 * theta2 / 4725)))
        β =
            2 / 3 +
            theta2 *
            (2 / 15 + (theta2 * (-4 / 105 + theta2 * (2 / 567 - 4 * theta2 / 22275))))
        γ = 4 / 3 + theta2 * (-2 / 15 + theta2 * (1 / 210 - theta2 / 11340))
    else
        theta3 = theta2 * theta
        α = (theta2 + theta * sint * cost - 2.0 * sint * sint) / theta3
        β = (2 * theta + 2 * theta * cost * cost - 4 * sint * cost) / theta3
        γ = 4 * (sint - theta * cost) / theta3
    end
    return FilonParams(α, β, γ, cost, sint)
end

@inline function filonCosSin3Points(
    p::FilonParams{T},
    x0::T,
    x1::T,
    x2::T,
    fcos0::T,
    fcos1::T,
    fcos2::T,
    fsin0::T,
    fsin1::T,
    fsin2::T,
    t::T,
)::T where {T}
    if x0 == x2
        return T(0)
    end

    h = (x2 - x0) / 2

    st0, ct0 = sincos(t * x0)
    st1 = p.sint * ct0 + p.cost * st0
    ct1 = p.cost * ct0 - p.sint * st0
    st2 = p.sint * ct1 + p.cost * st1
    ct2 = p.cost * ct1 - p.sint * st1

    s2n = fsin0 * st0 / 2
    s2n += fsin2 * st2 / 2

    s2nm1 = fsin1 * st1

    value = h * (p.α * (fsin0 * ct0 - fsin2 * ct2) + p.β * s2n + p.γ * s2nm1)
    c2n = fcos0 * ct0 / 2
    c2n += fcos2 * ct2 / 2

    c2nm1 = fcos1 * ct1
    value += h * (p.α * (fcos2 * st2 - fcos0 * st0) + p.β * c2n + p.γ * c2nm1)
    return value
end

function integrateSimpsonGG(f, a::T, b::T, tol::T, maxRecursionDepth::Int)::T where {T}
    if tol < eps(T)
        tol = eps(T)
    end

    if b <= a
        return Base.zero(T)
    end
    c = (a + b) / 2

    fa = f(a)
    fm = f(c)
    fb = f(b)
    yy0 = f(a + 0.9501 * (b - a))
    yy1 = f(a + 0.2311 * (b - a))
    yy2 = f(a + 0.6068 * (b - a))
    yy3 = f(a + 0.4860 * (b - a))
    yy4 = f(a + 0.8913 * (b - a))

    is = (b - a) / 8 * (fa + fm + fb + yy0 + yy1 + yy2 + yy3 + yy4)
    if is == 0
        is = b - a
    end
    is = is * tol / eps(T)

    result = adaptiveSimpsonAux(f, a, b, fa, fm, fb, is, maxRecursionDepth)
    return result
end

function adaptiveSimpsonAux(
    f,
    a::T,
    b::T,
    fa::T,
    fm::T,
    fb::T,
    is::T,
    bottom::Int,
)::T where {T}
    m = (b + a) / 2
    h = (b - a) / 4
    fml = f(a + h)
    fmr = f(b - h)
    i1 = h / 1.5 * (fa + 4 * fm + fb)
    i2 = h / 3 * (fa + 4 * (fml + fmr) + 2 * fm + fb)
    #i1 = (16 * i2 - i1) / 15
    if ((is + (i1 - i2)) == is) || (m <= a) || (b <= m)
        if ((m <= a) || (b <= m))
            #Interval contains no more machine numbers, required tolerance may not be met
        end
        return i2
    end

    if bottom <= 0
        # println("Number of maximum recursions reached")
        return i2
    end

    return adaptiveSimpsonAux(f, a, m, fa, fml, fm, is, bottom - 1) +
           adaptiveSimpsonAux(f, m, b, fm, fmr, fb, is, bottom - 1)
end

ext: jl


text:
export Bam, BamLoader

abstract AbstractBam

immutable BamLoader <: AbstractBam
    file::AbstractString
    header_chunk::Bytes
    refs::Vector{Tuple{String, i32}}
    handle::IO

    function BamLoader(file::AbstractString)
        f = open(file) |> ZlibInflateInputStream
        @assert f >> 4 == [0x42, 0x41, 0x4d, 0x01]

        l_text = f >> i32
        text   = f >> l_text
        n_ref  = f >> i32

        refs   = Vector{Tuple{String, i32}}(n_ref)
        for i in 1:n_ref
            l_name  = f >> i32
            name    = f >> (l_name - 1)
            l_ref   = f >>> 1 >> i32
            refs[i] = name, l_ref
        end

        new(file, text, refs, f)
    end
end

start(bl::BamLoader)            = nothing
next(bl::BamLoader, ::Void)     = Read(bl.handle), nothing
done(bl::BamLoader, ::Void)     = eof(bl.handle)
iteratorsize(::Type{BamLoader}) = Base.SizeUnknown()
eltype(::Type{BamLoader})       = Read

show(io::IO, bl::BamLoader) = show(io, "BamLoader($(bl.file))")

immutable Bam <: AbstractBam
    file::AbstractString
    header_chunk::Bytes
    refs::Vector{Tuple{String, i32}}
    reads::Vector{Read}
end

Bam(file::AbstractString) = let bl = BamLoader(file)
    Bam(bl.file, bl.header_chunk, bl.refs, collect(bl))
end

start(bam::Bam)           = start(bam.reads)
next(bam::Bam, x)         = next(bam.reads, x)
done(bam::Bam, x)         = done(bam.reads, x)
iteratorsize(::Type{Bam}) = iteratorsize(Vector{Read})
eltype(::Type{Bam})       = eltype(Vector{Read})

show(io::IO, bam::Bam) = show(io, "Bam($(bam.file))")


text:
"""
```
get_data(f::Fred, series::AbstractString; kwargs...)
```

Request one series using the FRED API.

### Arguments
- `f`: Fred connection object
- `series`: series mnemonic

### Optional Arguments
`kwargs...`: key-value pairs to be appended to the FRED request. Accepted keys include:

- `realtime_start`: the start of the real-time period as YYYY-MM-DD string
- `realtime_end`: the end of the real-time period as YYYY-MM-DD string
- `limit`: maximum number of results to return
- `offset`: non-negative integer
- `sort_order`: `"asc"`, `"desc"`
- `observation_start`: the start of the observation period as YYYY-MM-DD string
- `observation_end`: the end of the observation period as YYYY-MM-DD string
- `units`: one of `"lin"`, `"chg"`, `"ch1"`, `"pch"`, `"pc1"`, `"pca"`, `"cch"`, `"cca"`,
  `"log"`
- `frequency`: one of `"d"`, `"w"`, `"bw"`, `"m"`, `"q"`, `"sa"`, `"a"`, `"wef"`,
  `"weth"`, `"wew"`, `"wetu"`, `"wem"`, `"wesu"`, `"wesa"`, `"bwew"`, `"bwem"`
- `aggregation_method`: one of `"avg"`, `"sum"`, `"eop"`
- `output_type`: one of `1` (obsevations by real-time period), `2` (observations by vintage
  date, all observations), `3` (observations by vintage date, new and revised observations
  only), `4` (observations, initial release only)
- `vintage_dates`: vintage dates as comma-separated YYYY-MM-DD strings
"""
function get_data(f::Fred, series::AbstractString; kwargs...)
    # Validation
    validate_args!(kwargs)

    # Setup
    metadata_url = get_api_url(f) * "series"
    obs_url      = get_api_url(f) * "series/observations"
    api_key      = get_api_key(f)

    # Add query parameters
    metadata_params = Dict("api_key"   => api_key,
                           "file_type" => "json",
                           "series_id" => series)
    obs_params = copy(metadata_params)

    # Query observations. Expand query dict with kwargs. Do this first so we can use the
    # calculated realtime values for the metadata request.
    for (key, value) in kwargs
        obs_params[string(key)] = string(value)
    end
    obs_response = HTTP.request("GET", obs_url, []; query=obs_params)
    obs_json = JSON.parse(String(copy(obs_response.body)))

    # Parse observations
    realtime_start  = obs_json["realtime_start"]
    realtime_end    = obs_json["realtime_end"]
    transformation_short = obs_json["units"]

    df = parse_observations(obs_json["observations"])

    # Query metadata
    metadata_params["realtime_start"] = realtime_start
    metadata_params["realtime_end"] = realtime_end
    metadata_response = HTTP.request("GET", metadata_url, []; query=metadata_params)
    metadata_json = JSON.parse(String(copy(metadata_response.body)))
    # TODO catch StatusError and just return incomplete data to the caller

    # Parse metadata
    metadata_parsed = Dict{Symbol, AbstractString}()
    for k in ["id", "title", "units_short", "units", "seasonal_adjustment_short",
        "seasonal_adjustment", "frequency_short", "frequency", "notes"]
        try
            metadata_parsed[Symbol(k)] = metadata_json["seriess"][1][k]
        catch err
            metadata_parsed[Symbol(k)] = ""
            @warn("Metadata '$k' not returned from server.")
        end
    end

    function parse_last_updated(last_updated)
        return DateTime(ZonedDateTime(last_updated, FRED_DATE_FORMAT), OUTPUT_TZ_TYPE)
    end
    last_updated = parse_last_updated(
        metadata_json["seriess"][1]["last_updated"])

    # format notes field
    metadata_parsed[:notes] = strip(replace(replace(
        metadata_parsed[:notes], r"[\r\n]" => " "), r" +" => " "))

    return FredSeries(metadata_parsed[:id], metadata_parsed[:title],
                      metadata_parsed[:units_short], metadata_parsed[:units],
                      metadata_parsed[:seasonal_adjustment_short],
                      metadata_parsed[:seasonal_adjustment],
                      metadata_parsed[:frequency_short], metadata_parsed[:frequency],
                      realtime_start, realtime_end, last_updated, metadata_parsed[:notes],
                      transformation_short, df,
                      df) # deprecated
end

# obs is a vector, of which each element is a dict with four fields,
# - realtime_start
# - realtime_end
# - date
# - value
function parse_observations(obs::Vector)
    n_obs = length(obs)
    value = Vector{Float64}(undef, n_obs)
    date  = Vector{Date}(undef, n_obs)
    realtime_start = Vector{Date}(undef, n_obs)
    realtime_end = Vector{Date}(undef, n_obs)
    for (i, x) in enumerate(obs)
        try
            value[i] = parse(Float64, x["value"])
        catch err
            value[i] = NaN
        end
        date[i]           = Date(x["date"], "yyyy-mm-dd")
        realtime_start[i] = Date(x["realtime_start"], "yyyy-mm-dd")
        realtime_end[i]   = Date(x["realtime_end"], "yyyy-mm-dd")
    end
    return DataFrame(realtime_start=realtime_start, realtime_end=realtime_end,
                     date=date, value=value)
end

# Make sure everything is of the right format.
# kwargs is a vector of Tuple{Symbol, Any}.
isyyyymmdd(x) = occursin(r"^[0-9]{4}-[0-9]{2}-[0-9]{2}$", x)
function validate_args!(kwargs)
    d = Dict(kwargs)

    # dates
    for k in [:realtime_start, :realtime_end, :observation_start, :observation_end]
        if (v = pop!(d, k, nothing)) != nothing && !isyyyymmdd(v)
                error("$k: Invalid date format: $v")
        end
    end
    # limit and offset
    for k in [:limit, :offset]
        if (v = pop!(d, k, nothing)) != nothing &&
            ( !(typeof(v) <: Number ) || typeof(v) <: Number && !(v>0) )
                error("$k: Invalid format: $v")
        end
    end
    # units
    if (v = pop!(d, :units, nothing)) != nothing &&
        v ∉ ["lin", "chg", "ch1", "pch", "pc1", "pca", "cch", "log"]
            error("units: Invalid format: $v")
    end
    # frequency
    if (v = pop!(d, :frequency, nothing)) != nothing &&
        v ∉ ["d", "w", "bw", "m", "q", "sa", "a", "wef", "weth", "wew", "wetu", "wem",
             "wesu", "wesa", "bwew", "bwem"]
            error("frequency: Invalid format: $v")
    end
    # aggregation_method
    if (v = pop!(d, :aggregation_method, nothing)) != nothing &&
        v ∉ ["avg", "sum", "eop"]
            error("aggregation_method: Invalid format: $v")
    end
    # output_type
    if (v = pop!(d, :output_type, nothing)) != nothing &&
        v ∉ [1, 2, 3, 4]
            error("output_type: Invalid format: $v")
    end
    # vintage dates, and too early vintages
    if (v = pop!(d, :vintage_dates, nothing)) != nothing
        vds_arr = split(string(v), ",")
        vds_bad = map(x -> !isyyyymmdd(x), vds_arr)
        if any(vds_bad)
            error("vintage_dates: Invalid date format: $(vds_arr[vds_bad])")
        end
        vds_early = map(x -> x<EARLY_VINTAGE_DATE, vds_arr)
        if any(vds_early)
            @warn(:vintage_dates, ": Early vintage date, data might not exist: ",
                vds_arr[vds_early])
        end
    end
    # all remaining keys have unspecified behavior
    if length(d) > 0
        for k in keys(d)
            @warn(string(k), ": Bad key. Removed from query.")
            deleteat!(kwargs, findall(x -> x[1]==k, kwargs))
        end
    end
end


text:
Sys.iswindows() && error("The Shoco C library does not support Windows")

using BinaryProvider

const verbose = "--verbose" in ARGS
const prefix = Prefix(get(filter(!isequal("verbose"), ARGS), 1, joinpath(@__DIR__, "usr")))
products = [LibraryProduct(prefix, ["libshoco"], :libshoco)]

bin_prefix = "https://github.com/ararslan/ShocoBuilder/releases/download/v0.1.0"

downloads = Dict(
    FreeBSD(:x86_64) =>
        ("$bin_prefix/Shoco.v0.1.0.x86_64-unknown-freebsd11.1.tar.gz",
         "d6b192c4e965659cac273e499a586fa4f2094de141dac3dfd192e636f7a3c2ee"),
    Linux(:aarch64, :glibc) =>
        ("$bin_prefix/Shoco.v0.1.0.aarch64-linux-gnu.tar.gz",
         "d44588e4268f36aafa05ae9a1e2c40386d0289955c4fc192663c2d857ca2663c"),
    Linux(:aarch64, :musl) =>
        ("$bin_prefix/Shoco.v0.1.0.aarch64-linux-musl.tar.gz",
         "c21ca1c265e28adec6f5e4fd630892b1f64283987bbf3170ddddece3fbe2abf8"),
    Linux(:armv7l, :glibc, :eabihf) =>
        ("$bin_prefix/Shoco.v0.1.0.arm-linux-gnueabihf.tar.gz",
         "96a21441085150c4cfcb31e4a943e5fcf8055ddcca005c364a44efaa734d7937"),
    Linux(:armv7l, :musl, :eabihf) =>
        ("$bin_prefix/Shoco.v0.1.0.arm-linux-musleabihf.tar.gz",
         "7d55f6370dad1c7e38aff06675690913156daebe45c9eb53b913f6123ab70ad6"),
    Linux(:i686, :glibc) =>
        ("$bin_prefix/Shoco.v0.1.0.i686-linux-gnu.tar.gz",
         "0d8a07d3ae4070765284ef0b2c8d3a830c17875db3b6833550163e2f64c3cf6a"),
    Linux(:i686, :musl) =>
        ("$bin_prefix/Shoco.v0.1.0.i686-linux-musl.tar.gz",
         "b0b2af8b972a1a5d6afa50acc67c7a08d6b83c5ca745cc6408c10aab293ca878"),
    Linux(:powerpc64le, :glibc) =>
        ("$bin_prefix/Shoco.v0.1.0.powerpc64le-linux-gnu.tar.gz",
         "f55a47099cd38b0ad613f71990c94fd9beeb1e45704844b1615df4b910f7cfdc"),
    Linux(:x86_64, :glibc) =>
        ("$bin_prefix/Shoco.v0.1.0.x86_64-linux-gnu.tar.gz",
         "17de1ab56a576e458632c19aec5e0e3031f988730f7208944a9b13620f38902a"),
    Linux(:x86_64, :musl) =>
        ("$bin_prefix/Shoco.v0.1.0.x86_64-linux-musl.tar.gz",
         "e8207c469fa6e0ddcc80d9123bc642f5719e804c5ae1683264af9e184a4c0c23"),
    MacOS(:x86_64) =>
        ("$bin_prefix/Shoco.v0.1.0.x86_64-apple-darwin14.tar.gz",
         "0d6f1985638accdc55c95906e2f7c646f283dd7bb996448f5a90433b74f9f75f"),
)

unsatisfied = any(p->!satisfied(p, verbose=verbose), products)

if haskey(downloads, platform_key())
    url, hash = downloads[platform_key()]
    if unsatisfied || !isinstalled(url, hash, prefix=prefix)
        install(url, hash, prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    error("Your platform ($(triplet(platform_key()))) is not supported")
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
"""
    OPENINGS_URL::String
[URL](https://backend.econjobmarket.org/data/zz_public/json/Ads) for accessing the public data for job opening advertisements in JSON format from EconJobMarket.
"""
const OPENINGS_URL = "https://backend.econjobmarket.org/data/zz_public/json/Ads"
"""
    CATEGORIES_URL::String
[URL](https://backend.econjobmarket.org/data/zz_public/json/Categories) for accessing the public data for job opening categories in JSON format from EconJobMarket.
"""
const CATEGORIES_URL = "https://backend.econjobmarket.org/data/zz_public/json/Categories"
"""
    CATEGORIES::Vector{String}
Job opening advertisements categories.
```jldoctest
julia> EconJobMarket.CATEGORIES
34-element Array{String,1}:
 "Development; Growth"
 "Econometrics"
 "Economic History"
 "Environmental; Ag. Econ."
 "Experimental Economics"
 "Finance"
 "Industrial Organization"
 "International Finance/Macro"
 "International Trade"
 "Labor; Demographic Economics"
 "Law and Economics"
 "Macroeconomics; Monetary"
 "Microeconomics"
 "Public Economics"
 "Theory"
 "Behavioral Economics"
 ⋮
 "Urban; Rural; Regional Economics"
 "Health; Education; Welfare"
 "Business Economics"
 "Computational Economics"
 "Political Economy"
 "Accounting"
 "Decision Sciences"
 "Insurance"
 "Management, General"
 "Management, Health Care"
 "Management, Information Technology"
 "Marketing"
 "Operations Research"
 "Organizational Behavior"
 "Real Estate"
 "Statistics"
```
"""
const CATEGORIES = JSON3.read(read(joinpath(dirname(@__DIR__), "data", "Categories.json"))) |>
    (obj -> getproperty.(obj, :name))
"""
    categories_indicator(obj::AbstractString)::BitVector
Given a string, it provides indicators for whether the job ad listed each opening category.

# Example
```jldoctest
julia> EconJobMarket.categories_indicator("2,3")
34-element BitArray{1}:
 0
 1
 1
 0
 0
 0
 0
 0
 0
 0
 ⋮
 0
 0
 0
 0
 0
 0
 0
 0
 0
```
"""
categories_indicator(obj::AbstractString) =
    1:34 .∈ Ref([ parse(Int, elem.match) for elem ∈ eachmatch(r"\d+", obj) ])
"""
    parse_opening(node)::NamedTuple
For each opening, it parses the information.
"""
function parse_opening(node)
    posid = isnothing(node.posid) ? missing : parse(Int, node.posid)
    oid = isnothing(node.oid) ? missing : parse(Int, node.oid)
    adtitle = isnothing(node.adtitle) ? missing : node.adtitle
    position_type_id = isnothing(node.position_type) ? missing : node.position_type_id
    position_type = isnothing(node.position_type) ? missing : node.position_type
    adtext =
        reduce((x, y) -> "$x $y",
               nodeText(node) for node ∈ eachmatch(Selector("body"),
                                                   parsehtml(node.adtitle).root))
    startdate = isnothing(node.startdate) ? missing : Date(node.startdate)
    enddate = isnothing(node.enddate) ? missing : Date(node.enddate)
    country_code = isnothing(node.country_code) ? missing : node.country_code
    position_country = isnothing(node.position_country) ? missing : node.position_country
    department = isnothing(node.department) ? missing : node.department
    shortname = isnothing(node.shortname) ? missing : node.shortname
    address = isnothing(node.address) ? missing : node.address
    longitude = isnothing(node.longitude) ? missing : parse(Float64, node.longitude)
    latitude = isnothing(node.latitude) ? missing : parse(Float64, node.latitude)
    name = isnothing(node.name) ? missing : node.shortname
    cats = isnothing(node.categories) ? missing : categories_indicator(node.categories)
    (; zip(Tuple(append!([ Symbol(clean_name(elem)) for elem ∈ propertynames(node)
                           if clean_name(elem) ≠ :categories ],
                           clean_name.(CATEGORIES))),
           (posid, oid, adtitle, position_type_id, position_type,
            adtext, startdate, enddate, country_code, position_country,
            department, shortname, address, longitude, latitude, name,
            cats...))...)
end
"""
    fetch_ads(filepath::AbstractString; delim::Union{Char,AbstractString} = \'\\t\')
Writes a text file with the latest job opening advertisements from EJM.

# Example
```jldoctest
julia> data = joinpath(dirname(@__DIR__), "data", "ads.tsv") |>
       fetch_ads |>
       CSV.read;

julia> Tables.schema(data)
Tables.Schema:
 :posid                              Int64
 :oid                                Int64
 :adtitle                            String
 :position_type_id                   String
 :position_type                      String
 :adtext                             String
 :startdate                          Date
 :enddate                            Date
 :country_code                       String
 :position_country                   String
 :department                         String
 :shortname                          String
 :address                            Missing
 :longitude                          Float64
 :latitude                           Float64
 ⋮
 :business_economics                 Bool
 :computational_economics            Bool
 :political_economy                  Bool
 :accounting                         Bool
 :decision_sciences                  Bool
 :insurance                          Bool
 :management_general                 Bool
 :management_health_care             Bool
 :management_information_technology  Bool
 :marketing                          Bool
 :operations_research                Bool
 :organizational_behavior            Bool
 :real_estate                        Bool
 :statistics                         Bool
```
"""
fetch_ads(filepath::AbstractString; delim::Union{Char,AbstractString} = '\t') =
    OPENINGS_URL |>
    HTTP.get |>
    (response -> response.body) |>
    String |>
    JSON3.read |>
    (data -> CSV.write(filepath,
                       columntable(parse_opening(node) for node ∈ data),
                       delim = delim))


text:
using PointPatternStatistics
using Test
include("reference_implementation.jl")
@testset "PointPatternStatistics.jl" begin
    window = (x=(-1.3,1), y=(0,1.6))
    randu((a, b)) = a + rand() * (b-a)
    randw(window) = (randu(window.x), randu(window.y))
    xy = [randw(window) for _ in 1:300]
    r = range(0, 0.25, length=512)
    @testset "iterate" begin
        @test collect(PointPattern(xy, window)) == xy
    end

    @testset "pcf" begin
        est = pcf(PointPattern(xy, window), r)
        ref = pcfref(xy, window, r)
        @test sum(abs2, (est - ref)[2:end]) < 0.0001
    end

    @testset "K" begin
        est = Kest(PointPattern(xy, window), r)
        ref = Kref(xy, window, r)
        @test est ≈ ref
    end

    @testset "F" begin
        include("F.jl")
    end
    @testset "G" begin
        include("G.jl")
    end
    #include("envelope.jl")
    @testset "SSI" begin
        include("SSI.jl")
    end
end


text:
# sample function, its gradient & derivative definintion

foo(x) = sum(x, 2)
foo_grad(dzdy, x) = repmat(dzdy, 1, size(x, 2))
@specialdiff (y = foo(x)) (dz!dx = foo_grad(dz!dy, x))

function test_special(u, v)
    x = u .+ v
    y = foo(x)
    z = sum(1.0 .* y) # note: don't remove 1.0
end

test_compare(test_special; u=randn(3, 3), v=randn(3, 3))


text:
function from_vk_call(x::Spec, identifier = :x)
    prop = :($identifier.$(x.name))
    x.name == :pNext && return :(load_next_chain($prop, next_types...))
    x.name == :sType && return nothing
    jtype = hl_type(x)
    @match x begin
        # array pointer
        GuardBy(is_arr) => @match jtype begin
            :(Vector{$_}) => :(unsafe_wrap($jtype, $prop, $(len_expr(x, identifier)); own = true))
        end

        GuardBy(is_length) => nothing
        _ => from_vk_call(prop, x.type, jtype)
    end
end

function len_expr(x::Spec, identifier)
    @assert !isnothing(x.len)
    params = children(parent_spec(x))
    postwalk(x.len) do ex
        ex == :/ && return :÷
        ex isa Symbol && ex in params.name && return :($identifier.$ex)
        ex
    end
end

function from_vk_call(prop, t, jtype = hl_type(t))
    @match t begin
        :Cstring => :(unsafe_string($prop))
        GuardBy(in(spec_flags.name)) || GuardBy(in(spec_enums.name)) => prop
        GuardBy(in(getproperty.(filter(!isnothing, spec_flags.bitmask), :name))) => :($jtype(UInt32($prop)))
        GuardBy(in(spec_handles.name)) => :($(remove_vk_prefix(t))($prop))
        :(NTuple{$_,$T}) && if jtype ≠ :String end => broadcast_ex(from_vk_call(prop, ntuple_type(t)))
        if follow_constant(t) == jtype
        end => prop
        if is_hl(jtype) end => :($jtype($prop))
        _ => :(from_vk($jtype, $prop))
    end
end

function vk_call(x::Spec)
    var = wrap_identifier(x.name)
    jtype = idiomatic_julia_type(x)
    @match x begin
        ::SpecStructMember && if x.type == :VkStructureType && parent(x) ∈ keys(structure_types)
        end => :(structure_type($(parent(x))))
        ::SpecStructMember && if is_semantic_ptr(x.type)
        end => :(unsafe_convert($(x.type), $var))
        if is_fn_ptr(x.type)
        end => var
        GuardBy(is_size) && if x.requirement == POINTER_REQUIRED
        end => x.name # parameter converted to a Ref already
        GuardBy(is_opaque) => var
        GuardBy(is_length) => begin
            !x.autovalidity && @debug "Automatic validation was disabled for length parameter $x."
            @match x begin
                GuardBy(is_length_exception) || GuardBy(!is_inferable_length) => var
                _ => :(pointer_length($(wrap_identifier(first(arglen(x)))))) # Julia works with arrays, not pointers, so the length information can directly be retrieved from them
            end
        end
        GuardBy(is_pointer_start) => 0 # always set first* variables to 0, and the user should provide a (sub)array of the desired length
        if x.type ∈ spec_handles.name
        end => var # handled by unsafe_convert in ccall

        # constant pointer to a unique object
        if is_ptr(x.type) &&
           !is_arr(x) &&
           (x.is_constant || (func = func_by_name(x.func); func.type == FTYPE_QUERY && x ≠ last(children(func))))
        end => @match x begin
            if ptr_type(x.type) ∈ [spec_structs.name; spec_unions.name]
            end => var # handled by cconvert and unsafe_convert in ccall
            if x.requirement == OPTIONAL
            end => :($var == $(default(x)) ? $(default(x)) : Ref($var)) # allow optional pointers to be passed as C_NULL instead of a pointer to a 0-valued integer
            _ => :(Ref($var))
        end
        if x.type ∈ [spec_flags.name; spec_enums.name]
        end => var
        if x.type ∈ getproperty.(filter(!isnothing, spec_flags.bitmask), :name)
        end => :($(x.type)($var.val))
        if x.type ∈ extension_types
        end => var
        _ => @match jtype begin
            :String || :Bool || :(Vector{$et}) || if jtype == follow_constant(x.type)
            end => var # conversions are already defined
            if x.type in [spec_structs.name; spec_unions.name] && jtype == struct_name(x.type)
            end => :($var.vks)
            :(NTuple{Int($_N),$_T}) => var
            _ => :(to_vk($(x.type), $var)) # fall back to the to_vk function for conversion
        end
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
"""
    AbstractSystem

Abstract supertype for all system types.
"""
abstract type AbstractSystem end

"""
    statedim(s::AbstractSystem)

Returns the dimension of the state space of system `s`.
"""
statedim(::AbstractSystem) = nothing

"""
    stateset(s::AbstractSystem)

Returns the set of allowed states of system `s`.
"""
stateset(::AbstractSystem) = nothing

"""
    inputdim(s::AbstractSystem)

Returns the dimension of the input space of system `s`.
"""
inputdim(::AbstractSystem) = nothing

"""
    inputset(s::AbstractSystem)

Returns the set of allowed inputs of system `s`.
"""
inputset(::AbstractSystem) = nothing

"""
    noisedim(s::AbstractSystem)

Returns the dimension of the noise space of system `s`.
"""
noisedim(::AbstractSystem) = nothing

"""
    noiseset(s::AbstractSystem)

Returns the set of allowed noises of system `s`.
"""
noiseset(::AbstractSystem) = nothing

"""
    AbstractDiscreteSystem

Abstract supertype for all discrete system types.
"""
abstract type AbstractDiscreteSystem <: AbstractSystem end

"""
    AbstractContinuousSystem

Abstract supertype for all continuous system types.
"""
abstract type AbstractContinuousSystem <: AbstractSystem end

"""
    islinear(s::AbstractSystem)

Specifies if the dynamics of system `s` is specified by linear equations.

### Notes

We adopt the notion from [Section 2.7, 1]. For example, the system with inputs
``x' = f(t, x, u) = A x + B u`` is linear, since the function ``f(t, ⋅, ⋅)`` is
linear in ``(x, u)`` for each ``t ∈ \\mathbb{R}``. On the other hand,
``x' = f(t, x, u) = A x + B u + c`` is affine but not linear, since it is not
linear in ``(x, u)``.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`. So, if a system type allows an instance that
is not linear, it returns `false` by default. For example, polynomial systems
can be nonlinear; hence `islinear` returns `false`.

[1] Sontag, Eduardo D. *Mathematical control theory: deterministic finite dimensional
systems.* Vol. 6. Springer Science & Business Media, 2013.
"""
function islinear(::AbstractSystem) end

"""
    isaffine(s::AbstractSystem)

Specifies if the dynamics of system `s` is specified by affine equations.

### Notes

An affine system is the composition of a linear system and a translation.
See [`islinear(::AbstractSystem)`](@ref) for the notion of linear system adopted
in this library.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`. So, if a system type allows an instance that
is not affine, it returns `false` by default. For example, polynomial systems
can be nonlinear; hence `isaffine` is `false`.
"""
function isaffine(::AbstractSystem) end

"""
    ispolynomial(s::AbstractSystem)

Specifies if the dynamics of system `s` is specified by polynomial equations.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`. Hence, e.g. a `LinearContinuousSystem` is not
considered to be of polynomial type.
"""
function ispolynomial(::AbstractSystem) end

"""
    isblackbox(s::AbstractSystem)

Specifies if no specific structure is assumed for the dynamics of system `s`.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`.
"""
function isblackbox(::AbstractSystem) end

"""
    isnoisy(s::AbstractSystem)

Determines if the dynamics of system `s` contains a noise term `w`.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`.
"""
function isnoisy(::AbstractSystem) end

"""
    iscontrolled(s::AbstractSystem)

Determines if the dynamics of system `s` contains a control input `u`.

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`.
"""
function iscontrolled(::AbstractSystem) end

"""
    isconstrained(s::AbstractSystem)

Determines if the system `s` has constraints on the state, input and noise,
respectively (those that are available).

The result of this function only depends on the system type, not the value, and
can also be applied to `typeof(s)`.
"""
function isconstrained(::AbstractSystem) end

"""
    AbstractMap

Abstract supertype for all map types.
"""
abstract type AbstractMap end

"""
    outputdim(m::AbstractMap)

Returns the dimension of the output space of the map `m`.
"""
outputdim(::AbstractMap) = nothing

"""
    outputmap(s::SystemWithOutput)

Returns the output map of a system with output.
"""
outputmap(::AbstractSystem) = nothing

"""
    islinear(m::AbstractMap)

Specifies if the map `m` is linear or not.

### Notes

A map is *linear* if it preserves the operations of scalar multiplication and
vector addition.
"""
function islinear(::AbstractMap) end

"""
    isaffine(m::AbstractMap)

Specifies if the map `m` is affine or not.

### Notes

An affine map is the composition of a linear map and a translation.
See also [`islinear(::AbstractMap)`](@ref).
"""
function isaffine(::AbstractMap) end

"""
    apply(m::AbstractMap, args...)

Apply the rule specified by the map to the given arguments.
"""
function apply end

"""
    state_matrix(s::AbstractSystem)

Return the state matrix of an affine system.

### Notes

The state matrix is the matrix proportional to the state, e.g. the matrix `A`
in the linear continuous system ``x' = Ax``.
"""
state_matrix(::AbstractSystem) = nothing

"""
    input_matrix(s::AbstractSystem)

Return the input matrix of a system with linear input.

### Notes

The input matrix is the matrix proportional to the input, e.g. the matrix `B`
in the linear continuous system with input, ``x' = Ax + Bu``.
"""
input_matrix(::AbstractSystem) = nothing

"""
    noise_matrix(s::AbstractSystem)

Return the noise matrix of a system with linear noise.

### Notes

The noise matrix is the matrix proportional to the noise, e.g. the matrix `D`
in the linear system with noise, ``x' = Ax + Dw``.
"""
noise_matrix(::AbstractSystem) = nothing

"""
    affine_term(s::AbstractSystem)

Return the affine term in an affine system.

### Notes

The affine term is e.g. the vector ``c`` in the affine system ``x' = Ax + c``.
"""
affine_term(::AbstractSystem) = nothing


text:
using GeometryTypes

const AInt = Int32
const AFloat = Float32

include("/Users/maximilianbittens/.julia/v0.6/Unity/src/MeshReaders.jl")
include("/Users/maximilianbittens/.julia/v0.6/Unity/src/AnsysReader.jl")

using GeometryTypes
using ColorTypes
import Unity: UnityMesh, PyramidMesh, UnityVector3, UnityText
using Colors

function triangle_area(vertices::Vector{Point3f0}, faces::Vector{Face{3,AInt}})
	areas = Vector{AFloat}();
	for face in faces
		push!(areas,triangle_area(vertices[face]...))
	end
	return areas
end

function triangle_area(P1::Point, P2::Point, P3::Point)
	L = (
		sqrt(sum((P1-P2).^2)),
		sqrt(sum((P2-P3).^2)),
		sqrt(sum((P3-P1).^2))
		)
	s = sum(L)/2.0
	return sqrt( s*(s-L[1])*(s-L[2])*(s-L[3]) )
end

function centroids(vertices::Vector{Point3f0}, faces::Vector{Face{3,AInt}})
	#c = Vector{Point{3,Float64}}()
	c = Vector{Point{3,Float32}}()
	for face in faces
		vs = vertices[face]
		#push!(c,sum(convert(Tuple{Point{3,Float64},Point{3,Float64},Point{3,Float64}},vs))/3.)
		push!(c,sum(vs)/3.)
	end
	return c
end

#
function vertices_and_faces(mesh)
	vertices = Vector{Point3f0}()
	for i = 1:size(mesh.nodes,2)
		push!(vertices,Point3f0(mesh.nodes[:,i]))
	end
	faces = map(x->Face{3,AInt}(x[4:end]),filter(x->x[1:3]==[2,3,1],mesh.elements))
	return vertices,faces
end
#

function covariance_matrix_and_weighted_center(mesh)
	vertices, faces = vertices_and_faces(mesh)
	return covariance_matrix_and_weighted_center(vertices, faces)
end

function covariance_matrix_and_weighted_center(vertices::Vector{Point3f0}, faces::Vector{Face{3,AInt}})
	face_centroids = centroids(vertices,faces)
	info("face_centroids")
	display(face_centroids)
	tri_areas = triangle_area(vertices,faces)
	weighted_center = sum(face_centroids.*tri_areas)/sum(tri_areas)
	#info("weighted_center")
	#display(weighted_center)
	zero_centered_centroids = map((x,y)->(x-weighted_center)*sqrt(y),face_centroids,tri_areas)
	#info("zero_centeredweighted_center")
	#display(zero_centered_centroids)
	wcm = hcat(map(x->[x...],zero_centered_centroids)...)
	#info("wcm")
	#display(wcm)
	return (wcm*wcm')./(length(zero_centered_centroids)-1), weighted_center
end

#
function minarc(v1,coords=[Point3f0(1,0,0),Point3f0(0,1,0),Point3f0(0,0,1)])
	arc(v1,v2) = acos( dot(v1,v2)/(norm(v1)*norm(v2)))
	a = [ arc(v1,coords[1]),arc(v1,coords[2]),arc(v1,coords[3])]
	ind = a[1]<min(a[2],a[3]) ? 1 : (a[2]<min(a[1],a[3]) ? 2 : ( (a[3]<min(a[1],a[2]) ? 3 : error() )))
	return a[ind], ind
end

function pca_weighted(mesh)
	vertices, faces = vertices_and_faces(mesh)
	return pca_weighted(vertices,faces)
end
#

function pca_weighted(vertices::Vector{Point3f0}, faces::Vector{Face{3,AInt}})
	try
	cov_mat, wc = covariance_matrix_and_weighted_center(vertices,faces)
	#info("cov")
	#display(cov_mat)
	λ, ϕ = eigs(cov_mat,nev=2)
	#display(ϕ)
	#pca = align(Point3f0(ϕ[:,1]/norm(ϕ[:,1])), Point3f0(ϕ[:,2]/norm(ϕ[:,2])))
	pca = [Point3f0(ϕ[:,1]/norm(ϕ[:,1])),Point3f0(ϕ[:,2]/norm(ϕ[:,2])),Point3f0(cross(ϕ[:,1],ϕ[:,2])/norm(cross(ϕ[:,1],ϕ[:,2])))]
	return pca,wc
	catch
		display(faces)
		error()
	end
end

function mul!(res::Vector{Float32},mat::StaticArrays.SArray{Tuple{3,3},Float32,2,9},vec::Point3f0)
	fill!(res,0.0f0)
	for i = 1:3
		for j = 1:3
			res[i] += mat[i,j]*vec[j]
		end
	end
end

import StaticArrays

type OBB #Oriented Bounding Box
	pca::StaticArrays.SArray{Tuple{3,3},Float32,2,9}
	minb::Vector{Float32}
	maxb::Vector{Float32}
	OBB(pca,min,max) = new(pca,min,max)
	OBB(pca) = new(hcat(pca...),[Inf32,Inf32,Inf32],[-Inf32,-Inf32,-Inf32])
end

function Base.copy(obb::OBB)
	return OBB(copy(obb.pca),copy(obb.minb),copy(obb.maxb))
end

function oriented_bounding_box(vertices::Vector{Point3f0},pca::Vector{Point3f0})
	obb = OBB(pca)
	verts = map(x->inv(obb.pca)*x,vertices)
	obb.minb,obb.maxb = boundingbox(verts)
	return obb
end

function boundingbox(vertices::Vector{Point3f0})
	return Point3f0(minimum(map(x->x[1],vertices)),minimum(map(x->x[2],vertices)),minimum(map(x->x[3],vertices))),Point3f0(maximum(map(x->x[1],vertices)),maximum(map(x->x[2],vertices)),maximum(map(x->x[3],vertices)))
end

function boundingbox_vertices(minb,maxb)
	Δ = maxb-minb
	vertices = Point3f0[minb, minb+[Δ[1],0,0],minb+[0,Δ[2],0],minb+[0,0,Δ[3]],minb+[Δ[1],Δ[2],0],minb+[Δ[1],0,Δ[3]],minb+[0,Δ[2],Δ[3]],maxb]
	return vertices
end

function GeometryTypes.volume(obb::OBB)
	Δ = obb.maxb-obb.minb
	return reduce(*,Δ)
end

#
function pos(P::Union{Point3f0,Vector{Float32}},A,n)
	return dot(P-A,n)>0
end

function pos(P::Vector{Point3f0},A,n)
	return all(map(x->pos(x,A,n),P))
end
#

function inside(obb::OBB, x::Union{Point3f0,Vector{Float32}})
	#return all(obb.minb.<(x+Float32(eps()))) && all(obb.maxb.>(x-Float32(eps())))
	Δ = obb.maxb-obb.minb
	return all(obb.minb.<(x+Δ/100)) && all(obb.maxb.>(x-Δ/100))
end

function inside(obb::OBB, x::Vector{Point3f0})
	return sum(map(xx->inside(obb,xx),x))>1
end

function Base.isvalid(obb::OBB)
	return all(obb.minb .< obb.maxb)
end

function inside_vertices(obb::OBB, vertices)
	inds = [inside(obb,inv(obb.pca)*vert) for vert in vertices]
	return inds,vertices[inds]
end

function inside_faces(obb::OBB, vertices, faces)
	inds = [inside(obb,[map(x->inv(obb.pca)*x,vertices[face])...]) for face in faces]
	return inds,faces[inds]
end

function oriented_bounding_box(obb::OBB,vertices,faces,vinds,finds)
	isd_fcs_inds,isd_fcs = inside_faces(obb,vertices,faces[finds])
	faceverts = unique(vcat(convert(Vector{Vector{Int32}},faces[finds[isd_fcs_inds]])...))
	isd_verts_inds,isd_verts = inside_vertices(obb,vertices[faceverts])
	if length(faceverts) > 3 && length(isd_fcs) > 1
		pca, wc = pca_weighted(vertices,isd_fcs)
		#newobb = oriented_bounding_box(isd_verts,pca)
		newobb = oriented_bounding_box(vertices[faceverts],pca)
	else
		info("length(faceverts) ")
		display(length(faceverts))
		info("length(isd_fcs) ")
		display(length(isd_fcs))
		newobb = OBB(obb.pca)
		wc = Point3f0(0)
	end
	return newobb,wc,vinds,finds[isd_fcs_inds]
end

function subdivide(obb::OBB, wc, direction)
	obb1,obb2 = copy(obb),copy(obb)
	len = (obb.maxb[direction]-obb.minb[direction])
	Δ = (inv(obb.pca)*wc)[direction] - obb.minb[direction]
	#Δ = len/2
	obb1.maxb[direction] -= (len-Δ)
	obb2.minb[direction] += Δ
	return obb1,obb2
end

function GeometryTypes.volume(elem::Tuple{Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}},Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}}})
	return volume(elem[1][1])+volume(elem[2][1])
end

function get_children(obbs::Vector{Tuple{Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}},Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}}}})
	#return obbs[1]
	#return obbs[findmin(map(x->volume(x),obbs))[2]]
	val = findmin(map(x->volume(x),obbs[1:2]))
	if abs(val[1]) < Inf32
		return obbs[val[2]]
	else
		return obbs[3]
	end
end

function refine(obb::OBB,vertices,faces,wc,vinds,finds,obbtree)
	directions = sortperm(obb.maxb-obb.minb,rev=true)
	#directions = directions[[1]]
	obbs = Vector{Tuple{Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}},Tuple{OBB,Point3f0,Vector{Int32},Vector{Int32}}}}()
	for direction in directions
		obb1,obb2 = subdivide(obb,wc,direction)

		#isd_fcs_inds1,isd_fcs1 = inside_faces(obb1,vertices,faces[finds])
		#isd_verts_inds1,isd_verts1 = inside_vertices(obb1,vertices[vinds])
		#isd_fcs_inds2,isd_fcs2 = inside_faces(obb2,vertices,faces[finds])
		#isd_verts_inds2,isd_verts2 = inside_vertices(obb2,vertices[vinds])

		obb1tmp,wc1tmp,vinds1tmp,finds1tmp = oriented_bounding_box(obb1,vertices,faces,vinds,finds)
		obb2tmp,wc2tmp,vinds2tmp,finds2tmp = oriented_bounding_box(obb2,vertices,faces,vinds,finds)

		#obbtreeelement1 = OBBTreeElement(1,obb1,Nullable{OBBTreeElement}(),wc1tmp,vinds[isd_verts_inds1],finds[isd_fcs_inds1])
		#obbtreeelement2 = OBBTreeElement(1,obb2,Nullable{OBBTreeElement}(),wc2tmp,vinds[isd_verts_inds2],finds[isd_fcs_inds2])
		#println("OBB Mesh 1")
		#plot("OBB mesh 1", obbtree, obbtreeelement1, rand(RGBA{Float32}))
		#println("OBB Mesh 2")
		#plot("OBB mesh 2", obbtree, obbtreeelement2, rand(RGBA{Float32}))
		#println("OBB 1")
		#plot("Bounding Box 1", obbtreeelement1.obb, rand(RGBA{Float32}))
		#println("OBB 2")
		#plot("Bounding Box 2", obbtreeelement2.obb, rand(RGBA{Float32}))
		#println("done")
		#sleep(.1)

		obb1,wc1,vinds1,finds1 = obb1tmp,wc1tmp,vinds1tmp,finds1tmp
		obb2,wc2,vinds2,finds2 = obb2tmp,wc2tmp,vinds2tmp,finds2tmp

		push!(obbs,((obb1,wc1,vinds1,finds1),(obb2,wc2,vinds2,finds2)))
	end
	return get_children(obbs)
end

abstract type AbstractTreeElement end
type OBBTreeElement <: AbstractTreeElement
	level::Int32
	obb::OBB
	children::Nullable{Tuple{AbstractTreeElement,AbstractTreeElement}}
	parent::Nullable{AbstractTreeElement}
	wc::Point3f0
	vinds::Vector{Int32}
	finds::Vector{Int32}
	OBBTreeElement(level,obb,wc) = new(level,obb,Nullable{Tuple{AbstractTreeElement,AbstractTreeElement}}(),Nullable{AbstractTreeElement}(),wc,Vector{Int32}(),Vector{Int32}())
	OBBTreeElement(level,obb,parent,wc) = new(level,obb,Nullable{Tuple{AbstractTreeElement,AbstractTreeElement}}(),Nullable{AbstractTreeElement}(parent),wc,Vector{Int32}(),Vector{Int32}())
	OBBTreeElement(level,obb,parent,wc,vinds,finds) = new(level,obb,Nullable{Tuple{AbstractTreeElement,AbstractTreeElement}}(),Nullable{AbstractTreeElement}(parent),wc,vinds,finds)
end

type OBBTree
	maxdepth::Int32
	tree::Vector{Set{OBBTreeElement}}
	vertices::Vector{Point3f0}
	faces::Vector{Face{3,AInt}}
	OBBTree(maxdepth,vertices,faces) = new(maxdepth,Vector{Set{OBBTreeElement}}(maxdepth),vertices,faces)
end

function Base.show(io::IO, obt::OBBTree)
    print(io,"OBBTree")
end

function Base.show(io::IO, obt::OBBTreeElement)
    print(io,"OBBTreeElement Level $(obt.level)")
end

function GeometryTypes.volume(elem::OBBTreeElement)
	Δ = elem.obb.maxb-elem.obb.minb
	return reduce(*,Δ)
end

function Base.push!(obbt::OBBTree, elem::OBBTreeElement)
	push!(obbt.tree[elem.level],elem)
end

function intialize!(obbtree::OBBTree)
	for i = 1:obbtree.maxdepth; obbtree.tree[i] = Set{OBBTreeElement}(); end
	pca, wc = pca_weighted(obbtree.vertices,obbtree.faces)
	obb = oriented_bounding_box(obbtree.vertices,pca)
	elem1 = OBBTreeElement(1,obb,wc)
	elem1.vinds = collect(Int32(1):Int32(length(obbtree.vertices)))
	elem1.finds = collect(Int32(1):Int32(length(obbtree.faces)))
	push!(obbtree,elem1)
	#refine!(obbtree,elem1)
	return nothing
end

function refine!(obbtree::OBBTree, elem::OBBTreeElement)
	#((obb1,wc1,vinds1,finds1),(obb2,wc2,vinds2,finds2)) = refine(elem.obb,obbtree.vertices,obbtree.faces,elem.wc,elem.vinds,elem.finds)
	((obb1,wc1,vinds1,finds1),(obb2,wc2,vinds2,finds2)) = refine(elem.obb,obbtree.vertices,obbtree.faces,elem.wc,elem.vinds,elem.finds,obbtree)
	if isvalid(obb1) && isvalid(obb2)
		child1 = OBBTreeElement(elem.level+1,obb1,elem,wc1,vinds1,finds1)
		child2 = OBBTreeElement(elem.level+1,obb2,elem,wc2,vinds2,finds2)
		elem.children = Nullable{Tuple{OBBTreeElement,OBBTreeElement}}((child1,child2))
		push!(obbtree,child1)
		push!(obbtree,child2)
	else
		child1 = OBBTreeElement(elem.level+1,obb1,elem,wc1,vinds1,finds1)
		child2 = OBBTreeElement(elem.level+1,obb2,elem,wc2,vinds2,finds2)
	end
	return child1,child2
end

function get_highest_level_descendants(obbtr::OBBTree)
	ret = Set{OBBTreeElement}()
	tmp = Set{OBBTreeElement}()
	push!(tmp,first(obbtr.tree[1]))
	while !isempty(tmp)
		elem = pop!(tmp)
		if elem.children.hasvalue
			union!(tmp,elem.children.value)
		else
			push!(ret,elem)
		end
	end
	return ret
end

function Base.isvalid(elem::OBBTreeElement)
	return isvalid(elem.obb)
end

function Base.delete!(obbtree::OBBTree, child1::OBBTreeElement, child2::OBBTreeElement)
	if !child1.children.hasvalue && !child2.children.hasvalue
		delete!(obbtree.tree[child1.level],child1)
		delete!(obbtree.tree[child2.level],child2)
		child1.parent.value.children = Nullable{Tuple{OBBTreeElement,OBBTreeElement}}()
		child1.parent = Nullable{OBBTreeElement}()
		child2.parent = Nullable{OBBTreeElement}()
		return true
	else
		return false
	end
end

function adaptive_refine!(obbtr::OBBTree, tol::Float32)
	ret = get_highest_level_descendants(obbtr)
	while !isempty(ret)
		elem = pop!(ret)
		if !elem.children.hasvalue && elem.level < obbtr.maxdepth
			child1,child2 = refine!(obbtr,elem)
			println((volume(child1)+volume(child2))/volume(elem))
			if isvalid(child1) && isvalid(child2) && (volume(child1)+volume(child2))/volume(elem) < tol
				info("Refinement on Level: $(child1.parent.value.level)")
				push!(ret,child1)
				push!(ret,child2)
			else
				delete!(obbtr,child1,child2)
			end
		end
	end
end

function collect_parents!(ret::Set{OBBTreeElement},elem::OBBTreeElement)
	push!(ret,elem)
	if elem.parent.hasvalue
		collect_parents!(ret,elem.parent.value)
	end
	return nothing
end

function collect_parents(elem::OBBTreeElement)
	ret = Set{OBBTreeElement}()
	collect_parents!(ret,elem)
	return ret
end

function plot_parents(obbtree::OBBTree,leaf::OBBTreeElement)
	ret = collect_parents(leaf)
	for (i,elem) in enumerate(ret)
		plot("OBB fitted $i",elem.obb,RGBA{Float32}(rand(),rand(),rand(),.4))
	end
end

function plot(obbtr::OBBTree)
	ret = get_highest_level_descendants(obbtr)
	#colors = linspace(RGBA{Float32}(1,0,0,1),RGBA{Float32}(0,0,1,1),length(ret))
	colors = rand(RGB{Float32},length(ret))
	for (i,elem) in enumerate(ret)
		plot("OBB fitted $i",elem.obb,RGBA{Float32}(colors[i].r,colors[i].g,colors[i].b,.4))

		#sleep(.1)
		#plot("OBB wc $i",elem.obb.pca, elem.wc)

		plot("OBB mesh fitted $i", obbtr, elem, RGBA{Float32}(colors[i].r,colors[i].g,colors[i].b,.4))
	end
end

import Colors

function vertices_and_faces(vertices,faces,finds)
	ins_fcs = faces[finds]
	ins_fcs_new = similar(faces[finds])
	vertinds = convert(Vector{Vector{Int32}},ins_fcs)
	unq_vert_inds = unique(vcat(vertinds...))
	verts = vertices[unq_vert_inds]
	for i = 1:length(ins_fcs_new)
		ins_fcs_new[i] = map(y->findfirststop(x->x==ins_fcs[i][y],unq_vert_inds),1:3)
	end
	return verts,ins_fcs_new
end

function plot(tag::String,obbtree::OBBTree,elem::OBBTreeElement,color::RGBA{Float32})
	#test = [inside(elem.obb,[map(x->inv(elem.obb.pca)*x,obbtree.vertices[face])...]) for face in obbtree.faces]
	verts,fcs = vertices_and_faces(obbtree.vertices,obbtree.faces,elem.finds)
	#colors = [color for i = 1:length(obbtree.vertices)]
	colors = [color for i = 1:length(verts)]
	options = ["surface_shader = transparent"]
	unity_faces = vcat(map(xx->convert(Vector{UInt32},xx-1),fcs)...)
	utm = UnityMesh(tag,verts,Int32[],Int32[],unity_faces,colors,options)
	socket = connect(8052)
	write(socket,utm)
	close(socket)
	#sleep(.1)
end

function plot(tag::String,obb::OBB, color::RGBA{Float32}=RGBA{Float32}(1,0,0,1))
	vertices = boundingbox_vertices(obb.minb,obb.maxb)
	info("boundingbox size: ",volume(obb))
	lines =  [0,1,0,2,0,3,3,6,3,5,2,6,1,4,5,7,1,5,4,2,4,7,7,6]
	colors = Colors.linspace(RGBA{Float32}(1,1,1,1),color,8)
	verts = map(x->obb.pca*x,vertices)
	#println(obb.minb)
	#println(obb.maxb)
	#utm = UnityMesh(tag, verts, Int32[0,7], lines, Int32[], colors, [" = "])
	utm = UnityMesh(tag, verts, Int32[], lines, Int32[], colors, [" = "])
	socket = connect(8052)
	#println(socket)
	a = write(socket,utm)
	#println(a)
	close(socket)
end

function plot(points::Vector{Point3f0})
	red = RGBA{Float32}(1,0,0)
	unity_point_mesh = UnityMesh("Point Mesh",points,collect(0:length(points)-1), Int32[], Int32[], [red for i=1:length(points)], [" = "])
	socket = connect(8052)
	write(socket,unity_point_mesh)
	close(socket)
end

function plot(tag::String, pca::StaticArrays.SArray{Tuple{3,3},Float32,2,9}, weighted_center::GeometryTypes.Point{3,Float32})
	_pca = Point3f0[pca[:,i] for i = 1:3]
	plot(tag,_pca,weighted_center)
end

function plot(tag::String,pca::Vector{Point3f0}, weighted_center::Point3f0)
	black = RGBA{Float32}(0,0,0)
	red = RGBA{Float32}(1,0,0)
	green = RGBA{Float32}(0,1,0)
	blue = RGBA{Float32}(0,0,1)
	vertices = vcat([weighted_center],map(x->weighted_center+1*x,pca))
	unity_mesh = UnityMesh(tag,vertices, Int32[0,1,2,3], Int32[0,1,0,2,0,3], Int32[], [black, red, green, blue ], ["point_size = 0.02"])
	socket = connect(8052)
	write(socket,unity_mesh)
	close(socket)
end

function test()
	#vertices,faces = vertices_and_faces(mayo_mesh)
	vertices,_faces = vertices_and_faces(bone_mesh)
	pca, wc = pca_weighted(vertices,_faces)

	centered_vertices = map(x->x-wc,vertices)
	faces = Vector{GeometryTypes.Face{3,Int32}}()
	for face in _faces
		if all(map(x->x[3],centered_vertices[face]).>30)
			push!(faces,face)
		end
	end
	#centered_mesh = UnityMesh("Mayo Center",centered_vertices,Int32[],Int32[],mayo_triangles,mvayo_color,mayo_options)
	centered_mesh = UnityMesh("Mayo Center",centered_vertices,Int32[],Int32[],bone_triangles,bone_color,["surface_shader = transparent"])
	plot(centered_mesh)

	coords = [Point3f0(1,0,0),Point3f0(0,1,0),Point3f0(0,0,1)]
	plot("Coordinate System",coords,Point3f0(0,0,0))

	obbtree = OBBTree(5,centered_vertices,faces)
	intialize!(obbtree)
	adaptive_refine!(obbtree,1f0)
	plot(obbtree)
	return obbtree
end

function plot(um::UnityMesh)
	socket = connect(8052)
	write(socket,um)
	close(socket)
end

function findfirststop{T<:Function}(f::T,vec::Vector)
	for (i,v) in enumerate(vec)
	(f(v) == true) ? (return i) : (nothing)
	end
	return -1
end

function d(mesh)
	vertices,boundary_faces = vertices_and_faces(mesh)
	vertinds = convert(Vector{Vector{Int32}},boundary_faces)
	unq_vert_inds = unique(vcat(vertinds...))
	boundary_vertices = vertices[unq_vert_inds]
	for i = 1:length(boundary_faces)
		boundary_faces[i] = map(y->findfirststop(x->x==boundary_faces[i][y],unq_vert_inds),1:3)
	end
	return boundary_vertices,boundary_faces
end

function simpletest2()
	file1 = "/Users/maximilianbittens/.julia/v0.6/Unity/models/gmsh/box.msh"
	vertices,eldat = AnsysReader.read(file1)
	triangles = vcat(convert(Vector{Vector{Int32}},eldat)...)-Int32(1)
	options = ["surface_shader = transparent"]
	color = [RGBA{Float32}(1,0,0,1) for i = 1:length(vertices)]
	unity_mesh = UnityMesh("Mesh",vertices,Int32[],Int32[],triangles,color,options)
	plot(unity_mesh)
end

function simpletest()
	file1 = "/Users/maximilianbittens/.julia/v0.6/Unity/models/gmsh/box.msh"
	box_mesh = MeshReaders.load(MeshReaders.GmshReader, file1)
	unity_triangles = vcat(map(x->x[4:end],filter(x->x[1:3]==[2,3,1],box_mesh.elements))-1...)
	vertices = Point3f0[Point3f0(box_mesh.nodes[1,i],box_mesh.nodes[2,i],box_mesh.nodes[3,i]) for i = 1:size(box_mesh.nodes,2)]
	options = ["surface_shader = transparent","point_size = 0.02"]
	color = [RGBA{Float32}(1,1,1,1) for i = 1:length(vertices)]
	unity_mesh = UnityMesh("Mesh",vertices,collect(UInt32(0):UInt32(length(vertices)-1)),Int32[],unity_triangles,color,options)
	plot(unity_mesh)
	vertices,faces = vertices_and_faces(box_mesh)
	obbtree = OBBTree(5,vertices,faces)
	intialize!(obbtree)
	elem = first(obbtree.tree[1])
	plot("elem coords",elem.obb.pca, elem.wc)
	pca,wc = pca_weighted(vertices,faces)
	plot("elem coords2",pca, wc)
	#child1,child2 = refine!(obbtree,elem)
	plot(obbtree)
end

# load mesh
file1 = "/Users/maximilianbittens/.julia/v0.6/Unity/models/gmsh/Box.msh"
box_mesh = MeshReaders.load(MeshReaders.GmshReader, file1)
# get boundary
boundary_vertices,boundary_faces = d(box_mesh)
# get all vertices and faces
vertices,faces = vertices_and_faces(box_mesh)
# get weighted center
pca,wc = pca_weighted(vertices,faces)
#shift vertices and boundary_vertices
#vertices = map(x->x-wc,vertices)
#boundary_vertices = map(x->x-wc,boundary_vertices)

# plot mesh
unity_triangles = vcat(convert(Vector{Vector{Int32}},faces)...)-Int32(1)
options = ["surface_shader = transparent","point_size = 0.02"]
color = [RGBA{Float32}(.5,.5,.5,.3) for i = 1:length(vertices)]
unity_mesh = UnityMesh("Mesh",vertices,collect(UInt32(0):UInt32(length(vertices)-1)),Int32[],unity_triangles,color,options)
plot(unity_mesh)

# plot coordinate system
coords = [Point3f0(1,0,0),Point3f0(0,1,0),Point3f0(0,0,1)]
plot("Coordinate System",coords,Point3f0(0,0,0))

#pcac,wcc = pca_weighted(boundary_vertices,boundary_faces)
#plot("elem coordsc",pcac, wcc+Point3f0(.1,0,0))

obbtree = OBBTree(3,vertices,faces)
intialize!(obbtree)
elem = first(obbtree.tree[1])
#plot("elem coords",elem.obb.pca, elem.wc)
adaptive_refine!(obbtree,2.5f0)
plot(obbtree)
elem = first(filter(!isempty,obbtree.tree)[end])
#plot_parents(obbtree,elem)

#=
pca,wc = pca_weighted(boundary_vertices, boundary_faces[child1.finds])
obbnew = oriented_bounding_box(boundary_vertices[child1.vinds],pca)
obbtreeelem = OBBTreeElement(1,obbnew,Nullable{OBBTreeElement}(),wc,child1.vinds,child1.finds)

println("OBB Mesh new")
plot("OBB meshnew", obbtree, obbtreeelem, rand(RGBA{Float32}))
plot("Bounding Box new", obbtreeelem.obb, rand(RGBA{Float32}))
plot("elem coordsc",pca, wc)
#simpletest()
=#


text:
module bicycle
import Parameters: @with_kw
import MyterialColors: blue_grey_darker, blue_grey, cyan_dark

import jcontrol: closest_point_idx, Track, euclidean, movingaverage

export Bicycle, State

"""
Stores immutable geometric properties of the bike.
Also stores parameters for drawing.
"""
struct Bicycle
    # geometry
    l_f::Number     # distance from COM to front wheel | cm
    l_r::Number     # distance from rear wheel to COM | cm
    L::Number       # total length | cm
    width::Number   # 'width' of bike. Used to staty within track | cm

    # dynamics
    m::Number       # mass | g
    Iz::Number      # moment of angular inertia | Kg⋅m²
    c::Number       # corenring stiffness

    # for drawing
    wheel_length::Number
    wheel_lw::Number
    wheel_color::String
    front_wheel_color::String
    body_lw::Number
    body_color::String

    function Bicycle(;
        l_f::Number=3,
        l_r::Number=2,
        width::Number=2,
        m_f=12,
        m_r=11,
        c=4e3
        )

        # convert units g->Kg, cm->m
        mfKg = m_f # / 100
        mrKg = m_r # / 100
        lfM = l_f # / 100
        lrM = l_r # / 100

        # compute moment of angular inertia
        Iz = mfKg * lfM^2 + mrKg * lrM^2

        return new(
            l_f,
            l_r,
            l_f + l_r,
            width,
            m_f + m_r,
            Iz,
            c,
            0.8,                # wheel_length
            16,                 # wheel_lw
            blue_grey_darker,   # wheel_color
            cyan_dark,          # front_wheel_color
            8,                  # body_lw
            blue_grey,          # body_color
        )
    end
end

"""
Represents the state of the bicycle model at a moment in time.
Can be used to pass initial and final condistions to the control
model.
"""
@with_kw mutable struct State
    x::Number = 0  # position
    y::Number = 0
    θ::Number = 0  # orientation
    δ::Number = 0  # steering angle
    δ̇::Number = 0
    ω::Number = 0
    u::Number = 0  # velocity  | fpr DynamicsProblem its the longitudinal velocity component

    # KinematicsProblem only
    β::Number = 0  # slip angle

    # DynamicsProblem only
    v::Number = 0  # lateral velocity
    Fu::Number = 0  # forward force

    # track errors
    n::Number = 0
    ψ::Number = 0

    # track variables
    track_idx = 0
    s = 0
end

"""
    State(trial::Trial, frame::Int, track::Track)

Get `State` from experimental data at a frame
"""
function State(trial, frame::Int, track::Track; smoothing_window=1, kwargs...)
    if smoothing_window == 1
        _x, _y, _ω, _u, _θ = trial.x, trial.y, trial.ω, trial.u, trial.θ
    else
        _x = movingaverage(trial.x, smoothing_window)
        _y = movingaverage(trial.y, smoothing_window)
        _ω = movingaverage(trial.ω, smoothing_window)
        _u = movingaverage(trial.u, smoothing_window)
        _θ = movingaverage(trial.θ, smoothing_window)
    end

    # get track errors
    x, y = _x[frame], _y[frame]
    dist = sqrt.((track.X.-x).^2 + (track.Y.-y).^2)
    idx = argmin(dist)

    n = dist[idx]

    # get the sign of n right based on which side of the track the mouse is
    if idx > 2 && idx < length(track.X) - 2
        x1, y1, x2, y2 = track.X[idx-1], track.Y[idx-1], track.X[idx+1], track.Y[idx+1]
        d = (x - x1)*(y2 - y1) - (y - y1)*(x2 - x1)
        n *= sign(d)
    end

    ψ = track.θ[idx] - _θ[frame]
    # ψ = mod(abs(ψ), 2π)

    return State(; x=x, y=y, θ=_θ[frame], ω=_ω[frame], u=_u[frame], n=n, ψ=ψ, track_idx=idx, s=track.S[idx], kwargs...)
end

"""
    State(solution, frame::Int)

Get `State` from a forward problem solution at Δt from start.
"""
function State(solution, Δt::Float64)
    frame = (Int ∘ round)(Δt / solution.δt)

    return State(;
        x=solution.x[frame],
        y=solution.y[frame],
        θ=solution.θ[frame],
        δ=solution.δ[frame],
        δ̇=solution.δ̇[frame],
        ω=solution.ω[frame],
        u=solution.u[frame],
        β=solution.β[frame],
        v=solution.v[frame],
        n=solution.n[frame],
        ψ=solution.ψ[frame],
        Fu = solution.Fu[frame]
    )
end
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
using JLD
using BenchmarkTools
using Dates
using Base.Sys
using CellListsBenchmarks
using Plots

gettime(x) = getfield(x, :time)
gettime_nogc(x) = getfield(x, :time) - getfield(x, :gctime)

files = [
    "n100-d2-r0.01.jld",
    "n100-d3-r0.01.jld",
    "n200-d4-r0.01.jld",
    "n300-d5-r0.01.jld"
]
dir_cell_lists = "output/cell_lists_2021-04-09T08:28:29+0300/"

function foo(file)
    lc = load(joinpath(dir_cell_lists, file))["benchmarks"]
    d = getfield(lc[1], :d)
    r = getfield(lc[1], :r)
    ns = getfield.(lc, :n)
    tc = getfield.(lc, :trials)
    mc = [median.(t) for t in tc]
    meanc = [mean(gettime.(v)) for v in mc]
    return ns, meanc, d, r
end

vs = [foo(file) for file in files]
ns, y, d, r = vs[1]

p1 = plot(legend=:topleft)
m = ns .≤ 100
plot!(p1, ns[m], y[m], markershape=:circle, markersize=2, linestyle=:dash, label="d=$d")
for (ns2, y2, d2, r2) in vs[2:end]
    m2 = ns2 .≤ 100
    plot!(p1, ns2[m2], y2[m2], markershape=:circle, markersize=2, linestyle=:dash, label="d=$d2")
end

p2 = plot(legend=:topleft)
m = ns .≤ 100
plot!(p2, ns[m], y[m]./y[m], markershape=:circle, markersize=2, linestyle=:dash, label="d=$d")
for (ns2, y2, d2, r2) in vs[2:end]
    m2 = ns2 .≤ 100
    plot!(p2, ns2[m2], y2[m2]./y[m], markershape=:circle, markersize=2, linestyle=:dash, label="d=$d2")
end

directory = joinpath("figures", "dimensionality")
if !isdir(directory)
    mkpath(directory)
end

savefig(p1, joinpath(directory, "absolute.svg"))
savefig(p2, joinpath(directory, "ratio.svg"))

max_stars_repo_name: jaantollander/cell-lists-benchmarks

max_issues_repo_name: jaantollander/cell-lists-benchmarks

max_forks_repo_name: jaantollander/cell-lists-benchmarks


text:
using DualDecomposition
using JuMP, Ipopt, GLPK

const DD = DualDecomposition

const NS = 3  # number of scenarios
const probability = ones(3) / 3

const CROPS = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
const PURCH = 1:2 # set of crops to purchase (wheat and corn, resp.)
const SELL = 1:4  # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)

const Cost = [150 230 260]    # cost of planting crops
const Budget = 500            # budget capacity
const Purchase = [238 210]    # purchase price
const Sell = [170 150 36 10]  # selling price
const Yield = [3.0 3.6 24.0; 2.5 3.0 20.0; 2.0 2.4 16.0]
const Minreq = [200 240 0]    # minimum crop requirement

# This creates a Lagrange dual problem for each scenario s.
function create_scenario_model(s::Int64)
    m = Model(GLPK.Optimizer)
    @variable(m, 0 <= x[i=CROPS] <= 500, Int)
    @variable(m, y[j=PURCH] >= 0)
    @variable(m, w[k=SELL] >= 0)

    @objective(m, Min,
        probability[s] * sum(Cost[i] * x[i] for i=CROPS)
        + probability[s] * sum(Purchase[j] * y[j] for j=PURCH)
        - probability[s] * sum(Sell[k] * w[k] for k=SELL))

    @constraint(m, sum(x[i] for i=CROPS) <= Budget)
    @constraint(m, [j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
    @constraint(m, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
    @constraint(m, w[3] <= 6000)
    return m
end

# Create DualDecomposition instance.
algo = DD.LagrangeDual()

# Add Lagrange dual problem for each scenario s.
models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in 1:NS)
for s in 1:NS
    DD.add_block_model!(algo, s, models[s])
end

coupling_variables = Vector{DD.CouplingVariableRef}()
for s in 1:NS
    model = models[s]
    xref = model[:x]
    for i in CROPS
        push!(coupling_variables, DD.CouplingVariableRef(s, i, xref[i]))
    end
end

# Set nonanticipativity variables as an array of symbols.
DD.set_coupling_variables!(algo, coupling_variables)

# Lagrange master method
LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))

# Solve the problem with the solver; this solver is for the underlying bundle method.
DD.run!(algo, LM)


text:
module TestDeprecate

using Test
using VecMathOptInterface
const MOI = VecMathOptInterface

function test_deprecations_N()
    attr = MOI.VariablePrimal()
    @test_deprecated begin
        @test MOI._result_index_field(attr) == 1
    end
    for attr_type in (
        MOI.ObjectiveValue,
        MOI.DualObjectiveValue,
        MOI.VariablePrimal,
        MOI.ConstraintPrimal,
        MOI.ConstraintDual,
        MOI.ConstraintBasisStatus,
        MOI.PrimalStatus,
        MOI.DualStatus,
        MOI.NLPBlockDual,
    )
        attr = attr_type(1)
        @test_logs(
            (:warn, "Field attr.N is deprecated, use attr.result_index"),
            attr.N,
        )
    end
end

function test_deprecations_ScalarAffineTerm()
    x = MOI.VariableIndex(1)
    t = MOI.ScalarAffineTerm(1.0, x)
    @test_logs (:warn,) t.variable_index == x
end

function test_deprecations_ScalarQuadraticTerm()
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    t = MOI.ScalarQuadraticTerm(1.0, x, y)
    @test_logs (:warn,) t.variable_index_1 == x
    @test_logs (:warn,) t.variable_index_2 == y
end

function test_deprecations_ScalarQuadraticFunction()
    @test_logs(
        (:warn,),
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm{Int}[],
            MOI.ScalarQuadraticTerm{Int}[],
            0,
        ),
    )
    return
end

function test_deprecations_VectorQuadraticFunction()
    @test_logs(
        (:warn,),
        MOI.VectorQuadraticFunction(
            MOI.VectorAffineTerm{Int}[],
            MOI.VectorQuadraticTerm{Int}[],
            [0, 1],
        ),
    )
    return
end

function test_RawOptimizerAttribute()
    @test_logs (:warn,) MOI.RawParameter(:a) == MOI.RawOptimizerAttribute("a")
end

function test_default_copy_to()
    dest = MOI.Utilities.Model{Float64}()
    src = MOI.Utilities.Model{Float64}()
    @test_logs (:warn,) MOI.Utilities.default_copy_to(dest, src, true)
    return
end

function test_copy_to_copy_names()
    dest = MOI.Utilities.Model{Float64}()
    src = MOI.Utilities.Model{Float64}()
    @test_logs (:warn,) MOI.copy_to(dest, src; copy_names = true)
    return
end

function test_test_models_equal()
    dest = MOI.Utilities.Model{Float64}()
    src = MOI.Utilities.Model{Float64}()
    @test_logs(
        (:warn,),
        MOI.Utilities.test_models_equal(dest, src, String[], String[]),
    )
    return
end

function test_IndexMap()
    @test_logs (:warn,) MOI.IndexMap(1)
    return
end

function test_CleverDicts()
    K, V = MOI.VariableIndex, MOI.VariableIndex
    @test_logs (:warn,) MOI.Utilities.CleverDicts.CleverDict{K,V}(1)
    @test_logs (:warn,) MOI.Utilities.CleverDicts.CleverDict{K,V}(
        MOI.Utilities.CleverDicts.key_to_index,
        MOI.Utilities.CleverDicts.index_to_key,
        1,
    )
    return
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestDeprecate.runtests()


text:
# --------------
# Abstract types
# --------------

@doc """
# Interface

* reward(R, s, a)

## Optional

* reward(R, s, s, a)
* reward(R, s)

""" ->
abstract AbstractReward

abstract AbstractArrayReward <: AbstractReward

getindex(R::AbstractArrayReward, dims...) = getindex(R.array, dims...)

size(R::AbstractArrayReward, d) = size(R.array, d)

num_states(R::AbstractArrayReward) = size(R, 1)

# ----------
# Array type
# ----------

immutable ArrayReward{T<:Real,N} <: AbstractArrayReward
    array::Array{T,N}

    function ArrayReward(array::Array{T,N})
        N < 3 || (N == 3 && size(array, 1) == size(array, 2)) ||
            error("ArrayReward can not constructed with this Array.")
        new(array)
    end
end

ArrayReward{T,N}(array::Array{T,N}) = ArrayReward{T,N}(array)

# type construction helper
Reward(A::Array) = ArrayReward(A)

num_actions{T}(R::ArrayReward{T,1}) =
    error(string(typeof(R))*" does not know about actions.")
num_actions{T,N}(R::ArrayReward{T,N}) = size(R, N)

# 1 dimension
# -----------

reward{T}(R::ArrayReward{T,1}, state) = getindex(R, state)

reward{T}(R::ArrayReward{T,1}, state, action) = reward(R, state)

reward{T}(R::ArrayReward{T,1}, state, new_state, action) = reward(R, state)

# 2 dimensions
# ------------

reward{T}(R::ArrayReward{T,2}, state, action) = getindex(R, state, action)

reward{T}(R::ArrayReward{T,2}, state, new_state, action) = reward(R, state, action)

# 3 dimensions
# ------------

reward{T}(R::ArrayReward{T,3}, state, new_state, action) =
    getindex(R, state, new_state, action)

# -----------
# Sparse type
# -----------

immutable SparseReward{Tv<:Real,Ti} <: AbstractArrayReward
    array::SparseMatrixCSC{Tv,Ti}
end

Reward(A::SparseMatrixCSC) = SparseReward(A)

num_actions(R::SparseReward) = size(R, 2)

reward(R::SparseReward, state, action) = getindex(R, state, action)

reward(R::SparseReward, state, new_state, action) = reward(R, state, action)


text:
module TensorNetworkStatesFunctions
    include("./TensorNetworkStatesGates.jl")
    include("./MPSforQuantum.jl")
    using DataStructures
    using Reexport
    @reexport using .TensorNetworkStatesGates
    using .MPSforQuantum
    export apply!, restore, init_zero_MPS, gate_parser, mps_size

    function init_zero_MPS(N::Int64)
        mps = []
        for i in 1:N
            push!(mps, [reshape([1.0 + 0.0im], 1, 1), reshape([0.0 + 0.0im], 1, 1)])
        end
        return mps
    end

    function init_MPS(state::Array{ComplexF64,1})
        mps = MPS(state, param, 'l')
        return mps
    end

    # Restore prob. ampl. of a state from MPS
    function restore(arrs::Array{Any, 1}, n::Int64) # n is decimal representation of the state '0110...'
        N = Int64(size(arrs)[1])
        s = bitstring(n)[end - (N - 1):end] # 後ろからN bit目まで
        phys_idx = [convert(Int64, s[i]) - 48 for i = length(s):-1:1] # 後ろが最上位ビット
        return prod([arrs[i][phys_idx[i] + 1] for i = 1:length(phys_idx)])[1]
    end

    function mps_size(arrs::Array{Any, 1})
        d = OrderedDict()
        params = 0
        N = size(arrs)[1]
        for i = 1:N
            #println("array $i's size: ", size(arrs[i][1]))
            d["array $i's size"] = size(arrs[i][1])
            params += length(arrs[i][1]) * 2
        end
        #println("Num of parameters: ", params)
        d["Num of parameters"] = params
        #println("2^N: ", 2^N)
        d["2^N"] = 2^N
        return d
    end

    """
    X gate
    """
    function apply!(gate::X, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [0 1; 1 0])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    Y gate
    """
    function apply!(gate::Y, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [0 -1im; 1im 0])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    Z gate
    """
    function apply!(gate::Z, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [1 0; 0 -1])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    H gate
    """
    function apply!(gate::H, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [1 1; 1 -1] ./ sqrt(2))
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    S gate
    """
    function apply!(gate::S, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [1 0; 0 1im])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    T gate
    """
    function apply!(gate::T, mps::Array{Any, 1}, param::Number)
        target = gate._target
        mpo = convert(Array{ComplexF64,2}, [1 0; 0 _exp(0.25im * pi)])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    RX gate
    """
    function apply!(gate::RX, mps::Array{Any, 1}, param::Number)
        target = gate._target
        angpi = gate._theta * 0.5 / pi
        a = cospi(angpi)
        b = -1im*sinpi(angpi)
        mpo = convert(Array{ComplexF64,2}, [a b; b a])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    RZ gate
    """
    function apply!(gate::RZ, mps::Array{Any, 1}, param::Number)
        target = gate._target
        ang = gate._theta * 0.5
        mpo = convert(Array{ComplexF64,2}, [_exp(-1im*ang) 0; 0 _exp(1im*ang)])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    RY gate
    """
    function apply!(gate::RY, mps::Array{Any, 1}, param::Number)
        target = gate._target
        angpi = gate._theta * 0.5 / pi
        a = cospi(angpi)
        b = sinpi(angpi)
        mpo = convert(Array{ComplexF64,2}, [a -b; b a])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    U3 gate
    """
    function apply!(gate::U3, mps::Array{Any, 1}, param::Number)
        target = gate._target
        theta = gate._theta * 0.5
        thetapi = theta / pi
        phipi = gate._phi / pi
        lambdpi = gate._lambd / pi
        a = cospi(thetapi)
        b = -_exp(1im*lambpi)*sinpi(thetapi)
        c = _exp(1im*phipi)*sinpi(thetapi)
        d = _exp(1im*(phipi+lambpi))*cospi(thetapi)
        mpo = convert(Array{ComplexF64,2}, [a b; c d])
        mps_ = OneQubitGate(mps, mpo, target)
    end

    """
    CX gate
    """
    function apply!(gate::CX, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        _arr = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]
        mpo = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, mpo, control, target, param)
    end

    """
    CZ gate
    """
    function apply!(gate::CZ, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        _arr = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 -1]
        mpo = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, mpo, control, target, param)
    end

    """
    CY gate
    """
    function apply!(gate::CY, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        _arr = [1 0 0 0; 0 1 0 0; 0 0 0 -1im; 0 0 1im 0]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    """
    SWAP gate
    """
    function apply!(gate::SWAP, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        _arr = [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    """
    CRX gate
    """
    function apply!(gate::CRX, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        ang = gate._theta * 0.5
        angpi = ang / pi
        a = cospi(angpi)
        b = -1im*sinpi(angpi)
        _arr = [1 0 0 0; 0 1 0 0; 0 0 a b; 0 0 b a]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    """
    CRY gate
    """
    function apply!(gate::CRY, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        ang = gate._theta * 0.5
        angpi = ang / pi
        a = cospi(angpi)
        b = sinpi(angpi)
        _arr = [1 0 0 0; 0 1 0 0; 0 0 a -b; 0 0 b a]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    """
    CRZ gate
    """
    function apply!(gate::CRZ, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        ang = gate._theta * 0.5
        _arr = [1 0 0 0; 0 1 0 0; 0 0 _exp(-1im*ang) 0; 0 0 0 _exp(1im*ang)]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    """
    CP gate
    """
    function apply!(gate::CP, mps::Array{Any, 1}, param::Number)
        control = gate._control
        target = gate._target
        ang = gate._theta
        _arr = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 _exp(1im*ang)]
        arr = convert(Array{ComplexF64,2}, _arr)
        mps_ = TwoQubitGate(mps, arr, control, target, param)
    end

    function _exp(theta::ComplexF64)
        thetapi = theta / pi / im
        if thetapi == conj(thetapi)
            thetapi = convert(Float64, thetapi)
        else
            error("Invalid gate parameters.")
        end
        res = cospi(thetapi) + im * sinpi(thetapi)
        return res
    end

    function gate_parser(gate::T) where T <: TNS_TwoQubitGate
        c = gate._control
        t = gate._target
        swap_idx = []
        if t > c + 1
            for i in 1:(t - c - 1)
                push!(swap_idx, (c+i-1, c+i))
            end
            gate._control = t - 1
        elseif c > t + 1
            for i in 1:(c - t - 1)
                push!(swap_idx, (t+i-1, t+i))
            end
            gate._target = c - 1
        end
        parsed_gates = []
        for swap in swap_idx
            push!(parsed_gates, SWAP(swap[1], swap[2]))
        end
        push!(parsed_gates, gate)
        for swap in reverse(swap_idx)
            push!(parsed_gates, SWAP(swap[1], swap[2]))
        end
        return parsed_gates
    end

    function gate_parser(gate::T) where T <: TNS_OneQubitGate
        return [gate]
    end

end

max_stars_repo_name: ryuNagai/blueqat-julia

max_issues_repo_name: ryuNagai/blueqat-julia

max_forks_repo_name: ryuNagai/blueqat-julia


text:
# this function shows the behavior of the
# likelihood, the log-likelihood and the average
# log-likelihood as the sample size increases
using Distributions
function main(n)
p = 0.3
y = rand(Bernoulli(p),n)
f =  pdf(Bernoulli(p), y)
@show ℒ = prod(f)
@show lnℒ = sum(log.(f))
@show avglnℒ = (1/n)*sum(log.(f))
nothing
end

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using PyPlot
using Distributions
using Random
using LinearAlgebra
using Statistics
using Distributed

struct StaticWorld{A2<:Array{Float64,2}, F<:Float64, I<:Int64, A1<:Array{Float64,1}, AA<:Array{Array{Float64,2},1}, FN<:Function}
    #Sim Priors [position, object distance, object motion, velocity]
    A::A2
    C::A2
    muPrior::A2
    endT::F
    dt::F
    endI::I
    allT::A1
    allAs::AA
    Zs::FN
end
StaticWorld(A::Array{Float64,2}, C::Array{Float64,2}, muPrior::Array{Float64,2}, endT::Float64, dt::Float64) =
    StaticWorld(A, C, muPrior, endT, dt, Integer(ceil(endT/dt)))
StaticWorld(A::Array{Float64,2}, C::Array{Float64,2}, muPrior::Array{Float64,2}) =
    StaticWorld(A, C, muPrior, 500., .5, 1000)
StaticWorld(A::Array{Float64,2}, C::Array{Float64,2}, muPrior::Array{Float64,2}, endT::Float64, dt::Float64, endI::Int) =
    StaticWorld(A, C, muPrior, endT, dt, endI,
                Array{Float64,1}(0:dt:endI*dt),
                [A^t for t in 0:endI], (A::Array{Float64,2},Z::Array{Float64,2})->A*Z)

struct SimOpts{A<:Array{Float64,1}, F<:Float64, I<:Int64}
    sigmas::A
    a::F
    N::I
end
SimOpts(num::Int64) = SimOpts(ones(num,1))
SimOpts(static::StaticWorld) = SimOpts(ones(size(static.A,1)))
SimOpts(sigmas::Array{Float64,1}) = SimOpts(sigmas, 1., 500)
SimOpts(sigmas::Array{Float64,1}, N::Int) = SimOpts(sigmas, 1., N)
SimOpts(sigmas::Array{Float64,1}, a::Float64) = SimOpts(sigmas, a, 500)

struct PlotOpts{S<:String, A<:Array{Float64,1}, F<:Float64}
    label::S
    color::A
    width::F
    alpha::F
end
PlotOpts(label::String) = PlotOpts(label, [0,0,1], 2., .2)
PlotOpts(label::String, color::Array{Float64,1}) = PlotOpts(label, color, 2., .2)

function getKalman(static::StaticWorld, initVar::Array{Float64,2}, sigmaIn::Array{Float64,2})
    Vars = Array{Array{Float64,2}, 1}(undef, static.endI+1)
    Ks = Array{Array{Float64,2}, 1}(undef, static.endI+1)
    Vars[1] = initVar
    #Get filter
    for i = 2:static.endI+1
        Ks[i] = (static.A*Vars[i-1]*static.A'*static.C')/(static.C*static.A*Vars[i-1]*static.A'*static.C' .+ sigmaIn)
        Vars[i] = (I - Ks[i]*static.C)*static.A*Vars[i-1]*static.A'
    end
    return Ks, Vars
end

struct FullWorld{S<:StaticWorld, F<:Float64, A<:Array{Float64,2}, AA<:Array{Array{Float64,2},1}, FN<:Function, FN2<:Function}
    static::S
    a::F
    sigmaIn::A
    Ys::FN
    Ks::AA
    Vars::AA
    Z::FN2
end
FullWorld(static::StaticWorld, simOpts::SimOpts) =
    FullWorld(static::StaticWorld, simOpts.a,
    diagm(simOpts.a*simOpts.sigmas), (static.C*diagm(simOpts.sigmas)*static.C')/static.dt)
FullWorld(static::StaticWorld, a::Float64, initVar::Array{Float64,2}, sigmaIn::Array{Float64,2}) =
    FullWorld(static, a, sigmaIn,
    (Z::Array{Float64,2})->rand!(MvNormal(zeros(size(sigmaIn,1)), sigmaIn), similar(static.C*static.muPrior)) + static.C*Z,
    getKalman(static, initVar, sigmaIn)...,
    ()->(static.muPrior + sqrt.(diag(initVar)).*randn(size(static.muPrior))))

function runSim(kworld::FullWorld)
    # Zs = [A*Z for A in kworld.static.allAs]
    #Zs = A->A*Z
    # Z  = kworld.Z()
    # Zs = kworld.static.Zs(Z)
    Z = collect(Iterators.repeated(kworld.Z(),kworld.static.endI+1));
    Zs = kworld.static.Zs.(kworld.static.allAs, Z);
    Ys = kworld.Ys.(Zs);#[kworld.static.C*z .+ rand!(kworld.noise, similar(kworld.static.C*Z[1])) for z in Zs]

    Mus = Array{Array{Float64,2},1}(undef, kworld.static.endI+1)
    Mus[1] = kworld.static.muPrior
    for i = 2:kworld.static.endI+1
        Mus[i] = kworld.static.A*Mus[i-1] + kworld.Ks[i]*(Ys[i] - kworld.static.C*kworld.static.A*Mus[i-1]);
    end
    return hcat(se.(Mus-Zs)...)
end

function se(M::Array{Float64,2})
    return dropdims(mean(M.^2,dims=2), dims=2)#SQRT CHANGES METRIC FROM VARIANCE (JAN MEETING 2/3/20), note:summing across dimensions sums together variances (changed to averaging)
end

function plotshade(y::Array{Float64,1}, err::Array{Float64,1}, x::Array{Float64,1}, opts::PlotOpts)
    f = fill_between(x[1:length(y)], y+err, y-err,color=opts.color, alpha=opts.alpha, linestyle="-.", hatch="/")
    p = plot(x[1:length(y)],y,color=opts.color[:],linewidth = opts.width, label=opts.label) ## change color | linewidth to adjust mean line()
end

function plotMSE(MSE::Array{Float64,2}, mVars::Array{Float64,2}, allT::Array{Float64,1}, opts::PlotOpts)

    subplot(2, 2, 1)
    plot(allT, MSE[1,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
    plot(allT, mVars[1,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
    # plotshade(RMSE[1,:], eVars[1,:], allT, opts)
    # fill_between(allT, RMSE[1,:]+mVars[1,:], RMSE[1,:]-mVars[1,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
    xlabel("time")
    ylabel("mean square error")
    title("Position MSE")

    if size(MSE,1) > 2
        subplot(2, 2, 2)
        plot(allT, MSE[2,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
        plot(allT, mVars[2,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
        # plotshade(RMSE[2,:], eVars[2,:], allT, opts)
        # fill_between(allT, RMSE[2,:]+mVars[2,:], RMSE[2,:]-mVars[2,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
        xlabel("time")
        ylabel("mean square error")
        title("Object Distance MSE")
    end

    subplot(2, 2, 3)
    plot(allT, MSE[4,:], color=opts.color[:], linewidth=opts.width, label=opts.label)
    plot(allT, mVars[4,:], color=opts.color[:], linewidth=opts.width, linestyle=":")
    # plotshade(RMSE[end,:], eVars[end,:], allT, opts)
    # fill_between(allT, RMSE[end,:]+mVars[end,:], RMSE[end,:]-mVars[end,:], color=opts.color, alpha=opts.alpha, linestyle=":", hatch="|")
    xlabel("time")
    ylabel("mean square error")
    title("Velocity MSE")

    return
end

function runBatchSim(plotOpts::PlotOpts, static::StaticWorld, simOpts::SimOpts)
    kworld = FullWorld(static, simOpts)

    ses = Array{Array{Float64,2},1}(undef, simOpts.N)
    # pmap(n->ses[n] = runSim(kworld), 1:simOpts.N)
    @simd for n = 1:simOpts.N
       ses[n] = runSim(kworld)
    end
    MSE = mean(ses)
    # eVars = var(rses; mean=RMSE)
    mVars = hcat([diag(M) for M in kworld.Vars]...)#SQRT CHANGES METRIC FROM VARIANCE (JAN MEETING 2/3/20), note:summing across dimensions sums together variances
    plotMSE(MSE, mVars, kworld.static.allT, plotOpts)

    return ses
end

max_stars_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_issues_repo_name: rhoadesScholar/Drugowitsch-Rotation

max_forks_repo_name: rhoadesScholar/Drugowitsch-Rotation


text:
"""
    to_distribution(x)

Convert to a `Uniform` distribution from an `Interval`.
"""
to_distribution(x::Interval) = Uniform(x.lo, x.hi)
to_distribution(x::Distribution) = x
to_distribution(x::ReachInterval) = to_distribution(x.dat)
to_distribution(x) = x

"""
    to_interval(x)

Convert to a `Interval` from a `Distribution`.
"""
to_interval(x::Interval) = x
to_interval(x::Distribution) = Interval(_minimum(x), _maximum(x))
to_interval(x::ReachInterval) = x.dat
to_interval(x) = x

"""
    to_reach_interval(x)

Convert to a `ReachInterval` from an `Interval` or `Distribution`.
"""
to_reach_interval(x::Interval) = ReachInterval(x)
to_reach_interval(x::Distribution) = to_reach_interval(to_interval(x))
to_reach_interval(x::ReachInterval) = x
to_reach_interval(x) = RA.Singleton([x])
to_reach_interval(x::Integer) = to_reach_interval(float(x))

"""
    recursive_convert(x)

Recursively convert a variable `x` by applying a conversion function `f`.
Mostly to be used  with `to_distribution` and `to_interval`.
"""
recursive_convert(f, x) = f(x)
recursive_convert(f, x::NamedTuple) = map(x->recursive_convert(f, x), x)
recursive_convert(f, x::AbstractArray) = recursive_convert.(f, x)

# Internal maximum and minimum to work across all uncertainty types
_minimum(x) = minimum(x)
_minimum(x::Interval) = inf(x)
_minimum(x::ReachInterval) = min(x)

_maximum(x) = maximum(x)
_maximum(x::Interval) = sup(x)
_maximum(x::ReachInterval) = max(x)


text:
"""
    get_3D_field(product_year, tracer, period, resolution, field)

Downloads and returns the 3D field of a tracer from the World Ocean Atlas.

Availabe product years are 2009, 2013, and 2018.
(2018 more likely to work.)
Tracers are either "phosphate", "nitrate", or "silicate".
Resolution is either "5°", "1°", or "0.25°".
(1° most likely to work.)
Fields are "mean" or "objectively analyzed climatology".
(mean most likely to work.)
Note WOA's nomenclature should work too, e.g.,
"p" for "phosphate", "mn" for "mean", and so on.
"""
function get_3D_field(product_year, tracer, period, resolution, field)
    println("Getting the 3D field of WOA$(my_product_year(product_year)) $(my_averaging_period(period)) $(WOA_path_varname(tracer)) $(surface_grid_size(resolution)) data")
    ds = WOA_Dataset(product_year, tracer, period, resolution)
    field3D = ds[WOA_varname(tracer, field)][:][:, :, :, 1]
    println("  Rearranging data")
    field3D = permutedims(field3D, [2 1 3])
    return field3D
end

function get_gridded_3D_field(ds, tracer, field)
    println("  Reading NetCDF file")
    field3D = ds[WOA_varname(tracer, field)][:][:, :, :, 1]
    lon   = ds["lon"][:]   .|> Float64
    lat   = ds["lat"][:]   .|> Float64
    depth = ds["depth"][:] .|> Float64
    println("  Rearranging data")
    # Reorder the variable index order (lat <-> lon from WOA to OCIM)
    field3D = permutedims(field3D, [2 1 3])
    # Rearrange longitude range (from WOA data, which is -180:180, to 0:360)
    lon = mod.(lon, 360)
    lon_reordering = sortperm(lon)
    lon = lon[lon_reordering]
    field3D .= field3D[:, lon_reordering, :]
    return field3D, lat, lon, depth
end

"""
    mean_std_and_number_obs(ds, tracer)

Returns a DataFrame containing the following columns
lat, lon, depth, mean, std, nobs
"""
function mean_std_and_number_obs(ds, tracer)
    println("  Reading NetCDF file")
    lon   = mod.(ds["lon"][:] .|> Float64, 360)
    lat   = ds["lat"][:]   .|> Float64
    depth = ds["depth"][:] .|> Float64
    μ3D    = ds[WOA_varname(tracer, "mn")][:][:, :, :, 1]
    σ3D    = ds[WOA_varname(tracer, "sd")][:][:, :, :, 1]
    nobs3D = ds[WOA_varname(tracer, "dd")][:][:, :, :, 1]
    println("  Filtering missing data")
    CI = findall(@. !ismissing(μ3D) & !ismissing(nobs3D) & !iszero(nobs3D))
    lon1D   = lon[map(x -> x.I[1], CI)]
    lat1D   = lat[map(x -> x.I[2], CI)]
    depth1D = depth[map(x -> x.I[3], CI)]
    μ1D     = μ3D[CI] .|> Float64
    σ1D     = σ3D[CI] .|> Float64
    nobs1D  = nobs3D[CI] .|> Int64
    return DataFrame(
                     :Latitude => lat1D,
                     :Longitude => lon1D,
                     :Depth => depth1D,
                     :Mean => μ1D,
                     :Std => σ1D,
                     :n_obs => nobs1D,
                    )
end

function filter_gridded_3D_field(field3D, lat, lon, depth)
    # Find where there is data for both mean and std
    println("  Filtering data")
    CI = findall(x -> !ismissing(x) && !iszero(x), field3D) # filter out fill-values and 0's
    fieldvec = field3D[CI]
    latvec = lat[map(x -> x.I[1], CI)]
    lonvec = lon[map(x -> x.I[2], CI)]
    depthvec = depth[map(x -> x.I[3], CI)]
    return fieldvec, latvec, lonvec, depthvec, CI
end

get_unit(ds, tracer, field) = convert_to_Unitful(ds[WOA_varname(tracer, field)].attrib["units"])

function mean_and_variance_gridded_3d_field(grid::OceanGrid, field3D, lat, lon, depth)
    fieldvec, latvec, lonvec, depthvec, CI = filter_gridded_3D_field(field3D, lat, lon, depth)
    println("  Averaging data over each grid box")
    χ_3D = zeros(size(grid))
    σ²_3D = zeros(size(grid))
    n_3D = zeros(size(grid))
    # Use NearestNeighbors to bin into OceanGrid
    gridbox_centers = [ustrip.(vec(grid.lat_3D)) ustrip.(vec(grid.lon_3D)) ustrip.(vec(grid.depth_3D))] # TODO Maybe add option for using iwet instead of full vec
    gridbox_centers = permutedims(gridbox_centers, [2, 1])
    kdtree = KDTree(gridbox_centers)
    for i in eachindex(CI)
        idx = knn(kdtree, [lonvec[i], latvec[i], depthvec[i]], 1, true)[1][1]
        χ_3D[idx] += fieldvec[i]      # μ = Σᵢ μᵢ / n
        σ²_3D[idx] += fieldvec[i]^2   # σ² = Σᵢ μᵢ² / n - μ²
        n_3D[idx] += 1
    end
    χ_3D .= χ_3D ./ n_3D               # μ = Σᵢ μᵢ / n
    σ²_3D .= σ²_3D ./ n_3D .- χ_3D.^2  # σ² = Σᵢ μᵢ² / n - μ²
    # Enforce μ = 0 and σ = ∞ where no observations
    # (Instead of NaNs)
    println("  Setting μ = 0 and σ² = ∞ where no obs")
    χ_3D[findall(n_3D .== 0)] .= 0.0
    σ²_3D[findall(n_3D .== 0)] .= Inf
    println("  Setting a realistic minimum for σ²")
    meanχ = mean(χ_3D, weights(n_3D))
    σ²_3D .= max.(σ²_3D, 1e-4meanχ^2)
    return χ_3D, σ²_3D
end

function convert_to_SI_unit!(χ_3D, σ²_3D, ds, tracer, field)
    # Convert to SI units
    χ_unit = get_unit(ds, tracer, field)
    χ_3D .*= ustrip(upreferred(1.0χ_unit))
    σ²_3D .*= ustrip(upreferred(1.0χ_unit^2))
end

function fit_to_grid(grid::OceanGrid, product_year, tracer, period, resolution, field)
    ds = WOA_Dataset(product_year, tracer, period, resolution)
    field3D, lat, lon, depth = get_gridded_3D_field(ds, tracer, field)
    χ_3D, σ²_3D = mean_and_variance_gridded_3d_field(grid, field3D, lat, lon, depth)
    convert_to_SI_unit!(χ_3D, σ²_3D, ds, tracer, field)
    return χ_3D, σ²_3D
end

#==================================
Helper functions
==================================#

function WOA_Dataset(product_year, tracer, period, resolution; verbose=false)
    println("Registering WOA data with DataDeps")
    verbose && print(": ", url_WOA_THREDDS(product_year, tracer, period, resolution))
    register_WOA(product_year, tracer, period, resolution)
    ddstr = @datadep_str string(my_DataDeps_name(product_year, tracer, period, resolution),
                        "/",
                        WOA_NetCDF_filename(product_year, tracer, period, resolution))
    Dataset(ddstr)
end

function remove_DataDep(product_year, tracer, period, resolution)
    println("Removing WOA$(my_product_year(product_year)) $(my_averaging_period(period)) $(WOA_path_varname(tracer)) $(surface_grid_size(resolution)) data")
    nc_file = @datadep_str string(my_DataDeps_name(product_year, tracer, period, resolution), "/", WOA_NetCDF_filename(product_year, tracer, period, resolution))
    rm(nc_file; recursive=true, force=true)
end


text:
#=
lgn_l6l4:
- Julia version: 1.4.0
- Author: niallcullinane
- Date: 2020-06-07
=#