text:
"""A convenient macro copied from Mocha.jl that could be used to define structs
with default values and type checks. For example
```julia
@defstruct MyStruct Any (
  field1 :: Int = 0,
  (field2 :: AbstractString = "", !isempty(field2))
)
```
The macro will define a constructor without keyword arguments and `set` function.
We can use the `set` function to assign a field's value and check if the value is
valid.
"""
macro defstruct(name, fields)
  @assert(isa(name, Symbol) || isa(name, Expr),
  "name should be a Symbol")
  if isa(name, Expr)
    @assert(name.head == :(<:) && length(name.args) == 2 && all([isa(m) for m in name.args]),
    "name should be of form 'Name <: SuperType'")
    name, super_type = name.args
  else
    super_type = :Any
  end

  @assert(isa(fields, Expr) && fields.head == :tuple)

  field_names       = Vector{Symbol}(0)
  field_defs        = Vector{Expr}(0)
  field_defaults    = Vector{Expr}(0)
  field_asserts     = Vector{Pair{String, Expr}}(0)

  for field in fields.args
    if field.head == :tuple
      @assert(length(field.args) == 2, "length of tuple should be 2")
      push!(field_asserts, Pair(string(field.args[1].args[1].args[1]), field.args[2]))
      # push!(field_asserts, esc(field.args[2]))
      field = field.args[1]
    end
    @assert(field.head == :(=), string(field) * "should have default value")
    push!(field_defaults, esc(field))
    field = field.args[1]
    push!(field_defs, field)
    push!(field_names, field.args[1])
  end

  type_body = Expr(:block, field_defs...)
  struct_body = Expr(:block, field_defaults..., Expr(:call, :new, field_names...))
  asserts   = Expr(:call, :(Dict{String, Expr}), field_asserts...)
  set_ = esc(Symbol("set_" * string(name)))
  set = esc(Symbol("set" * string(name)))
  inner_func = Expr(:call, esc(:tmp), Expr(:parameters, Expr(:kw, :object, nothing), Expr(:kw, :(param::String), nothing), Expr(:kw, :(value::Any), nothing)))
  # outer_func = Expr(:(=), set, Expr(:call, set_))
  outer_func = "set(object::$(string(name)), param::String, value::Any) = set_" * string(name) * "()(object = object, param = param, value = value)"
  quote
    type $(name) <: $(super_type)
      $type_body

      function $(name)()
        $struct_body
      end
    end

    function $(set_)()
      $(esc(:asserts)) = $asserts
      $(inner_func) = begin
        if haskey(asserts, param)
          eval(Expr(:(=), Symbol(param), value))
          @assert(eval(asserts[param]), param * " is invalid")
        end;
        setfield!(object, Symbol(param), value)
      end
      return $(esc(:tmp))
    end

    # $outer_func
    eval(parse($(esc(outer_func))))
  end
end


text:
export Jacobi

"""
Jacobi Distribution

If E ~ Wishartp(I, n1) and H ~ Wishartp(I, n2) are independent,
then E(E + H)-1 has Jacobi(n1, n2, p) distribution.

If λi are the eigenvalues of EH-1 and μi are the Jacobi eigenvalues,
then μi = λi/(1 + λi) and  λi = μi/(1 - μi)
"""
struct Jacobi <: ContinuousMatrixDistribution
    beta::Integer
    n1::Integer
    n2::Integer
    p::Integer
end

function randreduced(d::Jacobi)
    randtridiagonal(d)
end

function randtridiagonal(d::Jacobi)

    a = d.beta/2 * (d.n1 - d.p + 1) - 1;
    b = d.beta/2 * (d.n2 - d.p + 1) - 1;

    alphaeven = [2*rand(Beta((2*d.p - k - 2)*d.beta/4 + a + 1, (2*d.p - k - 2) * d.beta/4 + b + 1)) - 1 for k in 0:(2 * d.p-2) if k%2 == 0]
    alphaodd = [2*rand(Beta((2*d.p - k - 1)*d.beta/4, (2*d.p - k - 3) * d.beta/4 + a + b + 2)) - 1 for k in 0:(2 * d.p-2) if k%2 == 1]

    alphaevenleft = [alphaeven; 0]
    alphaevenright = [0; alphaeven]
    alphaodd = [-1; alphaodd; -1]

    dv = (1 .- alphaodd) .* alphaevenleft - (1 .+ alphaodd).* alphaevenright
    dv = dv[1:(end-1)]

    ev = sqrt.((1 .- alphaodd[1:(end-1)]) .* (1 .- alphaeven.^2) .* (1 .+ alphaodd[2:end]))
    ev = ev[1:end-1]

    (SymTridiagonal(dv, ev) + 2I)/4
end

ext: jl


text:
using TextUnidecode
using Test

@testset "TextUnidecode.jl" begin
    @test unidecode("u") == "u"
    @test unidecode("uua") == "uua"
    @test unidecode(string(Char(0x10000))) == ""
    @test unidecode("ü") == "u"
    @test unidecode("ø") == "o"
    @test unidecode("😜") == ""
    @test unidecode("Ｈ") == "H"
    @test unidecode("南无阿弥陀佛") == "Nan Wu A Mi Tuo Fo"
end


text:
include("../src/lut/mc.jl")

connectivity_table = []
for i in eachindex(tri_table)
    a = zeros(Int, 15)
    for j = 1:15
        elt = findfirst(x -> x == tri_table[i][j], [_mc_verts[i]...])
        if elt == nothing || _mc_verts[i][elt] == 0
            a[j] = 0
        else
            a[j] = elt
        end
    end
    push!(connectivity_table, tuple(a...))
end

for elt in connectivity_table
    println(string(elt)*",")
end


text:
# This file is a part of SimilaritySearch.jl

using LinearAlgebra

function create_sequence(dim, sort)
    s = rand(1:10, dim)
    if sort
        sort!(s)
        s = unique(s)
    end

    s
end

include("testdb.jl")
include("testresults.jl")
include("testseq.jl")
include("testallknn.jl")
include("testneardup.jl")
include("testsearchgraph.jl")

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
# Working with the API
"""
    createtable!(name)
Creates a table in the database based on the O*NET API specification.
"""
function createtable!(name::AbstractString)
    url = "https://services.onetcenter.org/ws/database/info/$name"
    request = HTTP.get(url, headers)
    @assert request.status == 200
    schema = request.body |> String |> JSON3.read
    tblname = schema.table_id
    datatypes = mapreduce(varschema,
                          (x, y) -> "$x, $y",
                          schema.column)
    vals = mapreduce(x -> "\$$x",
                     (x, y) -> "$x, $y",
                     eachindex(schema.column))
    conn = dbconnect()
    execute(conn, "create table if not exists onet.$tblname($datatypes)")
    close(conn)
    "insert into onet.$tblname values ($vals);"
end

"""
    onettable(table)
Given a table node, it returns a row with the table information.
"""
function onettable(table)
    name = table.table_id
    request = HTTP.get(table.info, headers)
    @assert request.status == 200
    json = request.body |>
        String |>
        JSON3.read
    request = HTTP.get(table.rows, headers)
    @assert request.status == 200
    total = request.body |>
        String |>
        JSON3.read |>
        (obj -> getproperty(obj, :total))
    [ (table_name = name,
       table_description = table.description,
       variable_id = elem.column_id,
       variable_name = elem.title,
       variable_description = elem.description,
       n = total)
       for elem ∈ json.column ]
end

"""
    onettables!()
Uploads a `_metadata` table to the database.
"""
function onettables!()
    request = HTTP.get("https://services.onetcenter.org/ws/database/",
                       headers)
    @assert request.status == 200
    tables = request.body |> String |> JSON3.read |>
        (obj -> getproperty(obj, :table))
    data = onettable(tables[1])
    for table ∈ tables[2:end]
        append!(data, onettable(table))
    end
    vals = mapreduce(x -> "\$$x",
                     (x, y) -> "$x, $y",
                     eachindex(propertynames(output[1])))
    conn = dbconnect()
    execute(conn, """create table if not exists onet.onet_tables (
                     table_name text not null,
                     table_description text not null,
                     variable_id text not null,
                     variable_name text not null,
                     variable_description text not null,
                     n int not null
                     )
                  """)
    statement = "insert into onet.onet_tables values ($vals)"
    execute(conn, "begin;")
    load!(data, conn, statement)
    execute(conn, "commit;")
    close(conn)
end

"""
    processtable!(name::AbstractString)
Adds the table data to the database.
"""
function processtable!(name::AbstractString)
    url = "https://services.onetcenter.org/ws/database/rows/$name"
    request = HTTP.get(url, headers)
    @assert request.status == 200
    json = request.body |>
        String |>
        JSON3.read
    url = replace(json.link[1].href,
                  r"(?<=start=)\d+&end=\d+$" => "1&end=$(json.total)")
    request = HTTP.get(url, headers)
    @assert request.status == 200
    json = request.body |>
        String |>
        JSON3.read
    statement = createtable!(name)
    data = rowtable(json.row) |>
        (tbl -> transform(tbl,
                          Dict(cn => x -> something(x, missing)
                               for cn ∈ eachindex(tbl[1])))) |>
        rowtable
    isdir("/mnt/volume_nyc1_02/onet") || mkdir("/mnt/volume_nyc1_02/onet")
    CSV.write("data.csv", data)
    conn = dbconnect()
    execute(conn, "COPY $name from '/mnt/volume_nyc1_02/onet/data.csv' with (format csv, header)")
    close(conn)
    isfile("/mnt/volume_nyc1_02/onet/data.csv") &&
        rm("/mnt/volume_nyc1_02/onet/data.csv")
    true
end

"""
    onet!()
Runs the application to refresh all the data in the database.
"""
function onet!()
    conn = dbconnect()
    tbls = execute(conn, "select distinct table_name, n from onet.onet_tables") |>
        (tbl -> [(name = elem.table_name, count = elem.n) for elem ∈ tbl]) |>
        columntable
    close(conn)
    foreach(drop_duplicates!, tbls.name)
    to_add = filter(!isdone, tbls.name) |> sort!
    x = rowtable(tbls) |>
        (tbl -> filter(row -> row.name ∈ to_add, tbl)) |>
        (tbl -> sort!(tbl, by = last))
    foreach(processtable!, to_add)
    true
end


text:
"""
    geometric_factors(x, y, Dr, Ds)
    geometric_factors(x, y, z, Dr, Ds, Dt)

Compute metrics of mappings between "real" elements and reference elements,
outward pointing normals on faces of every elements, and Jacobian.

x,y,z are arrays of coordinates, and Dr, Ds, Dt are nodal differentiation matrices

Geometric terms in 3D are constructed to ensure satisfaction of free-stream
preservation using the curl-based construction of David Kopriva (2001).

"""

# 2D version
function geometric_factors(x, y, Dr, Ds)
    "Transformation and Jacobian"

    xr = Dr*x;   xs = Ds*x;
    yr = Dr*y;   ys = Ds*y;

    J = -xs.*yr + xr.*ys;
    rxJ =  ys;  sxJ = -yr;
    ryJ = -xs;  syJ =  xr;

    return rxJ, sxJ, ryJ, syJ, J
end

# 3D version. Filters = tuple of filtering matrices.
function geometric_factors(x, y, z, Dr, Ds, Dt, Filters=(I,I,I))

    xr = Dr*x;  xs = Ds*x;  xt = Dt*x
    yr = Dr*y;  ys = Ds*y;  yt = Dt*y
    zr = Dr*z;  zs = Ds*z;  zt = Dt*z

    Fr = (Dr*y).*z
    Fs = (Ds*y).*z
    Ft = (Dt*y).*z
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    rxJ = Dt*(Fs) - Ds*(Ft)
    sxJ = Dr*(Ft) - Dt*(Fr)
    txJ = Ds*(Fr) - Dr*(Fs)

    Fr = (Dr*x).*z
    Fs = (Ds*x).*z
    Ft = (Dt*x).*z
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    ryJ = -(Dt*(Fs) - Ds*(Ft))
    syJ = -(Dr*(Ft) - Dt*(Fr))
    tyJ = -(Ds*(Fr) - Dr*(Fs))

    Fr = (Dr*y).*x
    Fs = (Ds*y).*x
    Ft = (Dt*y).*x
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    rzJ = -(Dt*(Fs) - Ds*(Ft))
    szJ = -(Dr*(Ft) - Dt*(Fr))
    tzJ = -(Ds*(Fr) - Dr*(Fs))

    J = @. xr*(ys*zt-zs*yt) - yr*(xs*zt-zs*xt) + zr*(xs*yt-ys*xt)

    return rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


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

using ArraysOfArrays
using Test
using StatsBase
using Statistics

@testset "StatsBase support" begin
    VV = [rand(3) for i in 1:10]
    VV_aosa = ArrayOfSimilarArrays(VV)

    VA = [rand(2,3,3) for i in 1:10]
    VA_aosa = ArrayOfSimilarArrays(VA)

    w = FrequencyWeights(rand(10))

    array_cmp(A, B) = (A ≈ B) && (size(A) == size(B))

    # sum and mean for Vector{Vector} with weights currently fail with
    # the implementations in StatsBase. This should be considered a
    # bug in StatsBase, since Base and Statistics support sum and mean
    # for Vector{Vector} without weights. Also, adding products of vectors
    # and weights is perfectly natural, mathematically.

    _sum(A::AbstractVector{<:AbstractArray}, w::AbstractWeights) =
        sum(A .* w)

    _mean(A::AbstractVector{<:AbstractArray}, w::AbstractWeights) =
        _sum(A, w) ./ sum(w)

    @testset "sum and mean" begin

        @test array_cmp(@inferred(sum(VV_aosa, w)), _sum(VV, w))
        @test array_cmp(@inferred(sum(VA_aosa, w)), _sum(VA, w))

        @test array_cmp(@inferred(mean(VV_aosa, w)), _mean(VV, w))
        @test array_cmp(@inferred(mean(VA_aosa, w)), _mean(VA, w))
    end

    # Weighted var and std are currently not supported for Vector{Vector} by
    # StatsBase. This should be considered a bug in StatsBase, since
    # unweighted var and std for Vector{Vector} are supported by Statistics.

    function _var(A::AbstractVector{<:AbstractArray}, w::FrequencyWeights; corrected = true)
        wmean_A = _mean(A, w)
        wsum = sum(w)
        wsum_corr = corrected ? -1 : 0
        sum([(x .- wmean_A).^2 for x in A] .* w) ./ (wsum + wsum_corr)
    end

    _std(A::AbstractVector{<:AbstractArray}, w::AbstractWeights; corrected = true) =
        sqrt.(_var(A, w, corrected = corrected))

    @testset "var and std" begin
        @test array_cmp(@inferred(var(VV_aosa, w)), _var(VV_aosa, w))
        @test array_cmp(@inferred(var(VV_aosa, w, corrected = false)), _var(VV_aosa, w, corrected = false))
        @test array_cmp(@inferred(var(VA_aosa, w)), _var(VA_aosa, w))
        @test array_cmp(@inferred(var(VA_aosa, w, corrected = false)), _var(VA_aosa, w, corrected = false))

        @test array_cmp(@inferred(std(VV_aosa, w)), _std(VV_aosa, w))
        @test array_cmp(@inferred(std(VV_aosa, w, corrected = false)), _std(VV_aosa, w, corrected = false))
        @test array_cmp(@inferred(std(VA_aosa, w)), _std(VA_aosa, w))
        @test array_cmp(@inferred(std(VA_aosa, w, corrected = false)), _std(VA_aosa, w, corrected = false))
    end

    # For weighted cov of Vector{Vector}, StatsBase currently returns a vector
    # instead of a matrix, with `cov(VV, fill(1, 10)) != cov(VV)`.
    # This should be considered a bug in StatsBase.

    function _cov(A::AbstractVector{<:AbstractVector}, w::FrequencyWeights; corrected = true)
        wmean_A = _mean(A, w)
        wsum = sum(w)
        wsum_corr = corrected ? -1 : 0
        sum([[(A[i][j] - wmean_A[j]) * (A[i][k] - wmean_A[k]) * w[i] for j in eachindex(A[i]), k in eachindex(A[i])] for i in eachindex(A)]) ./ (wsum + wsum_corr)
    end

    @testset "cov" begin
        @test array_cmp(@inferred(cov(VV_aosa, w)), _cov(VV_aosa, w))
        @test array_cmp(@inferred(cov(VV_aosa, w, corrected = false)), _cov(VV_aosa, w, corrected = false))
    end

    # Weighted cor is currently not supported for Vector{Vector} by StatsBase.
    # This should be considered a bug in StatsBase, since unweighted cor
    # for Vector{Vector} is supported by Statistics.

    _cor(A::AbstractVector{<:AbstractVector}, w::AbstractWeights) = cov2cor(_cov(A, w), _std(A, w))

    @testset "cor" begin
        @test array_cmp(@inferred(cor(VV_aosa, w)), _cor(VV, w))
    end
end


text:
using CSV, DataFrames, InfluxDB2, Dates, Test, HTTP

# Lineprotocol
@testset "lineprotocol" begin
    testtable = CSV.File(joinpath(@__DIR__, "testtable.csv"))

    iobuffer = InfluxDB2.table2lineprotocol("testmeasurement", testtable; precision = :ms)
    correctlinep = "testmeasurement,c=hello b=2.0 1637496010000\ntestmeasurement,c=world b=4.0 1637496070000\ntestmeasurement,c=hello b=6.0 1637496130000\n"
    @test String(take!(iobuffer)) == correctlinep
end

if isfile(joinpath(@__DIR__, "influxsettings.txt"))
    token, org, bucket = readlines(joinpath(@__DIR__, "influxsettings.txt"))
    influx = InfluxServer("http://localhost:8086", org, token)
    measurementname = "juliatest"
    ndata = 10

    @testset "writing" begin
        # write without compression
        testtable = DataFrame(Dict("timestamp"=>[now(UTC)-Second(ndata)+Second(i) for i = 1:ndata], "f_Floatfield"=>randn(ndata), "f_Intfield"=>rand(1:10, ndata), "f_Stringfield"=>[rand(("String1", "String2", "String3")) for _ =1:ndata], "f_Boolfield"=>rand((false, true), ndata), "t_tag1"=>[rand(["hello", "world"]) for _ =1:ndata]))
        writeresp = writetable(influx, bucket, measurementname, testtable)
        @test writeresp.status == 204

        # write using gzip compression
        testtable2 = DataFrame(Dict("timestamp"=>[now(UTC)-Second(ndata*2)+Second(i) for i = 1:ndata], "f_Floatfield"=>randn(ndata), "f_Intfield"=>rand(1:10, ndata), "f_Stringfield"=>[rand(("String1", "String2", "String3")) for _ =1:ndata], "f_Boolfield"=>rand((false, true), ndata), "t_tag1"=>[rand(["hello", "world"]) for _ =1:ndata]))
        writeresp = writetable(influx, bucket, measurementname, testtable2; compression = :gzip)
        @test writeresp.status == 204
    end

    @testset "reading" begin
        # Query without compression
        data = simplequery(influx, bucket, measurementname, ["Floatfield", "Intfield", "Stringfield", "Boolfield"], (now(UTC)-Day(2), now(UTC)); tags = ["tag1"=>"hello", "tag2" =>"hi"])
        @test length(data) == 4 # we expect 4 tables, one for each field

        # Query foor gzip compressed data
        data = simplequery(influx, bucket, measurementname, ["Floatfield", "Intfield", "Stringfield", "Boolfield"], (now(UTC)-Day(2), now(UTC)); compression = :gzip)
        @test length(data) == 4
        @test size(data[1], 1) >= ndata

        querystring = "from(bucket: \"$bucket\")
        |>range(start: $(Dates.format(now(UTC)-Day(2), "yyyy-mm-ddTHH:MM:SS.sZ")), stop: $(Dates.format(now(UTC), "yyyy-mm-ddTHH:MM:SS.sZ")))
        |>filter(fn: (r) => r[\"_measurement\"] == \"$measurementname\")
        |>filter(fn: (r) =>r[\"_field\"] == \"Floatfield\" or r[\"_field\"] == \"Intfield\")
        |>filter(fn: (r) =>r[\"tag1\"] == \"hello\")
        |>drop(columns: [\"_measurement\", \"_start\", \"_stop\"])
        |>group(columns: [\"_field\"], mode: \"by\")
        |>sort(columns: [\"_time\"])"

        data = fluxquery(influx, querystring)
        @test length(data) == 2
    end
end


text:
using Documenter, RCall

makedocs(
    modules = [RCall],
    clean   = false,
    format   = Documenter.Formats.HTML,
    sitename = "RCall.jl",
    pages    = Any[
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Getting Started" => "gettingstarted.md",
        "Internal" => "internal.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaStats/RCall.jl.git",
    julia  = "0.4",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
# A basic functionality
include("basics.jl")

# An application for shared ballot
# include("synchronicballot.jl")


text:
#! /usr/bin/env julia
using ZenFab

if !isempty(ARGS) && ARGS[1] == "serve"
    cd(Pkg.dir("ZenFab", "src"))
    ZenFab.escher_serve()
end


text:
module Data

using CSV
using DataFrames
using JSON

export atom_data_df, feature_info, valenceshell_conf_df

atom_data_path = joinpath(@__DIR__, "..", "data", "pymatgen_atom_data.csv")
atom_data_df = DataFrame(CSV.File(atom_data_path))

feature_info_path = joinpath(@__DIR__, "..", "data", "feature_info.json")
feature_info = JSON.parsefile(feature_info_path)

valenceshell_conf_df = select(
    atom_data_df,
    :Symbol,
    "Electronic Structure" =>
        (
            x -> replace.(x, r"\[(.+)\]\.(?<valence>\w+)" => s"\g<valence>")
        ) => "Electronic Structure",
)

end


text:
# Perceived noise level, tone-corrected
function f_noy(settings::PyObject, spl)

    # Number of time steps
    n_t = size(spl)[1]

    # Get type of input vector
    T = eltype(spl)

    # Compute noy
    spl_a = [91. , 85.9, 87.3, 79. , 79.8, 76. , 74. , 74.9, 94.6,  1e8,  1e8, 1e8,  1e8,  1e8,  1e8,  1e8,  1e8,  1e8,  1e8,  1e8,  1e8,  1e8, 44.3, 50.7]
    spl_b = [64, 60, 56, 53, 51, 48, 46, 44, 42, 40, 40, 40, 40, 40, 38, 34, 32, 30, 29, 29, 30, 31, 34, 37]
    spl_c = [52, 51, 49, 47, 46, 45, 43, 42, 41, 40, 40, 40, 40, 40, 38, 34, 32, 30, 29, 29, 30, 31, 34, 37]
    spl_d = [49, 44, 39, 34, 30, 27, 24, 21, 18, 16, 16, 16, 16, 16, 15, 12,  9, 5,  4,  5,  6, 10, 17, 21]
    spl_e = [55, 51, 46, 42, 39, 36, 33, 30, 27, 25, 25, 25, 25, 25, 23, 21, 18, 15, 14, 14, 15, 17, 23, 29]
    m_b   = [0.043478, 0.04057 , 0.036831, 0.036831, 0.035336, 0.033333, 0.033333, 0.032051, 0.030675, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.02996 , 0.02996 , 0.02996 , 0.02996 , 0.02996 , 0.02996 , 0.02996 , 0.042285, 0.042285]
    m_c   = [0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103, 0.030103,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8,      1e8, 0.02996 , 0.02996 ]
    m_d   = [0.07952 , 0.06816 , 0.06816 , 0.05964 , 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.053013, 0.05964 , 0.053013, 0.053013, 0.047712, 0.047712, 0.053013, 0.053013, 0.06816 , 0.07952 , 0.05964 ]
    m_e   = [0.058098, 0.058098, 0.052288, 0.047534, 0.043573, 0.043573, 0.040221, 0.037349, 0.034859, 0.034859, 0.034859, 0.034859, 0.034859, 0.034859, 0.034859, 0.040221, 0.037349, 0.034859, 0.034859, 0.034859, 0.034859, 0.037349, 0.037349, 0.043573]

    # Generate function N(SPL) that is always 0 for any SPL
    #N = spl
    #N[findall( spl .< spl_d)] = zeros(T, (n_t, settings.N_f))[findall( spl .< spl_d)]
    ##N[findall(spl_a .<= spl          )] =       10 .^(m_c .* (spl - spl_c))[findall(spl_a .<= spl          )]
    #N[findall(spl_b .<= spl .<= spl_a)] =       10 .^(m_b .* (spl - spl_b))[findall(spl_b .<= spl .<= spl_a)]
    #N[findall(spl_e .<= spl .<= spl_b)] = 0.3 * 10 .^(m_e .* (spl - spl_e))[findall(spl_e .<= spl .<= spl_b)]
    #N[findall(spl_d .<= spl .<= spl_e)] = 0.1 * 10 .^(m_d .* (spl - spl_d))[findall(spl_d .<= spl .<= spl_e)]

    #N = spl.^ .- 1

    N = zeros(T, (n_t, settings.N_f))
    for i in range(1, n_t, step=1)
        for j in range(1, settings.N_f, step=1)
            if spl_a[j] <= spl[i,j]
                N[i,j] = 10^(m_c[j] * (spl[i,j] - spl_c[j]))
            elseif spl_b[j] <= spl[i,j] <= spl_a[j]
                N[i,j] = 10^(m_b[j] * (spl[i,j] - spl_b[j]))
            elseif spl_e[j] <= spl[i,j] <= spl_b[j]
                N[i,j] = 0.3 * 10^(m_e[j] * (spl[i,j] - spl_e[j]))
            elseif spl_d[j] <= spl[i,j] <= spl_e[j]
                N[i,j] = 0.1 * 10^(m_d[j] * (spl[i,j] - spl_d[j]))
            else
                N[i,j] = 0.
            end
        end
    end

    #N[findall(spl_a.* ones(T, (n_t, 1)) .<= spl)] = 10 .^(m_c.* ones(T, (n_t,1)) .* (spl - spl_c.* ones(T, (n_t,1))))[findall(spl_a.* ones(T, (n_t,1)) .<= spl)]
    #N[findall(spl_b.* ones(T, (n_t, 1)) .<= spl .<= spl_a.* ones(T, (n_t,1)) )] = 10 .^(m_b.* ones(T, (n_t,1)) .* (spl - spl_b.* ones(T, (n_t,1))))[findall(spl_b.* ones(T, (n_t,1)) .<= spl .<= spl_a.* ones(T, (n_t,1)) )]
    #N[findall(spl_e.* ones(T, (n_t, 1)) .<= spl .<= spl_b.* ones(T, (n_t,1)) )] = 0.3 * 10 .^(m_e.* ones(T, (n_t,1)) .* (spl - spl_e.* ones(T, (n_t,1))))[findall(spl_e.* ones(T, (n_t,1)) .<= spl .<= spl_b.* ones(T, (n_t,1)) )]
    #N[findall(spl_d.* ones(T, (n_t, 1)) .<= spl .<= spl_e.* ones(T, (n_t,1)) )] = 0.1 * 10 .^(m_d.* ones(T, (n_t,1)) .* (spl - spl_d.* ones(T, (n_t,1))))[findall(spl_d.* ones(T, (n_t,1)) .<= spl .<= spl_e.* ones(T, (n_t,1)) )]

    return N
end

function f_pnl(N)

    # Number of time steps
    n_t = size(N)[1]

    # Compute perceived noise level
    pnl = zeros(eltype(N), n_t)
    n_max = maximum(N, dims=2)
    n = n_max .+ 0.15*(sum(N, dims=2) .- n_max)
    pnl = 40. .+ 10. / log10(2) * log10.(n)

    # Remove all PNL below 0
    pnl = clamp.(pnl, 0, Inf)

    return reshape(pnl, (n_t, ))
end

function f_tone_corrections(settings::PyObject, spl)

    # Number of time steps
    n_t = size(spl)[1]

    # Get type of input vector
    T = eltype(spl)

    # Step 1: Compute the slope of SPL
    s = zeros(T, (n_t, settings.N_f))
    s[:,4:end] = @view(spl[:,4:end]) .- @view(spl[:,3:end-1])

    # Step 2: Compute the absolute value of the slope and compare to 5
    slope = zeros(T, (n_t, settings.N_f))
    slope[:,4:end] = @view(s[:,4:end]) .- @view(s[:,3:end-1])

    slope_large = zeros(T, (n_t, settings.N_f))
    idx_slope_large = findall(abs.(slope).>5)
    slope_large[idx_slope_large] = ones(size(spl))[idx_slope_large]

    # Step 3: Compute the encircled values of SPL
    spl_large = zeros(T, (n_t, settings.N_f))
    spl_large[findall((slope_large .== 1) .* (s.> 0) .* (s .- hcat(zeros(T, (n_t, )),s[:,1:end-1]) .> 0) )] = ones(size(spl))[findall( (slope_large .== 1) .* (s.>0) .* (s .- hcat(zeros(T, (n_t, )),s[:,1:end-1]) .> 0) )]
    interm = hcat(zeros(T, (n_t,1)), spl_large)
    interm[findall( (slope_large .== 1) .* (s.<=0) .* (hcat(zeros(T, (n_t, )),s[:,1:end-1]) .> 0) )] = ones(n_t, settings.N_f)[findall( (slope_large .== 1) .* (s.<=0) .* (hcat(zeros(T, (n_t, )),s[:,1:end-1]) .> 0) )]
    spl_large = interm[:,2:end]

    # Step 4: Compute new adjusted sound pressure levels SPL'
    spl_p = zeros(T, (n_t, settings.N_f))
    spl_p[:,24] = spl[:,22] + s[:,22]
    spl_p[:,2:23]= 0.5 * (spl[:,1:22] .+ spl[:,3:24])
    spl_p[findall(spl_large .==0)] = spl[findall(spl_large .==0)]

    # Step 5: Recompute the slope s'
    s_p = zeros(T, (n_t, settings.N_f + 1))
    s_p[:,4:end-1] = @view(spl_p[:,4:end]) .- @view(spl_p[:,3:end-1])
    s_p[:,3] = @view(s_p[:,4])
    # Compute 25th imaginary band
    s_p[:,25] = @view(s_p[:,24])

    # Step 6: Compute arithmetic average of the 3 adjacent slopes
    s_bar = zeros(T, (n_t, settings.N_f))
    s_bar[:,3:end-1] = 1. / 3. * (@view(s_p[:,3:end-2]) + @view(s_p[:,4:end-1]) + @view(s_p[:,5:end]))

    # Step 7: Compute final 1/3 octave-band sound pressure level
    spl_pp = zeros(T, (n_t, settings.N_f))
    spl_pp[:, 3] = spl[:,3]
    spl_pp[:, 4:end]  = cumsum(s_bar[:, 3:end-1], dims=2) .+  spl[:,3]

    # Step 8: Compute the difference between SPL and SPL_pp
    F = zeros(T, (n_t, settings.N_f))
    F[3:end] = spl[3:end] .- spl_pp[3:end]
    F[findall(F.<1.5)] = zeros(T, (n_t, settings.N_f))[findall(F.<1.5)]

    # Step 9: Compute the correction factor C
    c0  = zeros(T, (n_t, 2))
    c10 = zeros(T, (n_t, 8))
    c10[findall(1.5 .<= F[:,3:10] .< 3 )]    = @view(F[:,3:10][findall(1.5 .<= F[:,3:10] .< 3 )]) / 3. .- 0.5
    c10[findall(3.0 .<= F[:,3:10] .< 20)]    = @view(F[:,3:10][findall(3.0 .<= F[:,3:10] .< 20)]) / 6.
    c10[findall(20. .<= F[:,3:10]      )]    = @view(F[:,3:10][findall(20. .<= F[:,3:10]      )]) * 0. .+ (3 .+ 1/3.)
    c20 = zeros(T, (n_t, 10))
    c20[findall(1.5 .<= F[:,11:20] .< 3 )] = @view(F[:,11:20][findall(1.5 .<= F[:,11:20] .< 3 )])*2 / 3. .- 1.
    c20[findall(3.0 .<= F[:,11:20] .< 20)] = @view(F[:,11:20][findall(3.0 .<= F[:,11:20] .< 20)])/3.
    c20[findall(20. .<= F[:,11:20]      )] = @view(F[:,11:20][findall(20. .<= F[:,11:20]      )])*0. .+ (6 .+ 2/3.)
    cend = zeros(T, (n_t, 4))
    cend[findall(1.5 .<= F[:,21:end] .< 3 )] = @view(F[:,21:end][findall(1.5 .<= F[:,21:end] .< 3 )])/3. .- 0.5
    cend[findall(3.0 .<= F[:,21:end] .< 20)] = @view(F[:,21:end][findall(3.0 .<= F[:,21:end] .< 20)])/6.
    cend[findall(20. .<= F[:,21:end]      )] = @view(F[:,21:end][findall(20. .<= F[:,21:end]      )])*0. .+ (3. .+ 1/3.)
    C = hcat(c0, c10, c20, cend)

    return C
end

function f_pnlt(settings::PyObject, data::PyObject, spl)

    # Number of time steps
    n_t = size(spl)[1]

    # Compute noy
    T = eltype(spl)
    func_noy = LinearInterpolation((data.noy_spl, data.noy_freq), data.noy)

    N = func_noy.(spl, reshape(data.f, (1,settings.N_f)).*ones(T, (n_t,1)))
    # N = f_noy(settings, n_t, spl)

    # Compute pnl
    pnl = f_pnl(N)

    # Compute the correction factor C
    C = f_tone_corrections(settings, spl)

    # Step 10: Compute the largest of the tone correction
    if settings.TCF800
        c_max = maximum(C[:,14:end], dims=2)
    else
        c_max = maximum(C, dims=2)
    end

    # Compute PNLT
    pnlt = pnl + c_max

    return reshape(pnlt, (n_t,)), C
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
using Test
using Dates
using ScottishTaxBenefitModel
using .ModelHousehold: Household, Person, People_Dict, is_single,
    default_bu_allocation, get_benefit_units, get_head, get_spouse, search,
    pers_is_disabled, pers_is_carer
using .ExampleHouseholdGetter
using .Definitions
using .LegacyMeansTestedBenefits:
    tariff_income,
    LMTResults, is_working_hours, make_lmt_benefit_applicability,
    working_disabled, calc_allowances,
    apply_2_child_policy, calc_incomes
using .TimeSeriesUtils
using .Intermediate: MTIntermediate, make_intermediate

using .STBParameters: LegacyMeansTestedBenefitSystem, IncomeRules, HoursLimits,
    reached_state_pension_age, state_pension_age, AgeLimits
using .Results: LMTResults, LMTCanApplyFor
using .ExampleHelpers

## FIXME don't need both
lmt = LegacyMeansTestedBenefitSystem{Float64}()
sys = get_system( scotland=true )

@testset "Test Pension Ages" begin
    #
    # FIXME we need to jam a 'current' date on here or some of these will fail
    # next financial year
    #
    @test state_pension_age(sys.age_limits,Male) == 66

    @test state_pension_age(sys.age_limits,Female,2015) == 62
    @test state_pension_age(sys.age_limits,Female,2045) == 67
    @test state_pension_age(sys.age_limits,Female,2046) == 68
    @test state_pension_age(sys.age_limits,Male,2046) == 68

    @test reached_state_pension_age(sys.age_limits,66,Male)
    @test reached_state_pension_age(sys.age_limits,66,Female)
    @test ! reached_state_pension_age(sys.age_limits,62,Female, Date( 2022, 01, 01))
    # A person who's 67 now will have reached pension age
    # by then.
    @test reached_state_pension_age(sys.age_limits,67, Male, Date( 2046, 01, 01))
    # since this is financial year
    @test reached_state_pension_age(sys.age_limits,68, Male, 2046 )
    # The tests below have 'now' being 2021; tge `age_now` thing effectively jams in a Fixed
    # current date, by adding 1 to some of these ages next year and so on. Make Wholesale_trade_except_of_motor_vehicles_and_motorcycles
    # in June 2021, a 69 yo will have reached state pension age
    # but would have been 64 and so under state pension age (65) when the new state pension
    # was introduced in April 2016. Note all we have is years
    # for ages. So:
    # hadn't reached in April 2016
    # old style pension switch
    @test reached_state_pension_age(
        sys.age_limits,
        age_now(71),
        Male,
        sys.age_limits.savings_credit_to_new_state_pension )
    @test ! reached_state_pension_age(
        sys.age_limits,
        age_now(69),
        Male,
        sys.age_limits.savings_credit_to_new_state_pension )
    # reached pension age *now*
    @test reached_state_pension_age(
        sys.age_limits,
        age_now(69),
        Male )
    # 68 yo woman would have been 63 in 2016, so at state pension age (63 for women)
    @test reached_state_pension_age(
        sys.age_limits,
        age_now(69),
        Female,
        sys.age_limits.savings_credit_to_new_state_pension )
    # ..but a 67 yo would have been 62 so too young in '16 ...
    @test ! reached_state_pension_age(
        sys.age_limits,
        age_now(67),
        Female,
        sys.age_limits.savings_credit_to_new_state_pension )
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module ErrorFreeArith

export eftRecip, eftSquare, eftCube, eftSqrt, eftRecipSqrt, # see comment for eftSqrt
       eftSum2, eftSum2inOrder, eftDiff2, eftDiff2inOrder, eftProd2, eftDiv2,
       eftSum3, eftSum3inOrder, eftProd3, eftFMA, eftFMS,
       eftHorner

using Polynomials

if isdefined(:inlineErrorFreeArith) && (inlineErrorFreeArith == true)
   include("inlined.jl")
else
   include("funcall.jl")
end

end # module


text:
# Chance Sampled Monte Carlo Counterfactual Regret Minimization
# Johanson et al 2012 Algorithm 1
# http://modelai.gettysburg.edu/2013/cfr/cfr.pdf page 12 is also useful

# set up distributed processing

# add packages
using Pkg
Pkg.activate(pwd())
using Distributed
@everywhere using Random, StatsBase

# include other files
include("helper.jl")

# Define variables
N = 1:2 # number of players
I = 1:10 # information sets
A = [[1,2,3], [1,2], [1,2,3,4], [1,3,5], [1,2,6], [3,4,5], [4,6], [6], [4,5,6], [6]] # actions
actor = [1,2,1,2,1,2,1,1,1,1] # decision maker at each info set
ni = size(I,1) # number of info sets
na = sum.(A) # number of actions at each info set
regret = [zeros(na[i]) for i in 1:ni] # player 1 regret
strategy = [zeros(na[i]) for i in 1:ni] # cumulative strategies
σ1 = [fill(1 / na[i], na[i]) for i in 1:ni] # player 1 strategy
σ2 = [fill(1 / na[i], na[i]) for i in 1:ni] # player 2 strategy
σ = (σ1, σ2)

# sample tree
leaf1 = Node(1, [])
leaf2 = Node(3, [])
leaf3 = Node(4, [])
n = Node(2, [leaf2, leaf3])
head = Node(0, [leaf1, n])
nh = nnodes(head)
H, Z = gametree()

# constants
H, Z, C, P # all history nodes, terminal nodes, chance nodes, public nodes
S = [S[i] for i in N] # chance nodes only observable to player i
O = [O[i] for i in N] # chance nodes only observable to opponent of player i
# note: C = S1 ∪ O1 ∪ P and also C = S2 ∪ O2 ∪ P
Πi = [π1, π2] # contribution of player i to prob(h)
Πo = [πo1, πo2] # contribution of all other players than i to prob(h)
I1 = Vector{Vector{Node}}(undef, ni)
I2 = Vector{Vector{Node}}(undef, ni)
Iset = [I1, I2] # info sets
AI # dict from info set to action
Uz = [Uz1, Uz2] # utility of each leaf node for each player

function gametree
    # build h nodes into tree
    # record Z while doing this
    H, Z
end

function chance(h::Node, a::Int)
    1 / length(children(h)) # need to make this actually match public chance
end

function leafutility(h::Node, i::Int)
    Uz[i][idnode(h)]
end

function sampleoutcome(h::Node, πi::Float64, πo::Float64)
    if h ∈ C
        nextnodes = children(h)
        chance_weights = [chance(h, a) for a in actions]
        ha = sample(nextnodes, chance_weights)
        return updatetree!(ha, i, πi, πo)
    end
end

function infoset(h::Node)
    idinfo(h), Iset[idactor(h)][idinfo(h)]
end

function matchregret(I::Vector{Node})
    ind = findfirst(isequal(I), Iset)
    na = length(AI[ind]})
    σ = Vector{Float64}(undef, na)
    for a in AI[ind]
        denom = sum(b * max(r[ind][b], 0) for b in AI[ind])
        if denom > 0
            σ[a] = max(r[ind][a], 0) / denom
        else
            σ[a] = 1 / length(AI[ind])
        end
    end
    σ
end

function updateutility!(i::Int, ind::Int, h::Node, σ::Vector{Float64},
        u::Float64, πi::Float64, πo::Float64)
    na = length(AI[ind])
    m = Vector{Float64}(undef, na)
    for a in AI[ind]
        ind_ha = findfirst(x -> x == a, nextactions(h))
        ha = children(h)[ind_ha]
        if idactor(h) == i
            πip = σ[a] * πi
            up = updatetree(ha, i, πip, πo)
            m[a] = up
            u = u + σ[a] * up
        else
            πop = σ[a] * πo
            up = updatetree(ha, i, πi, πop)
            u = u + up
        end
    end
    m
end

function updateregret!(i::Int, ind::Int, σ::Vector{Float64},
        u::Float64, m::Vector{Float64}, πi::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}})
    if idactor(h) == i
        for a in AI[ind]
            r[ind][a] = r[ind][a] + m[a] - u
            s[ind][a] = s[ind][a] + πi * σ[a]
        end
    end
end

function updatetree(h::Node, i::Int, πi::Float64, πo::Float64)
    # πo is short for \vec{π}_{-i}
    h ∈ Z && return leafutility(h, i)
    sampleoutcome(h, πi, πo)
    indI, I = infoset(h)
    u = 0
    σ = matchregret(I)
    m = updateutility!(i, indI, h, σ, u, πi, πo)
    updateregret!(i, indI, σ, m, πi, r, s)
    u
end

function cfr(T::Int, ni::Int)
    r = [zeros(na[i]) for i in 1:ni] # regret
    s = [zeros(na[i]) for i in 1:ni] # cumulative strategies
    for t in 1:T
        for i in N
            u = updatetree!([], i, 1, 1, r, s) # utility
        end
    end
    u, r, s
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
"""
Fast multidimensional Chebyshev interpolation on a hypercube (Cartesian-product)
domain, using a separable (tensor-product) grid of Chebyshev interpolation points.

For domain upper and lower bounds `lb` and `ub`, and a given `order`
tuple, you would create an interpolator for a function `f` via:
```
using FastChebInterp
x = chebpoints(order, lb, ub) # an array of StaticVector
c = chebinterp(f.(x), lb, ub)
```
and then evaluate the interpolant for a point `y` (a vector)
via `c(y)`.

We also provide a function `chebgradient(c,y)` that returns a tuple `(c(y), ∇c(y))` of
the interpolant and its gradient at a point `y`.

The FastChebInterp package also supports complex and vector-valued functions `f`.  In
this case, `c(y)` returns a vector of interpolants, and one can use `chebjacobian(c,y)`
to compute the tuple `(c(y), J(y))` of the interpolant and its Jacobian matrix at `y`.
"""
module FastChebInterp

export chebpoints, chebinterp, chebinterp_v1, chebjacobian, chebgradient, chebregression

using StaticArrays
import FFTW

"""
    ChebPoly

A multidimensional Chebyshev-polynomial interpolation object.
Given a `c::ChebPoly`, you can evaluate it at a point `x`
with `c(x)`, where `x` is a vector (or a scalar if `c` is 1d).
"""
struct ChebPoly{N,T,Td<:Real}
    coefs::Array{T,N} # chebyshev coefficients
    lb::SVector{N,Td} # lower/upper bounds
    ub::SVector{N,Td} #    of the domain
end

function Base.show(io::IO, c::ChebPoly)
    print(io, "Chebyshev order ", map(i->i-1,size(c.coefs)), " polynomial on ",
          '[', c.lb[1], ',', c.ub[1], ']')
    for i = 2:length(c.lb)
        print(io, " × [", c.lb[i], ',', c.ub[i], ']')
    end
end
Base.ndims(c::ChebPoly) = ndims(c.coefs)

include("interp.jl")
include("regression.jl")
include("eval.jl")

end # module


text:
"""
    struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}

A ProductSpace is a tensor product space of `N` vector spaces of type
`S<:ElementarySpace`. Only tensor products between `ElementarySpace` objects of
the same type are allowed.
"""
struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}
    spaces::NTuple{N, S}
end
ProductSpace(spaces::NTuple{N,S}) where {S<:ElementarySpace, N} = ProductSpace{S,N}(spaces)
ProductSpace(spaces::Vararg{S,N}) where {S<:ElementarySpace, N} = ProductSpace{S,N}(spaces)

# Corresponding methods
#-----------------------
dim(P::ProductSpace{<:ElementarySpace,0}) = 1
dim(P::ProductSpace) = prod(dim, P.spaces)
dim(P::ProductSpace, n::Int) = dim(P.spaces[n], args...)

Base.indices(P::ProductSpace) = CartesianRange(map(indices, P.spaces))
Base.indices(P::ProductSpace, n::Int) = indices(P.spaces[n])

dual(P::ProductSpace{<:ElementarySpace,0}) = P
dual(P::ProductSpace) = ProductSpace(map(dual, reverse(P.spaces)))

# Base.conj(P::ProductSpace) = ProductSpace(map(conj, P.spaces))

function Base.show(io::IO, P::ProductSpace)
    spaces = P.spaces
    if length(spaces) == 1
        print(io,"ProductSpace")
    end
    print(io,"(")
    for i in 1:length(spaces)
        i==1 || print(io," ⊗ ")
        show(io, spaces[i])
    end
    print(io,")")
end

# more specific methods
sectors(P::ProductSpace) = _sectors(P, sectortype(P))
_sectors(P::ProductSpace{<:ElementarySpace, N}, ::Type{Trivial}) where {N} = (ntuple(n->Trivial(),Val{N}()),) # speed up sectors for ungraded spaces
_sectors(P::ProductSpace{<:ElementarySpace, N}, ::Type{<:Sector}) where {N} = product(map(sectors, P.spaces)...)
dim(P::ProductSpace{<:ElementarySpace, N}, sector::NTuple{N, Sector}) where {N} = prod(map(dim, P.spaces, sector))

Base.indices(P::ProductSpace{<:AbstractRepresentationSpace{G}, N}, sectors::NTuple{N, G}) where {G<:Sector, N} =
        CartesianRange(map(indices, P.spaces, sectors))

# Default construction from product of spaces
#---------------------------------------------
⊗(V1::S, V2::S) where {S<:ElementarySpace}= ProductSpace((V1, V2))
⊗(P1::ProductSpace{S}, V2::S) where {S<:ElementarySpace} = ProductSpace(tuple(P1.spaces..., V2))
⊗(V1::S, P2::ProductSpace{S}) where {S<:ElementarySpace} = ProductSpace(tuple(V1, P2.spaces...))
⊗(P1::ProductSpace{S}, P2::ProductSpace{S}) where {S<:ElementarySpace} = ProductSpace(tuple(P1.spaces..., P2.spaces...))
⊗(P::ProductSpace{S}, ::ProductSpace{S,0}) where {S<:ElementarySpace} = P1

# Functionality for extracting and iterating over spaces
#--------------------------------------------------------
Base.length(P::ProductSpace) = length(P.spaces)
Base.getindex(P::ProductSpace, n::Integer) = P.spaces[n]
Base.getindex(P::ProductSpace, r::AbstractVector{<:Integer}) = ProductSpace(P.spaces[r])
#
Base.start(P::ProductSpace) = start(P.spaces)
Base.next(P::ProductSpace, state) = next(P.spaces, state)
Base.done(P::ProductSpace, state) = done(P.spaces, state)

Base.iteratoreltype(P::ProductSpace) = Base.iteratoreltype(P.spaces)
Base.iteratorsize(P::ProductSpace) = Base.iteratorsize(P.spaces)


text:
struct DGModel{BL,G,NFND,NFD,GNF,AS,DS,D,MD}
  balancelaw::BL
  grid::G
  numfluxnondiff::NFND
  numfluxdiff::NFD
  gradnumflux::GNF
  auxstate::AS
  diffstate::DS
  direction::D
  modeldata::MD
end
function DGModel(balancelaw, grid, numfluxnondiff, numfluxdiff, gradnumflux;
                 auxstate=create_auxstate(balancelaw, grid),
                 diffstate=create_diffstate(balancelaw, grid),
                 direction=EveryDirection(), modeldata=nothing)
  DGModel(balancelaw, grid, numfluxnondiff, numfluxdiff, gradnumflux, auxstate,
          diffstate, direction, modeldata)
end

function (dg::DGModel)(dQdt, Q, ::Nothing, t; increment=false)
  bl = dg.balancelaw
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Nfp = Nq * Nqk
  nrealelem = length(topology.realelems)

  Qvisc = dg.diffstate
  auxstate = dg.auxstate

  FT = eltype(Q)
  nviscstate = num_diffusive(bl, FT)

  Np = dofs_per_element(grid)

  communicate = !(isstacked(topology) &&
                  typeof(dg.direction) <: VerticalDirection)

  aux_comm = update_aux!(dg, bl, Q, t)
  @assert typeof(aux_comm) == Bool

  ########################
  # Gradient Computation #
  ########################
  if communicate
    MPIStateArrays.start_ghost_exchange!(Q)
    if aux_comm
      MPIStateArrays.start_ghost_exchange!(auxstate)
    end
  end

  if nviscstate > 0

    @launch(device, threads=(Nq, Nq, Nqk), blocks=nrealelem,
            volumeviscterms!(bl, Val(dim), Val(N), dg.direction, Q.data,
                             Qvisc.data, auxstate.data, grid.vgeo, t, grid.D,
                             topology.realelems))

    if communicate
      MPIStateArrays.finish_ghost_recv!(Q)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    end

    @launch(device, threads=Nfp, blocks=nrealelem,
            faceviscterms!(bl, Val(dim), Val(N), dg.direction,
                           dg.gradnumflux, Q.data, Qvisc.data, auxstate.data,
                           grid.vgeo, grid.sgeo, t, grid.vmapM, grid.vmapP, grid.elemtobndy,
                           topology.realelems))

    if communicate
      MPIStateArrays.start_ghost_exchange!(Qvisc)
    end

    aux_comm = update_aux_diffusive!(dg, bl, Q, t)
    @assert typeof(aux_comm) == Bool

    if aux_comm
      MPIStateArrays.start_ghost_exchange!(auxstate)
    end
  end

  ###################
  # RHS Computation #
  ###################
  @launch(device, threads=(Nq, Nq, Nqk), blocks=nrealelem,
          volumerhs!(bl, Val(dim), Val(N), dg.direction, dQdt.data,
                     Q.data, Qvisc.data, auxstate.data, grid.vgeo, t,
                     grid.ω, grid.D, topology.realelems, increment))

  if communicate
    if nviscstate > 0
      MPIStateArrays.finish_ghost_recv!(Qvisc)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    else
      MPIStateArrays.finish_ghost_recv!(Q)
      if aux_comm
        MPIStateArrays.finish_ghost_recv!(auxstate)
      end
    end
  end

  @launch(device, threads=Nfp, blocks=nrealelem,
          facerhs!(bl, Val(dim), Val(N), dg.direction,
                   dg.numfluxnondiff,
                   dg.numfluxdiff,
                   dQdt.data, Q.data, Qvisc.data,
                   auxstate.data, grid.vgeo, grid.sgeo, t, grid.vmapM, grid.vmapP, grid.elemtobndy,
                   topology.realelems))

  # Just to be safe, we wait on the sends we started.
  if communicate
    MPIStateArrays.finish_ghost_send!(Qvisc)
    MPIStateArrays.finish_ghost_send!(Q)
  end
end

function init_ode_state(dg::DGModel, args...;
                        forcecpu=false,
                        commtag=888)
  device = arraytype(dg.grid) <: Array ? CPU() : CUDA()

  bl = dg.balancelaw
  grid = dg.grid

  state = create_state(bl, grid, commtag)

  topology = grid.topology
  Np = dofs_per_element(grid)

  auxstate = dg.auxstate
  dim = dimensionality(grid)
  N = polynomialorder(grid)
  nrealelem = length(topology.realelems)

  if !forcecpu
    @launch(device, threads=(Np,), blocks=nrealelem,
            initstate!(bl, Val(dim), Val(N), state.data, auxstate.data, grid.vgeo,
                     topology.realelems, args...))
  else
    h_state = similar(state, Array)
    h_auxstate = similar(auxstate, Array)
    h_auxstate .= auxstate
    @launch(CPU(), threads=(Np,), blocks=nrealelem,
      initstate!(bl, Val(dim), Val(N), h_state.data, h_auxstate.data, Array(grid.vgeo),
          topology.realelems, args...))
    state .= h_state
  end

  MPIStateArrays.start_ghost_exchange!(state)
  MPIStateArrays.finish_ghost_exchange!(state)

  return state
end

# fallback
function update_aux!(dg::DGModel, bl::BalanceLaw, Q::MPIStateArray, t::Real)
  return false
end

function update_aux_diffusive!(dg::DGModel, bl::BalanceLaw, Q::MPIStateArray, t::Real)
  return false
end

function indefinite_stack_integral!(dg::DGModel, m::BalanceLaw,
                                    Q::MPIStateArray,
                                    auxstate::MPIStateArray,
                                    t::Real)

  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  FT = eltype(Q)

  # do integrals
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_indefinite_stack_integral!(m, Val(dim), Val(N),
                                         Val(nvertelem),
                                         Q.data, auxstate.data,
                                         grid.vgeo, grid.Imat, 1:nhorzelem))
end

function reverse_indefinite_stack_integral!(dg::DGModel,
                                            m::BalanceLaw,
                                            Q::MPIStateArray,
                                            auxstate::MPIStateArray, t::Real)

  device = typeof(auxstate.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  FT = eltype(auxstate)

  # do integrals
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_reverse_indefinite_stack_integral!(m, Val(dim), Val(N),
                                                 Val(nvertelem),
                                                 Q.data, auxstate.data,
                                                 1:nhorzelem))
end

function nodal_update_aux!(f!, dg::DGModel, m::BalanceLaw, Q::MPIStateArray,
                           t::Real; diffusive=false)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  nrealelem = length(topology.realelems)

  Np = dofs_per_element(grid)

  ### update aux variables
  if diffusive
    @launch(device, threads=(Np,), blocks=nrealelem,
            knl_nodal_update_aux!(m, Val(dim), Val(N), f!,
                            Q.data, dg.auxstate.data, dg.diffstate.data, t,
                            topology.realelems))
  else
    @launch(device, threads=(Np,), blocks=nrealelem,
            knl_nodal_update_aux!(m, Val(dim), Val(N), f!,
                            Q.data, dg.auxstate.data, t,
                            topology.realelems))
  end
end

"""
    courant(local_courant::Function, dg::DGModel, m::BalanceLaw,
            Q::MPIStateArray, direction=EveryDirection())
Returns the maximum of the evaluation of the function `local_courant`
pointwise throughout the domain.  The function `local_courant` is given an
approximation of the local node distance `Δx`.  The `direction` controls which
reference directions are considered when computing the minimum node distance
`Δx`.
An example `local_courant` function is
    function local_courant(m::AtmosModel, state::Vars, aux::Vars,
                           diffusive::Vars, Δx)
      return Δt * cmax / Δx
    end
where `Δt` is the time step size and `cmax` is the maximum flow speed in the
model.
"""
function courant(local_courant::Function, dg::DGModel, m::BalanceLaw,
                 Q::MPIStateArray, Δt, direction=EveryDirection())
    grid = dg.grid
    topology = grid.topology
    nrealelem = length(topology.realelems)

    if nrealelem > 0
        N = polynomialorder(grid)
        dim = dimensionality(grid)
        Nq = N + 1
        Nqk = dim == 2 ? 1 : Nq
        device = grid.vgeo isa Array ? CPU() : CUDA()
        pointwise_courant = similar(grid.vgeo, Nq^dim, nrealelem)
        @launch(device, threads=(Nq, Nq, Nqk), blocks=nrealelem,
        Grids.knl_min_neighbor_distance!(Val(N), Val(dim), direction,
                                         pointwise_courant, grid.vgeo, topology.realelems))
        @launch(device, threads=(Nq*Nq*Nqk,), blocks=nrealelem,
                knl_local_courant!(m, Val(dim), Val(N), pointwise_courant,
                local_courant, Q.data, dg.auxstate.data,
                dg.diffstate.data, topology.realelems, direction, Δt))
        rank_courant_max = maximum(pointwise_courant)
    else
        rank_courant_max = typemin(eltype(Q))
    end

    MPI.Allreduce(rank_courant_max, max, topology.mpicomm)
end

function copy_stack_field_down!(dg::DGModel, m::BalanceLaw,
                                auxstate::MPIStateArray, fldin, fldout)

  device = typeof(auxstate.data) <: Array ? CPU() : CUDA()

  grid = dg.grid
  topology = grid.topology

  dim = dimensionality(grid)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  # do integrals
  nelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nelem, nvertelem)

  @launch(device, threads=(Nq, Nqk, 1), blocks=nhorzelem,
          knl_copy_stack_field_down!(Val(dim), Val(N), Val(nvertelem),
                                     auxstate.data, 1:nhorzelem, Val(fldin),
                                     Val(fldout)))
end

function MPIStateArrays.MPIStateArray(dg::DGModel, commtag=888)
  bl = dg.balancelaw
  grid = dg.grid

  state = create_state(bl, grid, commtag)

  return state
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using Distributions, ForwardDiff, Plots, LinearAlgebra

replics = 1000

x = randn(100)

y = 1.5*x + randn(100)

prior = Normal(0,1)

par = zeros(replics)

likelihood(y,x,par) = sum(logpdf.(Normal(0,1),y - x .*par))

loglike(y,x,par) = likelihood(y,x,par)+logpdf.(prior,par)

unif = Uniform(0,1)

M = I(9)

M_inv = inv(M)

moment_dist = Normal(0,1.5)

par[1] = rand(prior)

for i in 2:replics
    moment = rand(moment_dist)
    old_moment = moment
    par_aux = par[i-1]

    L = ceil(rand(Uniform(0,500)))
    leapfrog_step = 1/L

    for j in 1:L
        new_moment = moment + 1/2*leapfrog_step*ForwardDiff.derivative(pp->loglike(y,x,pp),par_aux)
         par_aux = par_aux .+ leapfrog_step*1*new_moment
         moment = new_moment + 1/2*leapfrog_step*ForwardDiff.derivative(pp->loglike(y,x,pp),par_aux)
    end

    num = loglike(y,x,par_aux) + logpdf(moment_dist,moment)#+ logpdf(kernel_novo,par[i-1])

    dem = loglike(y,x,par[i-1]) + logpdf(moment_dist,old_moment)#+ logpdf(kernel,novo_par)
    alpha = exp(min(0,num - dem))
    p = rand(unif)
    if alpha < p
        par[i] = par[i-1]
    else
        par[i] = par_aux
    end
    println(i)
end

histogram(par[100:1000])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DeploymentPropertiesExtended <: SwaggerModel
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    correlationId::Any # spec type: Union{ Nothing, String } # spec name: correlationId
    timestamp::Any # spec type: Union{ Nothing, String } # spec name: timestamp
    outputs::Any # spec type: Union{ Nothing, Any } # spec name: outputs
    providers::Any # spec type: Union{ Nothing, Vector{Provider} } # spec name: providers
    dependencies::Any # spec type: Union{ Nothing, Vector{Dependency} } # spec name: dependencies
    template::Any # spec type: Union{ Nothing, Any } # spec name: template
    templateLink::Any # spec type: Union{ Nothing, TemplateLink } # spec name: templateLink
    parameters::Any # spec type: Union{ Nothing, Any } # spec name: parameters
    parametersLink::Any # spec type: Union{ Nothing, ParametersLink } # spec name: parametersLink
    mode::Any # spec type: Union{ Nothing, String } # spec name: mode
    debugSetting::Any # spec type: Union{ Nothing, DebugSetting } # spec name: debugSetting

    function DeploymentPropertiesExtended(;provisioningState=nothing, correlationId=nothing, timestamp=nothing, outputs=nothing, providers=nothing, dependencies=nothing, template=nothing, templateLink=nothing, parameters=nothing, parametersLink=nothing, mode=nothing, debugSetting=nothing)
        o = new()
        validate_property(DeploymentPropertiesExtended, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(DeploymentPropertiesExtended, Symbol("correlationId"), correlationId)
        setfield!(o, Symbol("correlationId"), correlationId)
        validate_property(DeploymentPropertiesExtended, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(DeploymentPropertiesExtended, Symbol("outputs"), outputs)
        setfield!(o, Symbol("outputs"), outputs)
        validate_property(DeploymentPropertiesExtended, Symbol("providers"), providers)
        setfield!(o, Symbol("providers"), providers)
        validate_property(DeploymentPropertiesExtended, Symbol("dependencies"), dependencies)
        setfield!(o, Symbol("dependencies"), dependencies)
        validate_property(DeploymentPropertiesExtended, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(DeploymentPropertiesExtended, Symbol("templateLink"), templateLink)
        setfield!(o, Symbol("templateLink"), templateLink)
        validate_property(DeploymentPropertiesExtended, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(DeploymentPropertiesExtended, Symbol("parametersLink"), parametersLink)
        setfield!(o, Symbol("parametersLink"), parametersLink)
        validate_property(DeploymentPropertiesExtended, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(DeploymentPropertiesExtended, Symbol("debugSetting"), debugSetting)
        setfield!(o, Symbol("debugSetting"), debugSetting)
        o
    end
end # type DeploymentPropertiesExtended

const _property_map_DeploymentPropertiesExtended = Dict{Symbol,Symbol}(Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("correlationId")=>Symbol("correlationId"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("outputs")=>Symbol("outputs"), Symbol("providers")=>Symbol("providers"), Symbol("dependencies")=>Symbol("dependencies"), Symbol("template")=>Symbol("template"), Symbol("templateLink")=>Symbol("templateLink"), Symbol("parameters")=>Symbol("parameters"), Symbol("parametersLink")=>Symbol("parametersLink"), Symbol("mode")=>Symbol("mode"), Symbol("debugSetting")=>Symbol("debugSetting"))
const _property_types_DeploymentPropertiesExtended = Dict{Symbol,String}(Symbol("provisioningState")=>"String", Symbol("correlationId")=>"String", Symbol("timestamp")=>"String", Symbol("outputs")=>"Any", Symbol("providers")=>"Vector{Provider}", Symbol("dependencies")=>"Vector{Dependency}", Symbol("template")=>"Any", Symbol("templateLink")=>"TemplateLink", Symbol("parameters")=>"Any", Symbol("parametersLink")=>"ParametersLink", Symbol("mode")=>"String", Symbol("debugSetting")=>"DebugSetting")
Base.propertynames(::Type{ DeploymentPropertiesExtended }) = collect(keys(_property_map_DeploymentPropertiesExtended))
Swagger.property_type(::Type{ DeploymentPropertiesExtended }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DeploymentPropertiesExtended[name]))}
Swagger.field_name(::Type{ DeploymentPropertiesExtended }, property_name::Symbol) =  _property_map_DeploymentPropertiesExtended[property_name]

const _allowed_DeploymentPropertiesExtended_mode = ["Incremental", "Complete"]

function check_required(o::DeploymentPropertiesExtended)
    true
end

function validate_property(::Type{ DeploymentPropertiesExtended }, name::Symbol, val)
    if name === Symbol("mode")
        Swagger.validate_param(name, "DeploymentPropertiesExtended", :enum, val, _allowed_DeploymentPropertiesExtended_mode)
    end
end


text:
function init_boxes(f, x0::AbstractVector{T}, splits, u::AbstractVector{T}, v::AbstractVector{T}) where T
    xstar = copy(x0)
    fstar = f(xstar)
    root = box = Box{T}()
    n = length(x0)
    length(splits) == n || throw(DimensionMismatch("need one vector of splits per dimension, got $(length(splits)) for $n dimensions"))
    xtmp = copy(xstar)
    for i = 1:n
        xtmp = ipcopy!(xtmp, xstar)
        xsplit = splits[i]
        L = length(xsplit)
        L >= 3 || error("there have to be at least 3 evaluations per coordinate")
        issorted(xsplit) || error("splits must be in increasing order")
        # Evaluate f along the splits, keeping track of the best
        fsplit = Vector{typeof(fstar)}(uninitialized, L)
        fmin, idxmin = oftype(fstar, Inf), 0
        for l = 1:L
            xtmp = replacecoordinate!(xtmp, i, xsplit[l])
            ftmp = f(xtmp)
            if ftmp < fmin
                fmin = ftmp
                idxmin = l
            end
            fsplit[l] = ftmp
        end
        box.splitdim = i
        box.xvalues = copy(xsplit)
        box.fvalues = fsplit
        idxmin == 0 && error("function was not finite at any evaluation point")
        # Create the child boxes
        overhang_left = xsplit[1] > u[i]
        idx = Int(!overhang_left)
        level = box.level + 1
        while idx <= L
            if idx == 0
                # beyond the edge of xsplit
                Box{T}(box, level, u[i])  # we can drop this because it gets stored in the parent
            elseif idx == L
                if xsplit[end] < v[i]
                    Box{T}(box, level, v[i]) # also beyond edge
                end
            else
                if fsplit[idx] < fsplit[idx+1]
                    # Bigger box (with lower level) is closer to the smaller f value
                    edge = xsplit[idx] + (xsplit[idx+1]-xsplit[idx]) * ((sqrt(5)-1)/2)
                    Box{T}(box, level, edge)
                    Box{T}(box, level+1, edge)
                else
                    edge = xsplit[idx+1] - (xsplit[idx+1]-xsplit[idx]) * ((sqrt(5)-1)/2)
                    Box{T}(box, level+1, edge)
                    Box{T}(box, level, edge)
                end
            end
            idx += 1
        end
        # Update xstar
        if fmin < fstar
            fstar = fmin
            xstar = replacecoordinate!(xstar, i, xsplit[idxmin])
        end
        # Pick the best box for splitting for the next coordinate
        im, i0, ip = idxmin == 1 ? (1, 2, 3) :
                     idxmin == L ? (L-2, L-1, L) :
                     (idxmin-1, idxmin, idxmin+1)
        xqmin, _ = qfit_min((xsplit[im], fsplit[im]), (xsplit[i0], fsplit[i0]), (xsplit[ip], fsplit[ip]))
        cindex = xqmin < xsplit[im] ? (assert(idxmin == 1); 1) : # at left edge
                 xqmin > xsplit[ip] ? (assert(idxmin == L); length(box.children)) :  # at right edge
                 xqmin < xsplit[i0] ? 2*idxmin-2+overhang_left : 2*idxmin-1+overhang_left
        box = box.children[cindex]
    end
    xstar, fstar, box
end

function init_sweep_queues(lastbox::Box{T}, smax) where T
    root = get_root(lastbox)
    pqs = [PriorityQueue{Box{T},T}() for s = 1:smax]
    for box in visit_leaves(root)
        pqs[box.level][box] = minimum(box)
    end
    pqs
end

function run_sweep(pqs)
    for pq in pqs
        isempty(pq) && continue
        box = dequeue!(pq)
        if !split!(box)
            if box.level < length(pqs)
                box.level += 1
                pqs[box.level][box] = minimum(box)
            end
        end
    end
end


text: Base.@kwdef mutable struct model1 <: ReactiveModel


text:
# Tests the routine phisalpha jacobian.  This routine
# computes the force gradient correction after Dehnen & Hernandez (2017).

import NbodyGradient: phisalpha!, init_nbody

# Next, try computing three-body Keplerian Jacobian:

@testset "phisalpha" begin
n = 3
H = [3,1,1]
t0 = 7257.93115525
h  = 0.05
tmax = 600.0
dlnq = big(1e-15)

elements = readdlm("elements.txt", ',')

m = zeros(n)
x0 = zeros(3, n)
v0 = zeros(3, n)
alpha = 2.0

# Define which pairs will have impulse rather than -drift+Kepler:
pair = zeros(Bool, n, n)
for k = 1:n
    m[k] = elements[k,1]
end
m0 = copy(m)

init = ElementsIC(t0, H, elements)
ic_big = ElementsIC(big(t0), H, big.(elements))
x0, v0, _ = init_nbody(init)

# Tilt the orbits a bit:
x0[2,1] = 5e-1 * sqrt(x0[1,1]^2 + x0[3,1]^2)
x0[2,2] = -5e-1 * sqrt(x0[1,2]^2 + x0[3,2]^2)
x0[2,3] = -5e-1 * sqrt(x0[1,2]^2 + x0[3,2]^2)
v0[2,1] = 5e-1 * sqrt(v0[1,1]^2 + v0[3,1]^2)
v0[2,2] = -5e-1 * sqrt(v0[1,2]^2 + v0[3,2]^2)
v0[2,3] = -5e-1 * sqrt(v0[1,2]^2 + v0[3,2]^2)

# Save to state
s = State(init)
s.x .= x0; s.v .= v0
s.pair .= pair

# Take a step:
ahl21!(s,h)
x0 .= s.x; v0 .= s.v

# Now, copy these to compute Jacobian (so that I don't step
# x0 & v0 forward in time):
x = copy(s.x); v = copy(s.v); m = copy(s.m)
# Compute jacobian exactly:
d = Derivatives(Float64, s.n)
phisalpha!(s,d,h,alpha)
d.jac_phi .+= Matrix{Float64}(I, size(d.jac_phi))

# Now compute numerical derivatives, using BigFloat to avoid
# round-off errors:
jac_step_num = zeros(BigFloat, 7 * n, 7 * n)
# Save these so that I can compute derivatives numerically:
xsave = big.(x0)
vsave = big.(v0)
msave = big.(m0)
hbig = big(h)
abig = big(alpha)
sbig = deepcopy(State(ic_big))
sbig.x .= xsave; sbig.v .= vsave; sbig.m .= msave
# Carry out step using BigFloat for extra precision:
phisalpha!(sbig,hbig,abig)
# Copy back to saves
xsave .= sbig.x; vsave .= sbig.v; msave .= sbig.m
# Compute numerical derivatives wrt time:
dqdt_num = zeros(BigFloat, 7 * n)
# Initial positions, velocities & masses:
hbig = big(h)
dq = dlnq * hbig
hbig += dq
sbig = deepcopy(State(ic_big))
sbig.x .= big.(x0); sbig.v .= big.(v0); sbig.m .= big.(m0)
phisalpha!(sbig,hbig,abig)

# Now x & v are final positions & velocities after time step
for i in 1:n, k in 1:3
    dqdt_num[(i - 1) * 7 +  k] = (sbig.x[k,i] - xsave[k,i]) / dq
    dqdt_num[(i - 1) * 7 + 3 + k] = (sbig.v[k,i] - vsave[k,i]) / dq
end

hbig = big(h)
# Vary the initial parameters of planet j:
for j = 1:n
    # Vary the initial phase-space elements:
    for jj = 1:3
        # Initial positions, velocities & masses:
        sbig = deepcopy(State(ic_big))
        sbig.x .= big.(x0); sbig.v .= big.(v0); sbig.m .= big.(m0)
        dq = dlnq * sbig.x[jj,j]
        if sbig.x[jj,j] != 0.0
            sbig.x[jj,j] +=  dq
        else
            dq = dlnq
            sbig.x[jj,j] = dq
        end
        phisalpha!(sbig, hbig, abig)
        # Now x & v are final positions & velocities after time step
        for i = 1:n
            for k = 1:3
                jac_step_num[(i - 1) * 7 +  k,(j - 1) * 7 + jj] = (sbig.x[k,i] - xsave[k,i]) / dq
                jac_step_num[(i - 1) * 7 + 3 + k,(j - 1) * 7 + jj] = (sbig.v[k,i] - vsave[k,i]) / dq
            end
        end

        sbig = deepcopy(State(ic_big))
        sbig.x .= big.(x0); sbig.v .= big.(v0); sbig.m .= big.(m0)
        dq = dlnq * sbig.v[jj,j]
        if sbig.v[jj,j] != 0.0
            sbig.v[jj,j] +=  dq
        else
            dq = dlnq
            sbig.v[jj,j] = dq
        end
        phisalpha!(sbig, hbig, abig)

        for i = 1:n
            for k = 1:3
                jac_step_num[(i - 1) * 7 +  k,(j - 1) * 7 + 3 + jj] = (sbig.x[k,i] - xsave[k,i]) / dq
                jac_step_num[(i - 1) * 7 + 3 + k,(j - 1) * 7 + 3 + jj] = (sbig.v[k,i] - vsave[k,i]) / dq
            end
        end
    end
    # Now vary mass of planet:
    sbig = deepcopy(State(ic_big))
    sbig.x .= big.(x0); sbig.v .= big.(v0); sbig.m .= big.(m0)
    dq = sbig.m[j] * dlnq
    sbig.m[j] += dq
    phisalpha!(sbig, hbig, abig)
    for i = 1:n
        for k = 1:3
            jac_step_num[(i - 1) * 7 +  k,j * 7] = (sbig.x[k,i] - xsave[k,i]) / dq
            jac_step_num[(i - 1) * 7 + 3 + k,j * 7] = (sbig.v[k,i] - vsave[k,i]) / dq
        end
        # Mass unchanged -> identity
        jac_step_num[7 * i,7 * i] = big(1.0)
    end
end

# Now, compare the results:
dqdt_num = convert(Array{Float64,1}, dqdt_num)
@test isapprox(d.jac_phi, jac_step_num;norm=maxabs)
@test isapprox(d.dqdt_phi, dqdt_num;norm=maxabs)
end


text:
"""
    CapsuleCapsuleCollision

    collision between two capsules

    origin_parent: radius of contact for parent body
    origin_child: radius of contact for child body
    orientation_parent: radius of contact for parent body
    orientation_child: radius of contact for child body
    radius_parent: radius of contact for parent body
    radius_child: radius of contact for child body
    height_parent: radius of contact for parent body
    height_child: radius of contact for child body
"""
mutable struct CapsuleCapsuleCollision{T,O,I,OI} <: Collision{T,O,I,OI}
    origin_parent::SVector{I,T}
    origin_child::SVector{I,T}
    orientation_parent::Quaternion{T}
    orientation_child::Quaternion{T}
    radius_parent::T
    radius_child::T
    height_parent::T
    height_child::T
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, collision::CapsuleCapsuleCollision)
#     summary(io, collision)
#     println(io, "")
#     println(io, "origin_parent:      "*string(collision.origin_parent))
#     println(io, "origin_child:       "*string(collision.origin_child))
#     println(io, "orientation_parent: "*string(collision.orientation_parent))
#     println(io, "orientation_child:  "*string(collision.orientation_child))
#     println(io, "radius_parent:      "*string(collision.radius_parent))
#     println(io, "radius_child:       "*string(collision.radius_child))
#     println(io, "height_parent:      "*string(collision.height_parent))
#     println(io, "height_child:       "*string(collision.height_child))
# end

# distance
function distance(collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    return norm(pp - pc, 2)
end

function ∂distance∂x(gradient::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    X = ∂norm∂x(d) * (∂contact_point∂x(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂x(:child, jacobian, collision, xp, qp, xc, qc))

    return X
end

function ∂distance∂q(gradient::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # contact points
    pp = contact_point(:parent, collision, xp, qp, xc, qc)
    pc = contact_point(:child,  collision, xp, qp, xc, qc)

    # difference
    d = pp - pc

    X = ∂norm∂x(d) * (∂contact_point∂q(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂q(:child, jacobian, collision, xp, qp, xc, qc))

    return X
end

# contact point in world frame
function contact_point(relative::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # call mehrotra

    if relative == :parent
        return collision.ip.z[0 .+ (1:3)]
    elseif relative == :child
        return collision.ip.z[3 .+ (1:3)]
    end
end

function ∂contact_point∂x(relative::Symbol, jacobian::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
    # call mehrotra

    if relative == :parent
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    elseif relative == :child
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    end

    return X
end

function ∂contact_point∂q(relative::Symbol, jacobian::Symbol, collision::CapsuleCapsuleCollision, xp, qp, xc, qc)
# call mehrotra

    if relative == :parent
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    elseif relative == :child
        if jacobian == :parent
            X = collision.ip.δz[]
        elseif jacobian == :child
            X = collision.ip.δz[]
        end
    end

    return X
end


text:
using aoc
using Test

@testset "aoc.jl" begin
    @test aoc.solve(aoc.day1) == 858496
    @test aoc.solve(aoc.day2) == 396
    @test aoc.solve(aoc.day3) == 169
end


text:
using StatsBase
using LinearAlgebra, Random, SparseArrays, Test

@testset "StatsBase.Weights" begin
weight_funcs = (weights, aweights, fweights, pweights)

# Construction
@testset "$f" for f in weight_funcs
    @test isa(f([1, 2, 3]), AbstractWeights{Int})
    @test isa(f([1., 2., 3.]), AbstractWeights{Float64})
    @test isa(f([1 2 3; 4 5 6]), AbstractWeights{Int})

    @test isempty(f(Float64[]))
    @test size(f([1, 2, 3])) == (3,)

    w  = [1., 2., 3.]
    wv = f(w)
    @test eltype(wv) === Float64
    @test length(wv) === 3
    @test values(wv) === w
    @test sum(wv) === 6.0
    @test !isempty(wv)

    b  = trues(3)
    bv = f(b)
    @test eltype(bv) === Bool
    @test length(bv) === 3
    @test values(bv) === b
    @test sum(bv)    === 3
    @test !isempty(bv)

    ba = BitArray([true, false, true])
    sa = sparsevec([1., 0., 2.])

    @test sum(ba, wv) === 4.0
    @test sum(sa, wv) === 7.0
end

@testset "$f, setindex!" for f in weight_funcs
    w = [1., 2., 3.]
    wv = f(w)

    # Check getindex & sum
    @test wv[1] === 1.
    @test sum(wv) === 6.
    @test values(wv) == w

    # Test setindex! success
    @test (wv[1] = 4) === 4             # setindex! returns original val
    @test wv[1] === 4.                  # value correctly converted and set
    @test sum(wv) === 9.                # sum updated
    @test values(wv) == [4., 2., 3.]    # Test state of all values

    # Test mulivalue setindex!
    wv[1:2] = [3., 5.]
    @test wv[1] === 3.
    @test wv[2] === 5.
    @test sum(wv) === 11.
    @test values(wv) == [3., 5., 3.]   # Test state of all values

    # Test failed setindex! due to conversion error
    w = [1, 2, 3]
    wv = f(w)

    @test_throws InexactError wv[1] = 1.5   # Returns original value
    @test wv[1] === 1                       # value not updated
    @test sum(wv) === 6                     # sum not corrupted
    @test values(wv) == [1, 2, 3]           # Test state of all values
end

@testset "$f, isequal and ==" for f in weight_funcs
    x = f([1, 2, 3])

    y = f([1, 2, 3]) # same values, type and parameters
    @test isequal(x, y)
    @test x == y

    y = f([1.0, 2.0, 3.0]) # same values and type, different parameters
    @test isequal(x, y)
    @test x == y

    if f != fweights # same values and parameters, different types
        y = fweights([1, 2, 3])
        @test !isequal(x, y)
        @test x != y
    end

    x = f([1, 2, NaN]) # isequal and == treat NaN differently
    y = f([1, 2, NaN])
    @test isequal(x, y)
    @test x != y

    x = f([1.0, 2.0, 0.0]) # isequal and == treat ±0.0 differently
    y = f([1.0, 2.0, -0.0])
    @test !isequal(x, y)
    @test x == y
end

## wsum
x = [6., 8., 9.]
w = [2., 3., 4.]
p = [1. 2. ; 3. 4.]
q = [1., 2., 3., 4.]

@test wsum(Float64[], Float64[]) === 0.0
@test wsum(x, w) === 72.0
@test wsum(p, q) === 29.0

## wsum along dimension
@test wsum(x, w, 1) == [72.0]

x  = rand(6, 8)
w1 = rand(6)
w2 = rand(8)

@test size(wsum(x, w1, 1)) == (1, 8)
@test size(wsum(x, w2, 2)) == (6, 1)

@test wsum(x, w1, 1) ≈ sum(x .* w1, dims = 1)
@test wsum(x, w2, 2) ≈ sum(x .* w2', dims = 2)

x = rand(6, 5, 4)
w1 = rand(6)
w2 = rand(5)
w3 = rand(4)

@test size(wsum(x, w1, 1)) == (1, 5, 4)
@test size(wsum(x, w2, 2)) == (6, 1, 4)
@test size(wsum(x, w3, 3)) == (6, 5, 1)

@test wsum(x, w1, 1) ≈ sum(x .* w1, dims = 1)
@test wsum(x, w2, 2) ≈ sum(x .* w2', dims = 2)
@test wsum(x, w3, 3) ≈ sum(x .* reshape(w3, 1, 1, 4), dims = 3)

v = view(x, 2:4, :, :)

@test wsum(v, w1[1:3], 1) ≈ sum(v .* w1[1:3], dims = 1)
@test wsum(v, w2, 2)      ≈ sum(v .* w2', dims = 2)
@test wsum(v, w3, 3)      ≈ sum(v .* reshape(w3, 1, 1, 4), dims = 3)

## wsum for Arrays with non-BlasReal elements
x = rand(1:100, 6, 8)
w1 = rand(6)
w2 = rand(8)

@test wsum(x, w1, 1) ≈ sum(x .* w1, dims = 1)
@test wsum(x, w2, 2) ≈ sum(x .* w2', dims = 2)

## wsum!
x = rand(6)
w = rand(6)

r = ones(1)
@test wsum!(r, x, w, 1; init=true) === r
@test r ≈ [dot(x, w)]

r = ones(1)
@test wsum!(r, x, w, 1; init=false) === r
@test r ≈ [dot(x, w) + 1.0]

x = rand(6, 8)
w1 = rand(6)
w2 = rand(8)

r = ones(1, 8)
@test wsum!(r, x, w1, 1; init=true) === r
@test r ≈ sum(x .* w1, dims = 1)

r = ones(1, 8)
@test wsum!(r, x, w1, 1; init=false) === r
@test r ≈ sum(x .* w1, dims = 1) .+ 1.0

r = ones(6)
@test wsum!(r, x, w2, 2; init=true) === r
@test r ≈ sum(x .* w2', dims = 2)

r = ones(6)
@test wsum!(r, x, w2, 2; init=false) === r
@test r ≈ sum(x .* w2', dims = 2) .+ 1.0

x = rand(8, 6, 5)
w1 = rand(8)
w2 = rand(6)
w3 = rand(5)

r = ones(1, 6, 5)
@test wsum!(r, x, w1, 1; init=true) === r
@test r ≈ sum(x .* w1, dims = 1)

r = ones(1, 6, 5)
@test wsum!(r, x, w1, 1; init=false) === r
@test r ≈ sum(x .* w1, dims = 1) .+ 1.0

r = ones(8, 1, 5)
@test wsum!(r, x, w2, 2; init=true) === r
@test r ≈ sum(x .* w2', dims = 2)

r = ones(8, 1, 5)
@test wsum!(r, x, w2, 2; init=false) === r
@test r ≈ sum(x .* w2', dims = 2) .+ 1.0

r = ones(8, 6)
@test wsum!(r, x, w3, 3; init=true) === r
@test r ≈ sum(x .* reshape(w3, (1, 1, 5)), dims = 3)

r = ones(8, 6)
@test wsum!(r, x, w3, 3; init=false) === r
@test r ≈ sum(x .* reshape(w3, (1, 1, 5)), dims = 3) .+ 1.0

## the sum and mean syntax
a = reshape(1.0:27.0, 3, 3, 3)

@testset "Sum $f" for f in weight_funcs
    @test sum([1.0, 2.0, 3.0], f([1.0, 0.5, 0.5])) ≈ 3.5
    @test sum(1:3, f([1.0, 1.0, 0.5]))             ≈ 4.5

    for wt in ([1.0, 1.0, 1.0], [1.0, 0.2, 0.0], [0.2, 0.0, 1.0])
        @test sum(a, f(wt), 1)  ≈ sum(a.*reshape(wt, length(wt), 1, 1), dims = 1)
        @test sum(a, f(wt), 2)  ≈ sum(a.*reshape(wt, 1, length(wt), 1), dims = 2)
        @test sum(a, f(wt), 3)  ≈ sum(a.*reshape(wt, 1, 1, length(wt)), dims = 3)
    end
end

@testset "Mean $f" for f in weight_funcs
    @test mean([1:3;], f([1.0, 1.0, 0.5])) ≈ 1.8
    @test mean(1:3, f([1.0, 1.0, 0.5]))    ≈ 1.8

    for wt in ([1.0, 1.0, 1.0], [1.0, 0.2, 0.0], [0.2, 0.0, 1.0])
        @test mean(a, f(wt), dims=1) ≈ sum(a.*reshape(wt, length(wt), 1, 1), dims=1)/sum(wt)
        @test mean(a, f(wt), dims=2) ≈ sum(a.*reshape(wt, 1, length(wt), 1), dims=2)/sum(wt)
        @test mean(a, f(wt), dims=3) ≈ sum(a.*reshape(wt, 1, 1, length(wt)), dims=3)/sum(wt)
        @test_throws ErrorException mean(a, f(wt), dims=4)
    end
end

# Quantile fweights
@testset "Quantile fweights" begin
    data = (
        [7, 1, 2, 4, 10],
        [7, 1, 2, 4, 10],
        [7, 1, 2, 4, 10, 15],
        [1, 2, 4, 7, 10, 15],
        [0, 10, 20, 30],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
        [30, 40, 50, 60, 35],
        [2, 0.6, 1.3, 0.3, 0.3, 1.7, 0.7, 1.7],
        [1, 2, 2],
        [3.7, 3.3, 3.5, 2.8],
        [100, 125, 123, 60, 45, 56, 66],
        [2, 2, 2, 2, 2, 2],
        [2.3],
        [-2, -3, 1, 2, -10],
        [1, 2, 3, 4, 5],
        [5, 4, 3, 2, 1],
        [-2, 2, -1, 3, 6],
        [-10, 1, 1, -10, -10],
    )
    wt = (
        [3, 1, 1, 1, 3],
        [1, 1, 1, 1, 1],
        [3, 1, 1, 1, 3, 3],
        [1, 1, 1, 3, 3, 3],
        [30, 191, 9, 0],
        [10, 1, 1, 1, 9],
        [10, 1, 1, 1, 900],
        [1, 3, 5, 4, 2],
        [2, 2, 5, 0, 2, 2, 1, 6],
        [1, 1, 8],
        [5, 5, 4, 1],
        [30, 56, 144, 24, 55, 43, 67],
        [1, 2, 3, 4, 5, 6],
        [12],
        [7, 1, 1, 1, 6],
        [1, 0, 0, 0, 2],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 2, 1],
        [0, 1, 1, 1, 1],
    )
    p = [0.0, 0.25, 0.5, 0.75, 1.0]
    function _rep(x::AbstractVector, lengths::AbstractVector{Int})
        res = similar(x, sum(lengths))
        i = 1
        for idx in 1:length(x)
            tmp = x[idx]
            for kdx in 1:lengths[idx]
                res[i] = tmp
                i += 1
            end
        end
        return res
    end
    # quantile with fweights is the same as repeated vectors
    for i = 1:length(data)
        @test quantile(data[i], fweights(wt[i]), p) ≈ quantile(_rep(data[i], wt[i]), p)
    end
    # quantile with fweights = 1  is the same as quantile
    for i = 1:length(data)
        @test quantile(data[i], fweights(fill!(similar(wt[i]), 1)), p) ≈ quantile(data[i], p)
    end

    # Issue #313
    @test quantile([1, 2, 3, 4, 5], fweights([0,1,2,1,0]), p) ≈ quantile([2, 3, 3, 4], p)
    @test quantile([1, 2], fweights([1, 1]), 0.25) ≈ 1.25
    @test quantile([1, 2], fweights([2, 2]), 0.25) ≈ 1.0

    # test non integer frequency weights
    quantile([1, 2], fweights([1.0, 2.0]), 0.25) == quantile([1, 2], fweights([1, 2]), 0.25)
    @test_throws ArgumentError quantile([1, 2], fweights([1.5, 2.0]), 0.25)

    @test_throws ArgumentError quantile([1, 2], fweights([1, 2]), nextfloat(1.0))
    @test_throws ArgumentError quantile([1, 2], fweights([1, 2]), prevfloat(0.0))
end

@testset "Quantile aweights, pweights and weights" for f in (aweights, pweights, weights)
    data = (
        [7, 1, 2, 4, 10],
        [7, 1, 2, 4, 10],
        [7, 1, 2, 4, 10, 15],
        [1, 2, 4, 7, 10, 15],
        [0, 10, 20, 30],
        [1, 2, 3, 4, 5],
        [1, 2, 3, 4, 5],
        [30, 40, 50, 60, 35],
        [2, 0.6, 1.3, 0.3, 0.3, 1.7, 0.7, 1.7],
        [1, 2, 2],
        [3.7, 3.3, 3.5, 2.8],
        [100, 125, 123, 60, 45, 56, 66],
        [2, 2, 2, 2, 2, 2],
        [2.3],
        [-2, -3, 1, 2, -10],
        [1, 2, 3, 4, 5],
        [5, 4, 3, 2, 1],
        [-2, 2, -1, 3, 6],
        [-10, 1, 1, -10, -10],
    )
    wt = (
        [1, 1/3, 1/3, 1/3, 1],
        [1, 1, 1, 1, 1],
        [1, 1/3, 1/3, 1/3, 1, 1],
        [1/3, 1/3, 1/3, 1, 1, 1],
        [30, 191, 9, 0],
        [10, 1, 1, 1, 9],
        [10, 1, 1, 1, 900],
        [1, 3, 5, 4, 2],
        [2, 2, 5, 1, 2, 2, 1, 6],
        [0.1, 0.1, 0.8],
        [5, 5, 4, 1],
        [30, 56, 144, 24, 55, 43, 67],
        [0.1, 0.2, 0.3, 0.4, 0.5, 0.6],
        [12],
        [7, 1, 1, 1, 6],
        [1, 0, 0, 0, 2],
        [1, 2, 3, 4, 5],
        [0.1, 0.2, 0.3, 0.2, 0.1],
        [1, 1, 1, 1, 1],
    )
    quantile_answers = (
        [1.0, 4.0, 6.0, 8.0, 10.0],
        [1.0, 2.0, 4.0, 7.0, 10.0],
        [1.0, 4.75, 7.5, 10.4166667, 15.0],
        [1.0, 4.75, 7.5, 10.4166667, 15.0],
        [0.0, 2.6178010, 5.2356021, 7.8534031, 20.0],
        [1.0, 4.0, 4.3333333, 4.6666667, 5.0],
        [1.0, 4.2475, 4.4983333, 4.7491667, 5.0],
        [30.0, 37.5, 44.0, 51.25, 60.0],
        [0.3, 0.7, 1.3, 1.7, 2.0],
        [1.0, 2.0, 2.0, 2.0, 2.0],
        [2.8, 3.15, 3.4, 3.56, 3.7],
        [45.0, 62.149253, 102.875, 117.4097222, 125.0],
        [2.0, 2.0, 2.0, 2.0, 2.0],
        [2.3, 2.3, 2.3, 2.3, 2.3],
        [-10.0, -2.7857143, -2.4285714, -2.0714286, 2.0],
        [1.0, 2.0, 3.0, 4.0, 5.0],
        [1.0, 1.625, 2.3333333, 3.25, 5.0],
        [-2.0, -1.3333333, 0.5, 2.5, 6.0],
        [-10.0, -10.0, -10.0, 1.0, 1.0]
    )
    p = [0.0, 0.25, 0.5, 0.75, 1.0]

    Random.seed!(10)
    for i = 1:length(data)
        @test quantile(data[i], f(wt[i]), p) ≈ quantile_answers[i] atol = 1e-5
        for j = 1:10
            # order of p does not matter
            reorder = sortperm(rand(length(p)))
            @test quantile(data[i], f(wt[i]), p[reorder]) ≈ quantile_answers[i][reorder] atol = 1e-5
        end
        for j = 1:10
            # order of w does not matter
            reorder = sortperm(rand(length(data[i])))
            @test quantile(data[i][reorder], f(wt[i][reorder]), p) ≈ quantile_answers[i] atol = 1e-5
        end
    end
    # All equal weights corresponds to base quantile
    for v in (1, 2, 345)
        for i = 1:length(data)
            w = f(fill(v, length(data[i])))
            @test quantile(data[i], w, p) ≈ quantile(data[i], p) atol = 1e-5
            for j = 1:10
                prandom = rand(4)
                @test quantile(data[i], w,  prandom) ≈ quantile(data[i], prandom) atol = 1e-5
            end
        end
    end
    # test zeros are removed
    for i = 1:length(data)
        @test quantile(vcat(1.0, data[i]), f(vcat(0.0, wt[i])), p) ≈ quantile_answers[i] atol = 1e-5
    end
    # Syntax
    v = [7, 1, 2, 4, 10]
    w = [1, 1/3, 1/3, 1/3, 1]
    answer = 6.0
    @test quantile(data[1], f(w), 0.5)    ≈  answer atol = 1e-5
end

@testset "Median $f" for f in weight_funcs
    data = [4, 3, 2, 1]
    wt = [0, 0, 0, 0]
    @test_throws ArgumentError median(data, f(wt))
    @test_throws ArgumentError median(Float64[], f(Float64[]))
    wt = [1, 2, 3, 4, 5]
    @test_throws ArgumentError median(data, f(wt))
    if VERSION >= v"1.0"
        @test_throws MethodError median([4 3 2 1 0], f(wt))
        @test_throws MethodError median([[1 2] ; [4 5] ; [7 8] ; [10 11] ; [13 14]], f(wt))
    end
    data = [1, 3, 2, NaN, 2]
    @test isnan(median(data, f(wt)))
    wt = [1, 2, NaN, 4, 5]
    @test_throws ArgumentError median(data, f(wt))
    data = [1, 3, 2, 1, 2]
    @test_throws ArgumentError median(data, f(wt))
    wt = [-1, -1, -1, -1, -1]
    @test_throws ArgumentError median(data, f(wt))
    wt = [-1, -1, -1, 0, 0]
    @test_throws ArgumentError median(data, f(wt))

    data = [4, 3, 2, 1]
    wt = [1, 2, 3, 4]
    @test median(data, f(wt)) ≈ quantile(data, f(wt), 0.5) atol = 1e-5
end

@testset "Mismatched eltypes" begin
    @test round(mean(Union{Int,Missing}[1,2], weights([1,2])), digits=3) ≈ 1.667
end

@testset "Exponential Weights" begin
    @testset "Usage" begin
        θ = 5.25
        λ = 1 - exp(-1 / θ)     # simple conversion for the more common/readable method
        v = [λ*(1-λ)^(1-i) for i = 1:4]
        w = Weights(v)

        @test round.(w, digits=4) == [0.1734, 0.2098, 0.2539, 0.3071]

        @testset "basic" begin
            @test eweights(1:4, λ) ≈ w
        end

        @testset "1:n" begin
            @test eweights(4, λ) ≈ w
        end

        @testset "indexin" begin
            v = [λ*(1-λ)^(1-i) for i = 1:10]

            # Test that we should be able to skip indices easily
            @test eweights([1, 3, 5, 7], 1:10, λ) ≈ Weights(v[[1, 3, 5, 7]])

            # This should also work with actual time types
            t1 = DateTime(2019, 1, 1, 1)
            tx = t1 + Hour(7)
            tn = DateTime(2019, 1, 2, 1)

            @test eweights(t1:Hour(2):tx, t1:Hour(1):tn, λ) ≈ Weights(v[[1, 3, 5, 7]])
        end
    end

    @testset "Empty" begin
        @test eweights(0, 0.3) == Weights(Float64[])
        @test eweights(1:0, 0.3) == Weights(Float64[])
        @test eweights(Int[], 1:10, 0.4) == Weights(Float64[])
    end

    @testset "Failure Conditions" begin
        # λ > 1.0
        @test_throws ArgumentError eweights(1, 1.1)

        # time indices are not all positive non-zero integers
        @test_throws ArgumentError eweights([0, 1, 2, 3], 0.3)

        # Passing in an array of bools will work because Bool <: Integer,
        # but any `false` values will trigger the same argument error as 0.0
        @test_throws ArgumentError eweights([true, false, true, true], 0.3)
    end
end

end # @testset StatsBase.Weights


text:
struct OfflineFile end

struct KM3NETDAQSnapshotHit <: UnROOT.CustomROOTStruct
    dom_id::Int32
    channel_id::UInt8
    t::Int32
    tot::UInt8
end
function UnROOT.readtype(io, T::Type{KM3NETDAQSnapshotHit})
    T(UnROOT.readtype(io, Int32), read(io, UInt8), read(io, Int32), read(io, UInt8))
end
function UnROOT.interped_data(rawdata, rawoffsets, ::Type{Vector{KM3NETDAQSnapshotHit}}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, KM3NETDAQSnapshotHit, skipbytes=10)
end

struct KM3NETDAQTriggeredHit <: UnROOT.CustomROOTStruct
    dom_id::Int32
    channel_id::UInt8
    t::Int32
    tot::UInt8
    trigger_mask::UInt64
end
UnROOT.packedsizeof(::Type{KM3NETDAQTriggeredHit}) = 24  # incl. cnt and vers

function UnROOT.readtype(io, T::Type{KM3NETDAQTriggeredHit})
    dom_id = UnROOT.readtype(io, Int32)
    channel_id = read(io, UInt8)
    tdc = read(io, Int32)
    tot = read(io, UInt8)
    cnt = read(io, UInt32)
    vers = read(io, UInt16)
    trigger_mask = UnROOT.readtype(io, UInt64)
    T(dom_id, channel_id, tdc, tot, trigger_mask)
end

function UnROOT.interped_data(rawdata, rawoffsets, ::Type{Vector{KM3NETDAQTriggeredHit}}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, KM3NETDAQTriggeredHit, skipbytes=10)
end

struct KM3NETDAQEventHeader
    detector_id::Int32
    run::Int32
    frame_index::Int32
    UTC_seconds::UInt32
    UTC_16nanosecondcycles::UInt32
    trigger_counter::UInt64
    trigger_mask::UInt64
    overlays::UInt32
end
packedsizeof(::Type{KM3NETDAQEventHeader}) = 76

function UnROOT.readtype(io::IO, T::Type{KM3NETDAQEventHeader})
    skip(io, 18)
    detector_id = UnROOT.readtype(io, Int32)
    run = UnROOT.readtype(io, Int32)
    frame_index = UnROOT.readtype(io, Int32)
    skip(io, 6)
    UTC_seconds = UnROOT.readtype(io, UInt32)
    UTC_16nanosecondcycles = UnROOT.readtype(io, UInt32)
    skip(io, 6)
    trigger_counter = UnROOT.readtype(io, UInt64)
    skip(io, 6)
    trigger_mask = UnROOT.readtype(io, UInt64)
    overlays = UnROOT.readtype(io, UInt32)
    T(detector_id, run, frame_index, UTC_seconds, UTC_16nanosecondcycles, trigger_counter, trigger_mask, overlays)
end
function UnROOT.interped_data(rawdata, rawoffsets, ::Type{KM3NETDAQEventHeader}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, KM3NETDAQEventHeader, jagged=false)
end

struct OnlineEvent
    header::KM3NETDAQEventHeader
    snapshot_hits::Vector{KM3NETDAQSnapshotHit}
    triggered_hits::Vector{KM3NETDAQTriggeredHit}
end

struct OnlineFile
    fobj::UnROOT.ROOTFile
    headers::Vector{KM3NETDAQEventHeader}
    snapshot_hits
    triggered_hits

    function OnlineFile(filename::AbstractString)
        customstructs = Dict(
            "KM3NETDAQ::JDAQEvent.snapshotHits" => Vector{KM3NETDAQSnapshotHit},
            "KM3NETDAQ::JDAQEvent.triggeredHits" => Vector{KM3NETDAQTriggeredHit},
            "KM3NETDAQ::JDAQEvent.KM3NETDAQ::JDAQEventHeader" => KM3NETDAQEventHeader
        )
        fobj = UnROOT.ROOTFile(filename, customstructs=customstructs)

        new(fobj,
            fobj["KM3NET_EVENT/KM3NET_EVENT/KM3NETDAQ::JDAQEventHeader"],
            fobj["KM3NET_EVENT/KM3NET_EVENT/snapshotHits"],
            fobj["KM3NET_EVENT/KM3NET_EVENT/triggeredHits"])
    end
end
Base.getindex(f::OnlineFile, idx::Integer) = OnlineEvent(
    f.headers[idx], f.snapshot_hits[idx], f.triggered_hits[idx]
)
Base.length(f::OnlineFile) = length(f.headers)

Base.close(f::OnlineFile) = close(f.fobj)

function read_headers(f::OnlineFile)
    data, offsets = UnROOT.array(f.fobj, "KM3NET_EVENT/KM3NET_EVENT/KM3NETDAQ::JDAQEventHeader"; raw=true)
    UnROOT.splitup(data, offsets, KM3NETDAQEventHeader; jagged=false)
end

function read_snapshot_hits(f::OnlineFile)
    data, offsets = UnROOT.array(f.fobj, "KM3NET_EVENT/KM3NET_EVENT/snapshotHits"; raw=true)
    UnROOT.splitup(data, offsets, KM3NETDAQSnapshotHit, skipbytes=10)
end

function read_triggered_hits(f::OnlineFile)
    data, offsets = UnROOT.array(f.fobj, "KM3NET_EVENT/KM3NET_EVENT/triggeredHits"; raw=true)
    UnROOT.splitup(data, offsets, KM3NETDAQTriggeredHit, skipbytes=10)
end

# Triggers
is3dmuon(e::DAQEvent) = nthbitset(Trigger.JTRIGGER3DMUON, e.trigger_mask)
is3dshower(e::DAQEvent) = nthbitset(Trigger.JTRIGGER3DSHOWER, e.trigger_mask)
ismxshower(e::DAQEvent) = nthbitset(Trigger.JTRIGGERMXSHOWER, e.trigger_mask)
isnb(e::DAQEvent) = nthbitset(Trigger.JTRIGGERNB, e.trigger_mask)
@deprecate is_3dmuon is3dmuon
@deprecate is_3dshower is3dshower
@deprecate is_mxshower ismxshower

# Online DAQ readout

function Base.read(s::IO, ::Type{T}; legacy=false) where T<:DAQEvent
    length = read(s, Int32)
    type = read(s, Int32)
    version = Int16(0)
    !legacy && (version = read(s, Int16))
    det_id = read(s, Int32)
    run_id = read(s, Int32)
    timeslice_id = read(s, Int32)
    timestamp = read(s, Int32)
    ticks = read(s, Int32)
    trigger_counter = read(s, Int64)
    trigger_mask = read(s, Int64)
    overlays = read(s, Int32)

    n_triggered_hits = read(s, Int32)
    triggered_hits = Vector{TriggeredHit}()
    sizehint!(triggered_hits, n_triggered_hits)
    triggered_map = Dict{Tuple{Int32, UInt8, Int32, UInt8}, Int64}()
    @inbounds for i ∈ 1:n_triggered_hits
        dom_id = read(s, Int32)
        channel_id = read(s, UInt8)
        time = bswap(read(s, Int32))
        tot = read(s, UInt8)
        trigger_mask = read(s, Int64)
        triggered_map[(dom_id, channel_id, time, tot)] = trigger_mask
        push!(triggered_hits, TriggeredHit(dom_id, channel_id, time, tot, trigger_mask))
    end

    n_hits = read(s, Int32)
    hits = Vector{Hit}()
    sizehint!(hits, n_hits)
    @inbounds for i ∈ 1:n_hits
        dom_id = read(s, Int32)
        channel_id = read(s, UInt8)
        time = bswap(read(s, Int32))
        tot = read(s, UInt8)
        key = (dom_id, channel_id, time, tot)
        triggered = false
        if haskey(triggered_map, key)
            triggered = true
        end
        push!(hits, Hit(channel_id, dom_id, time, tot, triggered))
    end

    T(det_id, run_id, timeslice_id, timestamp, ticks, trigger_counter, trigger_mask, overlays, n_triggered_hits, triggered_hits, n_hits, hits)
end


text:
# wrappers

function cudnnSoftmaxForward(x::DenseCuArray{T,4}, y::DenseCuArray{T,4}=x;
                             algo=CUDNN_SOFTMAX_FAST, # or CUDNN_SOFTMAX_ACCURATE
                             mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                             alpha=1.0, beta=0.0) where T
    cudnnSoftmaxForward(handle(), algo, mode,
                        scalingParameter(T, alpha), TensorDesc(x), x,
                        scalingParameter(T, beta ), TensorDesc(y), y)
    return y
end

function cudnnSoftmaxBackward(y::DenseCuArray{T,4}, dy::DenseCuArray{T,4}, dx::DenseCuArray{T,4}=dy;
                              algo=CUDNN_SOFTMAX_FAST, # or CUDNN_SOFTMAX_ACCURATE
                              mode=CUDNN_SOFTMAX_MODE_INSTANCE, # or CUDNN_SOFTMAX_MODE_CHANNEL
                              alpha=1.0, beta=0.0) where T
    cudnnSoftmaxBackward(handle(), algo, mode,
                         scalingParameter(T, alpha), TensorDesc(y), y,
                         TensorDesc(dy), dy,
                         scalingParameter(T, beta ), TensorDesc(dx), dx)
    return dx
end


text:
module RemoteREPL

export connect_repl, serve_repl, @remote

const DEFAULT_PORT = 27754
const PROTOCOL_MAGIC = "RemoteREPL"
const PROTOCOL_VERSION = UInt32(1)

const STDOUT_PLACEHOLDER = Symbol("#RemoteREPL_STDOUT_PLACEHOLDER")

include("tunnels.jl")
include("server.jl")
include("client.jl")

end


text:
@testset "206.reverse-linked-list.jl" begin
    @test reverse_list(ListNode{Int}([1, 2, 3, 4, 5])) == ListNode{Int}([5, 4, 3, 2, 1])
    @test reverse_list(ListNode{Int}([1, 2])) == ListNode{Int}([2, 1])
    @test reverse_list(nothing) === nothing
end


text:
"""
测试单带三角格子，tJ
"""

using Plots
using MARY_fRG.Interactions
using MARY_fRG.Fermi, MARY_fRG.Fermi.Surface
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

function run()
    model = common_triangle_lattice(-2.0)
    Γ4 = ECGamma4(
        model, 6.85, 18, 50
    )
    Γ4.V[1, 1, 1, 1, :, :, :] .+= triangular_system_heisenberg(
        Γ4.patches[1], Γ4.model.kadd, 1.0
    )
    lval = 0.
    lstep = 0.01
    for idx = 1:1:301
        bubb_pp, bubb_fs, bubb_nfs, bubb_ex, bubb_nex = all_bubble_ec_mt(Γ4, lval)
        dl = dl_ec_mt(Γ4, bubb_pp, bubb_fs, bubb_nfs, bubb_ex, bubb_nex)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma4"*string(idx))
        end
        lval += lstep
    end
end

run()


text:
# patch for ColorVectorSpace 0.9
# for CVS < 0.9, we can just use the fallback solution in Distances
if isdefined(ImageCore.ColorVectorSpace, :⊙)
    # Because how abs2 calculated in color vector space is ambiguious, abs2(::RGB) is un-defined
    # since ColorVectorSpace 0.9
    # https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/131
    _abs2(c::Colorant) = mapreducec(abs2, +, 0, c)
    _abs2(x::Number) = abs2(x)

    _abs(c::Colorant) = mapreducec(abs, +, 0, c)
    _abs(x::Number) = abs(x)
else
    _abs2(c::Union{Colorant, Number}) = abs2(c)
    _abs(x::Number) = abs(x)
end