text:
using .JSort
import Base: read, read!, show, length, sum
using DataFrames
import PyPlot; const plt = PyPlot
using LaTeXStrings

macro tostring(arg)
    string(arg)
end

abstract type AbstractDetector end
abstract type AbstractDetector1D <: AbstractDetector end
abstract type AbstractDetector2D <: AbstractDetector end

length(d::AbstractDetector1D) = length(d.ref[])

function read(path::AbstractString, V::Type{W}; T=UInt32) where W <: AbstractDetector
    detector = open(path, "r") do io
        read(io, V, T=T)
    end
end

function read(io::IO, V::Type{<:AbstractDetector1D}; T=UInt32)
    len = stat(io).size / sizeof(T) |> Int
    data = Vector{Float32}(undef, len)
    for i in 1:len
        data[i] = read(io, T) |> Float32
    end
    pattern = r".*chn(?<c>\d+).*"
    m = match(pattern, io.name)
    chn = if isnothing(m) 0 else parse(Int, m[:c]) end
    V(data, chn)
end

function read(io::IO, V::Type{W} where W <: AbstractDetector2D; T=UInt32)
    # Assumes binary format int UInt32
    len = stat(io).size / 2sizeof(T) |> Int
    data = Matrix{Float32}(undef, len, 2)
    for i in 1:len
        data[i, 1] = read(io, T) |> Float32
        data[i, 2] = read(io, T) |> Float32
    end

    pattern = r".*b(?<b>\d)f(?<f>\d).*"
    m = match(pattern, io.name)
    if isnothing(m)
        pattern = r".*f(?<f>\d).*"
        m = match(pattern, io.name)
        if isnothing(m)
            f = b = 0
        else
            f = parse(Int, m[:f])
            b = 0
        end
    else
        f, b = parse.(Int, (m[:f], m[:b]))
    end
    V(data[:, 1], data[:, 2], f, b)
end

function write(path::AbstractString, V::AbstractDetector)
    open(path, "w") do io
        write(io, V)
    end
end

function write(io::IO, V::AbstractDetector1D)
    write(io, V.ref[])
end

function write(io::IO, V::AbstractDetector2D)
    for i in eachindex(V.Xref[])
        write(io, V.Xref[][i], V.Yref[][i])
    end
end

function histogram(X::AbstractDetector2D; nbins=1000)
    edges = histogram_edges(X.Xref[], X.Yref[], nbins=nbins)
    histogram(X, edges...), edges
end

function histogram(X::AbstractDetector2D, xedges, yedges)
    #histogram(X.Xref[], X.Yref[], xedges, yedges)
    histogram(X.Xref[] .+ rand(length(X.Xref[])),
              X.Yref[] .+ rand(length(X.Yref[])),
              xedges, yedges)
end

function histogram(X::AbstractDetector2D, edges)
    histogram(X.Xref[], X.Yref[], edges...)
end

function histogram(X::AbstractDetector1D, edges::AbstractRange)
    # Why do I need to add a random number to remove histogram artifacts?
    histogram(X.ref[] .+ rand(length(X.ref[])), edges)
    #histogram(X.ref[], edges)
end
function histogram(X::AbstractDetector1D; nbins=1000)
    edges = range(minimum(X.ref[]), stop=maximum(X.ref[]), length=nbins)
    histogram(X, edges), edges
end

function (calib::AbstractCalibrator1D)(d::T; inplace::Bool=true)::T where {T<:AbstractDetector1D}
    poly = calib.poly[d.channel]
    if inplace
        d.ref[] .= poly.(d.ref[])
        return d
    else
        d = deepcopy(d)
        d.ref[] .= poly.(d.ref[])
        return d
    end
end

function show(io::IO, V::AbstractDetector1D)
    println("Channel: $(V.channel)")
    show(io, DataFrame(e=V.ref[]), allrows=false)
end
## Concrete implementations

struct ParticleDetector <: AbstractDetector2D
    e::Vector{Float32}
    Δe::Vector{Float32}
    front::Int
    back::Int
    Xref::Ref{Vector{Float32}}
    Yref::Ref{Vector{Float32}}
    function ParticleDetector(a::AbstractVector, b::AbstractVector, c::Number, d::Number)
        new(a, b, c, d, Ref(a), Ref(b))
    end
end
ParticleDetector(path::AbstractString; T=UInt32) = read(path, ParticleDetector, T=T)

# function write(io::IO, d::ParticleDetector)
#     for i in 1:length(d.e)
#         write(io, d.e[i], d.Δe[i])
#     end
# end

function show(io::IO, detector::AbstractDetector2D)
    println("Front: $(detector.front)\nBack: $(detector.back)")
    show(DataFrame(e=detector.Xref[], Δe=detector.Yref[]), allrows=false)
end

function (calib::ParticleCalibrator)(d::ParticleDetector)
    poly = calib.poly[d.front, d.back]
    if poly.var == :e
        map!(poly, d.e, d.e)
    elseif poly.var == :Δe
        map!(poly, d.Δe, d.Δe)
    else
        throw("Unable to deduce how to apply calibration. Specify `e` or `Δe`")
    end
end

function (calib::ParticleCalibrator)(d::ParticleDetector, var::Symbol)
    poly = calib.poly[d.front, d.back]
    if var == :e
        map!(poly, d.e, d.e)
    elseif var == :Δe
        map!(poly, d.Δe, d.Δe)
    else
        throw("Variable must be `e` or `Δe`, not $var")
    end
end

function calibrateall(path::AbstractString; epath="coefficients_e.txt", Δepath="coefficients_de.txt")
    j(x) = joinpath(path, x)
    ce = ParticleCalibrator(epath, var=:e)
    cde = ParticleCalibrator(Δepath, var=:Δe)

    for f in 1:8
        totale = Float32[]
        totalde = Float32[]
        for b in 1:8
            detector = ParticleDetector(j("edeb$(b)f$(f).bin"))
            ce(detector)
            cde(detector)
            append!(totale, detector.e)
            append!(totalde, detector.Δe)
        end
        write(j("edef$(f).bin"), ParticleDetector(totale, totalde, f, 0))
    end
end

function plot(d::ParticleDetector; ax=nothing)
    (fig, ax) = if isnothing(ax)
        plt.subplots()
    else
        ax.figure, ax
    end
    mat, edges = histogram(d, nbins=1000)
    #ax.scatter(d.e, d.Δe, marker=",", s=1, alpha=0.5)
    ax.pcolormesh(edges[1], edges[2], mat |> transpose .|> log10)
    ax.set_xlabel(L"E\quad [keV?]")
    ax.set_ylabel(L"\Delta E\quad [keV?]")
    ax.set_title("f$(d.front)b$(d.back)")
end

### Gamma Detector

struct GammaDetector <: AbstractDetector1D
    e::Vector{Float32}
    channel::Int
    ref::Ref{Vector{Float32}}
    GammaDetector(e::AbstractVector) = new(e, 0, Ref(e))
    GammaDetector(e::AbstractVector, chn::Int) = new(e, chn, Ref(e))
end
GammaDetector(path::AbstractString; T=Float32) = read(path, GammaDetector, T=T)
GammaDetector() = GammaDetector(Float32[])

function plot(D::AbstractDetector1D; ax=nothing, nbins=1000)
    (fig, ax) = if isnothing(ax)
        plt.subplots()
    else
        ax.figure, ax
    end
    hist, edges = histogram(D, nbins=nbins)
    ax.step(edges, hist)
    ax.set_yscale("log")
    ax.set_title("Channel $(D.channel)")
end

### Time 'Detector'

struct TimeDetector <: AbstractDetector1D
    t::Vector{Float32}
    channel::Int
    ref::Ref{Vector{Float32}}
    TimeDetector(t::AbstractVector) = new(t, 0, Ref(t))
    TimeDetector(t::AbstractVector, chn::Int) = new(t, chn, Ref(t))
end
TimeDetector(path::AbstractString; T=Float32) = read(path, TimeDetector, T=T)
TimeDetector() = TimeDetector(Float32[])

### Channel 'Detector'
struct ChannelDetector{T} <: AbstractDetector
    data::Vector{T}
end

function ChannelDetector(name::AbstractString, path="."; T=Float32)
    pattern = Regex("^" * name * "_chn(?<chn>\\d+).*")
    ChannelDetector(pattern, path, T=T)
end

function ChannelDetector(pattern::Regex, path="."; T=Float32)
    data = nothing
    for fname in readdir(path)
        m = match(pattern, fname)
        if !isnothing(m)
            name = splitext(basename(fname))[1]
            types = parsetype(name)
            Tcalib = types[1]
            if isnothing(data)
                data = Tcalib[Tcalib() for i in 1:33]
            end
            chn = parse(Int, m[:chn])
            data[chn] = Tcalib(joinpath(path, fname), T=T)
        end
    end
    ChannelDetector(data)
end

function show(io::IO, d::ChannelDetector)
    for chn in eachindex(d.data)
        println(io, "Channel $chn:\t$(length(d.data[chn]))")
    end
end

function plot(D::ChannelDetector; ax=nothing, nbins=1000)
    (fig, ax) = if isnothing(ax)
        plt.subplots()
    else
        ax.figure, ax
    end
    hists = Vector{Int}[]
    edges = histogram_edges(D.data[1].ref[], nbins=nbins)
    chns = goodchannels(D)
    for chn in chns
        hist = histogram(D.data[chn], edges)
        push!(hists, hist)
    end
    hists = cat(hists'..., dims=1)
    ax.pcolormesh(edges, 1:length(chns)+1, hists .|> log10)
    #ax.matshow(hists .|> log10)
    #ax.pcolormesh(1:length(edges), 1:length(chns), hists .|> log10)
    ax.set_yticks(1.5:length(chns)+1-0.5)
    ax.set_yticklabels(chns)
    ax
end

function goodchannels(D::ChannelDetector; deviationlimit=10)
    good_lengths = [chn for (chn, d) in enumerate(D.data) if length(d) > 100 && chn != 32 && chn != 30]
    lengths = length.(D.data[chn] for chn in good_lengths)
    means = mean.(D.data[chn].ref[] for chn in good_lengths)
    M = mean(means)
    deviation = @. 1/(lengths - 1)*√((means - M)^2)
    good_lengths[deviation .< deviationlimit]
end

function histogram(D::ChannelDetector; nbins=1000)
    chns = goodchannels(D)
    edges = histogram_edges(D.data[1].ref[], nbins=nbins)
    [histogram(D.data[chn], edges) for chn in chns], edges
end

function histogram(D::ChannelDetector, edges::AbstractRange)
    [histogram(D.data[chn], edges) for chn in goodchannels(D)]
end

function align(detector::ChannelDetector; nbins=1000, kwargs...)
    d = deepcopy(detector)
    calib, invcalib = align!(d, nbins=nbins; kwargs...)
    d, calib, invcalib
end

function align!(detector::ChannelDetector{GammaDetector}; nbins=1000, kwargs...)
    signal = pyimport("scipy.signal")
    smoother(x) = signal.savgol_filter(x, 51, 4)

    chns = goodchannels(detector)
    hist, edges = histogram(detector, nbins=nbins)
    coeff, invcoeff = alignspectra(hist, edges; highsmoother=smoother, kwargs...)
    def = if (length(coeff[1]) == 2)
        [0.0 1.0] else [0.0 1.0 0.0] end
    coefficients = repeat(def, outer=(33, 1))
    invcoefficients = copy(coefficients)

    coefficients[chns, :] .= hcat(coeff...)'
    calib = GammaCalibrator(coefficients)
    calib(detector)
    invcalib = GammaCalibrator(coefficients)
    calib, invcalib
end

function align!(detector::ChannelDetector{TimeDetector}; nbins=1000, kwargs...)
    chns = goodchannels(detector)
    hist, edges = histogram(detector, nbins=nbins)
    coeff = aligntime(hist, edges; kwargs...)
    coefficients = repeat([0.0 1.0], outer=(33, 1))
    coefficients[chns, :] .= hcat(coeff...)'

    calib = TimeCalibrator(coefficients)
    calib(detector)
    calib, calib
end

function (c::AbstractCalibrator1D)(d::ChannelDetector; inplace::Bool=true)::ChannelDetector
    if inplace
        map(c, d.data)
        return d
    else
        d = deepcopy(d)
        map(c, d.data)
        return d
    end
end

function sum(d::ChannelDetector{T}) where T
    T(cat([det.ref[] for det in d.data]...; dims=1))
end

getindex(X::ChannelDetector, i) = getindex(X.data, i)
getindex(X::ChannelDetector, I...) = getindex(X.data, I...)
setindex!(X::ChannelDetector, i) = setindex!(X.data, i)
setindex!(X::ChannelDetector, I...) = setindex!(X.data, I...)

### General loading
struct Detector <: AbstractDetector
    data::Array{Float32}
    meta::Dict{Symbol, Int}
end
function Detector(path::AbstractString; T=UInt32)
    data, meta = readres(path, T=T)
    meta = Dict(i => j for (i, j) in meta if !isnothing(j))
    Detector(data, meta)
end

function plot(D::Detector; ax=nothing, nbins=1000)
    if ndims(D.data) == 1
        plot1D(D, ax=ax, nbins=nbins)
    elseif ndims(D.data) == 2
        plot2D(D, ax=ax, nbins=nbins)
    end
end

function plot1D(D::Detector; ax=nothing, nbins=1000)
    (fig, ax) = if isnothing(ax)
        plt.subplots()
    else
        ax.figure, ax
    end
    hist, edges = histogram(D.data, nbins=nbins)
    ax.step(edges, hist)
    ax.set_title(maketitle(D.meta))
end

function plot2D(D::Detector; ax=nothing, nbins=1000)
    (fig, ax) = if isnothing(ax)
        plt.subplots()
    else
        ax.figure, ax
    end
    mat, edges = histogram(D.data, nbins=nbins)
    ax.pcolormesh(edges[1], edges[2], mat |> transpose .|> log10)
    ax.set_title(maketitle(D.meta))
end

function readres(path::AbstractString; T=UInt32)
    open(path, "r") do io
        ndim = filedimensions(io.name)
        data = if ndim == 1
            readres1D(io, T=T)
        elseif ndim == 2
            readres2D(io, T=T)
        else
            error("Impossible dimension $ndim")
        end
        data, parsefilename(path)
    end
end

function readres1D(io::IO; T=UInt32)::Vector{Float32}
    len = stat(io).size / sizeof(T) |> Int
    data = Vector{Float32}(undef, len)
    for i in 1:len
        data[i] = read(io, T) |> Float32
    end
    data
end

function readres2D(io::IO; T=UInt32)::Matrix{Float32}
    len = stat(io).size / 2sizeof(T) |> Int
    data = Matrix{Float32}(undef, len, 2)
    for i in 1:len
        data[i, 1] = read(io, T) |> Float32
        data[i, 2] = read(io, T) |> Float32
    end
    data
end

function parsefilename(fname::AbstractString)
    fb_pattern = r".*b(?<b>\d)f(?<f>\d).*"
    chn_pattern = r".*chn(?<chn>\d+).*"
    f_pattern = r".*f(?<f>\d).*"
    m = match(fb_pattern, fname)
    f = b = chn = nothing
    if isnothing(m)
        m = match(f_pattern, fname)
        if !isnothing(m)
            f = parse(Int, m[:f])
        end
    else
        f, b = parse.(Int, (m[:f], m[:b]))
    end
    m = match(chn_pattern, fname)
    if !isnothing(m)
        chn = parse(Int, m[:chn])
    end
    (:f => f, :b => b, :chn => chn)
end

function filedimensions(fname::AbstractString)::Int
    if occursin("edeb", fname)
        return 2
    end
    underscores = count("_", fname)
    all = if occursin(r"\d", fname) 1 else 0 end
    underscores - all + 1
end

function maketitle(dict::Dict{Symbol, Int})
    @show dict
    title = ""
    k = keys(dict)
    if :f in k
        title *= "f$(dict[:f])"
    end
    if :b in k
        title *= " b$(dict[:b])"
    end
    if :chn in k
        title *= "Channel $(dict[:chn])"
    end
    title
end

function parsetype(fname::AbstractString)::Vector{DataType}
    tokens = split(fname, "_")
    types = DataType[]
    for token in tokens
        if token == "eg"
            push!(types, GammaDetector)
        elseif token == "t" || token == "tc"
            push!(types, TimeDetector)
        end
    end
    types
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
let
    x, y = variables(["x", "y"])

    #@test (x*y).alpha == [1 1]
    #@test (x*y).c == [1]

    #@test (1 + 2*x*y).alpha == [0 0; 1 1]
    #@test (1 + 2*x*y).c == [1; 2]

    #@test ( (x+y)*(x-y) ).alpha == [2 0; 0 2]
    #@test ( (x+y)*(x-y) ).c == [1; -1]

    #@test ( (1+x+y)^2 * y ).deg == 3
    #@test ( (1+x+y)^2 * y ).alpha == [0 1; 1 1; 0 2; 2 1; 1 2; 0 3]
    #@test ( (1+x+y)^2 * y ).c == [1; 2; 2; 1; 2; 1]

    @test ( x*(x^2-y) == x^3 - x*y)
    @test ( (x^2-y)*x == x^3 - x*y)
    @test ( (x^2-y)*(x+y+x*y) == x^3 + x^2*y + x^3*y - x*y - y^2 - x*y^2)
    @test ( (x+y+x*y)*(x^2-y) == x^3 + x^2*y + x^3*y - x*y - y^2 - x*y^2)

    @test ( [x, y]*[x, y]' == [ x^2 x*y; x*y y^2] )

    @test ( [1 0; 0 2]*[x, y] == [x, 2*y])
    @test ( dot([x,y], [x,y]) == 1.0*x^2 + y^2 )

    @test ( dot([1//2*x, y, 1], [2, 1, 1]) == x + y + 1 )

    @test ( x + (Polyopt.Poly(0) + Polyopt.Poly(0)) == x )
end

let
    x, = variables(["x"])
    @test ( (1+x)^2 == 1+2*x+x^2 )
end

let
    x, y, z = variables(["x", "y", "z"])
    @test Polyopt.monomials(3,[x,y,z]) == sort([1,x,y,z,x^2,x*y,x*z,y^2,y*z,z^2,x^3,x^2*y,x^2*z,x*y^2,x*y*z,x*z^2,y^3,y^2*z,y*z^2,z^3])
end

ext: jl


text:
@testset "Statistics" begin

    @testset "Statistics: Pearson Correlation" begin
        a = 1:10
        b = 1:10
        @test pearson_correlation(a, b) == 1

        x = [12,15,17,21,18]
        y = [6,9,-1,3,1]

        @test pearson_correlation(x, y) ==  -0.5201361507569889
    end

    @testset "Statistics: Variance" begin
        a = 1:10

        @test variance(a) == 9.166666666666666

        a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

        @test variance(a) == 9.166666666666666

        a = [12, 15, 17, 20]

        @test variance(a) == 11.333333333333334

    end

end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
import PyCall
import POMDPs: MDP, Policy

#
# Atari environment used to play Open AI Classic Control games: https://gym.openai.com/envs/#classic_control
#
# PS. It also supports multi-dimensional states, but it wont solve them efficiently
# as they will get converted to 1 dimensional array and will be used in a fully
# connected layer.
#
struct ClassicControlMDP <: MDP{Vector{Float64}, Int64}
    env::PyCall.PyObject
    state_space_size::Int64
    action_space::Vector{Int64}
end

#
# Initializer function that is using PyObject with Classic Control environment
# to initialize the variables
#
function ClassicControlMDP(env::PyCall.PyObject)

    states = prod(env[:observation_space][:shape])
    actions = Vector{Int64}(0:(env[:action_space][:n] - 1))

    return ClassicControlMDP(env, states, actions)
end

#
# We define number of function required to play the game
#
states(mdp::ClassicControlMDP) = mdp.state_space_size
actions(mdp::ClassicControlMDP) = mdp.action_space
n_states(mdp::ClassicControlMDP) = mdp.state_space_size
n_actions(mdp::ClassicControlMDP) = length(mdp.action_space)
initial_state(mdp::ClassicControlMDP, rng::AbstractRNG = MersenneTwister(0)) = mdp.env[:reset]()[:] .* 1.0
isterminal(mdp::ClassicControlMDP, s::Vector{Float64}) = iszero(s)
action_index(mdp::ClassicControlMDP, a::Int64) = a

#
# The only function we can use from OpenAIGym is getting a reward when taking a
# particular action in a specific step
#
function generate_sr(mdp::ClassicControlMDP, s::Vector{Float64}, a::Int64)

    state, reward, done = mdp.env[:step](a)

    if done
        return (state[:] .* 0., reward)
    end

    return (state[:] * 1., reward)
end

# Play a single episode
function run_experiment(env::ClassicControlMDP, policy::Policy; max_frame_iterations::Int64 = 100, keep_history = false)

    reward::Float64 = .0
    previous_state::Vector{Float64} = initial_state(env)

    #
    # We might decide to keep the history for some of our learning algorithms
    # therefore we need to initialize lists of a fixed size
    #
    list_size = if keep_history max_frame_iterations else zero(max_frame_iterations) end

    states = zeros(Float64, env.state_space_size, list_size)
    new_states = zeros(Float64, env.state_space_size, list_size)
    actions = zeros(Int64, list_size)
    rewards = zeros(Float64, list_size)

    frames_played = 0

    for i=1:max_frame_iterations

        current_action_i = action(policy, previous_state)
        current_action = action_index(env, current_action_i)
        current_state, current_reward = generate_sr(env, previous_state, current_action)
        reward += current_reward
        frames_played += 1

        #
        # In case we keep history lets track it
        #
        if keep_history
            states[:, i] = previous_state
            new_states[:, i] = current_state
            actions[i] = current_action_i
            rewards[i] = current_reward
        end

        if isterminal(env, current_state)
            break
        end

        previous_state = current_state

    end

    #
    # In case we keep history we should also cut in based on a number of frames required
    #
    if keep_history
        return states[:, 1:frames_played], actions[1:frames_played], reward, new_states[1:frames_played], rewards[1:frames_played]
    else
        return states, actions, reward, new_states
    end
end

#
# Executes policy on an environment for a number of repeats with a
# limitation on how many actions/ frames the episode can handle
#
function execute_policy(env::ClassicControlMDP, policy::Policy, experiment_repeats::Int64, max_frame_iterations::Int64; experiment_index = nothing, verbose = false)

    if experiment_index != nothing && verbose == true
        println("Policy execution #$experiment_index")
    end

    return mean(map(x -> run_experiment(env, policy, max_frame_iterations)[3], zeros(experiment_repeats)))
end

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
using OrdinaryDiffEq, DiffEqParamEstim, Test, NLopt, BlackBoxOptim
ms_f = function (du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3.0 * u[2] + u[1]*u[2]
end
ms_u0 = [1.0;1.0]
tspan = (0.0,10.0)
ms_p = [1.5,1.0]
ms_prob = ODEProblem(ms_f,ms_u0,tspan,ms_p)
t = collect(range(0, stop=10, length=200))
data = Array(solve(ms_prob,Tsit5(),saveat=t,abstol=1e-12,reltol=1e-12))
bound = Tuple{Float64, Float64}[(0, 10),(0, 10),(0, 10),(0, 10),
                                (0, 10),(0, 10),(0, 10),(0, 10),
                                (0, 10),(0, 10),(0, 10),(0, 10),
                                (0, 10),(0, 10),(0, 10),(0, 10),(0, 10),(0, 10)]

ms_obj = multiple_shooting_objective(ms_prob,Tsit5(),L2Loss(t,data);discontinuity_weight=1.0,abstol=1e-12,reltol=1e-12)
result = bboptimize(ms_obj;SearchRange = bound, MaxSteps = 21e3)
@test result.archive_output.best_candidate[end-1:end] ≈ [1.5,1.0] atol = 2e-1

priors = [Truncated(Normal(1.5,0.1),0,2),Truncated(Normal(1.0,0.1),0,1.5)]
ms_obj1 = multiple_shooting_objective(ms_prob,Tsit5(),L2Loss(t,data);priors=priors,discontinuity_weight=1.0,abstol=1e-12,reltol=1e-12)
result = bboptimize(ms_obj1;SearchRange = bound, MaxSteps = 21e3)
@test result.archive_output.best_candidate[end-1:end] ≈ [1.5,1.0] atol = 2e-1


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

JLLWrappers.@generate_main_file_header("HiGHS")
JLLWrappers.@generate_main_file("HiGHS", UUID("8fd58aa0-07eb-5a78-9b36-339c94fd15ea"))
end  # module HiGHS_jll


text:
immutable TransposedVector{T,V<:AbstractVector} <: AbstractMatrix{T}
    vec::V
    function TransposedVector(v::V)
        check_types(T,v)
        new(v)
    end
end

@pure check_types{T1,T2}(::Type{T1},::AbstractVector{T2}) = transpose_type(T1) === transpose_type(T2) ? nothing : error("Element type mismatch. Tried to create a `TransposedVector{$T}` from a $(typeof(v))")

# The element type is transformed as transpose is recursive
@inline transpose_type{T}(::Type{T}) = promote_op(transpose, T)

# Constructors that take a vector
@inline TransposedVector{T}(vec::AbstractVector{T}) = TransposedVector{transpose_type(T),typeof(vec)}(vec)
@inline (::Type{TransposedVector{T}}){T}(vec::AbstractVector{T}) = TransposedVector{T,typeof(vec)}(vec)

# Constructors that take a size and default to Array
@inline (::Type{TransposedVector{T}}){T}(n::Int) = TransposedVector{T}(Vector{transpose_type(T)}(n))
@inline (::Type{TransposedVector{T}}){T}(n1::Int, n2::Int) = n1 == 1 ? TransposedVector{T}(Vector{transpose_type(T)}(n2)) : error("TransposedVector expects 1×N size, got ($n1,$n2)")
@inline (::Type{TransposedVector{T}}){T}(n::Tuple{Int}) = TransposedVector{T}(Vector{transpose_type(T)}(n[1]))
@inline (::Type{TransposedVector{T}}){T}(n::Tuple{Int,Int}) = n[1] == 1 ? TransposedVector{T}(Vector{transpose_type(T)}(n[2])) : error("TransposedVector expects 1×N size, got $n")

# similar()
@inline similar(tvec::TransposedVector) = TransposedVector(similar(tvec.vec))
@inline similar{T}(tvec::TransposedVector, ::Type{T}) = TransposedVector(similar(tvec.vec, transpose_type(T)))
# There is no resizing similar() because it would be ambiguous if the result were a Matrix or a TransposedVector

# Basic methods
@inline transpose(vec::AbstractVector) = TransposedVector(vec)
@inline ctranspose{T}(vec::AbstractVector{T}) = TransposedVector(conj(vec))
@inline ctranspose{T<:Real}(vec::AbstractVector{T}) = TransposedVector(vec)

@inline transpose(tvec::TransposedVector) = tvec.vec
@inline ctranspose{T}(tvec::TransposedVector{T}) = conj(tvec.vec)
@inline ctranspose{T<:Real}(tvec::TransposedVector{T}) = tvec.vec

# Some overloads from Base
@inline transpose(r::Range) = TransposedVector(r)
@inline transpose(bv::BitVector) = TransposedVector(bv)
@inline transpose(sv::SparseVector) = TransposedVector(sv)

@inline ctranspose(r::Range) = TransposedVector(conj(r)) # is there such a thing as a complex range?
@inline ctranspose(bv::BitVector) = TransposedVector(bv)
@inline ctranspose(sv::SparseVector) = TransposedVector(conj(sv))

# Strictly, these are unnecessary but will make things stabler if we introduce
# a "view" for conj(::AbstractArray)
@inline conj(tvec::TransposedVector) = TransposedVector(conj(tvec.vec))
@inline conj{T<:Real}(tvec::TransposedVector{T}) = tvec

# AbstractArray interface
@inline length(tvec::TransposedVector) =  length(tvec.vec)
@inline size(tvec::TransposedVector) = (1, length(tvec.vec))
@inline size(tvec::TransposedVector, d) = ifelse(d==2, length(tvec.vec), 1)
linearindexing{V<:TransposedVector}(::Union{V,Type{V}}) = LinearFast()

@propagate_inbounds getindex(tvec::TransposedVector, i) = transpose(tvec.vec[i])
@propagate_inbounds setindex!(tvec::TransposedVector, v, i) = setindex!(tvec.vec, transpose(v), i)

# Cartesian indexing is distorted by getindex
# Furthermore, Cartesian indexes don't have to match shape, apparently!
@inline function getindex(tvec::TransposedVector, i::CartesianIndex)
    @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(tvec.vec)[1] && check_tail_indices(i.I...))
        throw(BoundsError(tvec, i.I))
    end
    @inbounds return transpose(tvec.vec[i.I[2]])
end
@inline function setindex!(tvec::TransposedVector, v, i::CartesianIndex)
    @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(tvec.vec)[1] && check_tail_indices(i.I...))
        throw(BoundsError(tvec, i.I))
    end
    @inbounds tvec.vec[i.I[2]] = transpose(v)
end

@propagate_inbounds getindex(tvec::TransposedVector, ::CartesianIndex{0}) = getindex(tvec)
@propagate_inbounds getindex(tvec::TransposedVector, i::CartesianIndex{1}) = getindex(tvec, i.I[1])

@propagate_inbounds setindex!(tvec::TransposedVector, v, ::CartesianIndex{0}) = setindex!(tvec, v)
@propagate_inbounds setindex!(tvec::TransposedVector, v, i::CartesianIndex{1}) = setindex!(tvec, v, i.I[1])

@inline check_tail_indices(i1, i2) = true
@inline check_tail_indices(i1, i2, i3, is...) = i3 == 1 ? check_tail_indices(i1, i2, is...) :  false

# Some conversions
convert(::Type{AbstractVector}, tvec::TransposedVector) = tvec.vec
convert{V<:AbstractVector}(::Type{V}, tvec::TransposedVector) = convert(V, tvec.vec)
convert{T,V}(::Type{TransposedVector{T,V}}, tvec::TransposedVector) = TransposedVector(convert(V, tvec.vec))

# helper function for below
@inline to_vec(tvec::TransposedVector) = transpose(tvec)
@inline to_vec(x::Number) = x
@inline to_vecs(tvecs...) = (map(to_vec, tvecs)...)

# map
@inline map(f, tvecs::TransposedVector...) = TransposedVector(map(f, to_vecs(tvecs...)...))

# broacast
@inline broadcast(f, tvecs::Union{Number,TransposedVector}...) = TransposedVector(broadcast(f, to_vecs(tvecs...)...))
# Other combinations default to Array...


text:
using SeisIO, JLD2, Wavelets, Statistics, DSP, FFTW, Plots, GLM, DataFrames

include("wavelet_methods.jl")

function test_WXS(finame::String, foname::String, InputDict::Dict, type::String; plot::Bool=true)
    valData = jldopen(finame)
    dvVlist = valData["info/dvVlist"]
    noiselist = valData["info/noiselist"]

    fmin    = InputDict["fmin"]
    fmax    = InputDict["fmax"]
    tmin    = InputDict["tmin"]
    tmax    = InputDict["tmax"]
    dt      = InputDict["dt"]
    dj      = InputDict["dj"]
    s0      = InputDict["s0"]
    J       = InputDict["J"]
    wvn     = InputDict["wvn"]
    sig     = InputDict["sig"]
    fbands  = InputDict["fbands"]
    unwrap  = InputDict["unwrap"]
    mintime = InputDict["mintime"]
    winflag = InputDict["winflag"]
    nwindow = InputDict["nwindow"]
    siglvl  = InputDict["siglvl"]
    twindow = [tmin, tmax]
    fwindow = [fmin, fmax]

    # output file
    f_out = jldopen(foname, "w")

    for dvV in dvVlist
        println("dv/v: $dvV")
        for noiselvl in noiselist
            signals = valData["$type/$dvV.$noiselvl"]
            ref = signals["ref"]
            cur = signals["cur"]

            npts = length(cur)
            tvec = (collect(0:npts-1) .* dt) .+ mintime

            # Wavelet cross-spectrum on discrete frequency bands
            (freqbands, dvv, err) = wxs_freqbands(cur, ref, tvec, twindow, fbands, dj, s0, J, wvn=wvn, unwrapflag=unwrap, sig=sig, siglvl=siglvl, windowflag=winflag, nwindow=nwindow)
            f_out["bands/$dvV/$noiselvl"] = [freqbands, dvv, err]
            p=PlotlyJS.plot(dvv, mode="lines+markers")
            display(p)
            readline()

            # Wavelet cross-spectrum for each frequency
            (freqs, dvv, err) = wxs_allfreq(cur, ref, tvec, twindow, fwindow, dj, s0, J, wvn=wvn, unwrapflag=unwrap, sig=sig, siglvl=siglvl, windowflag=winflag, nwindow=nwindow)
            f_out["full/$dvV/$noiselvl"] = [freqbands, dvv, err]
            p=PlotlyJS.plot(freqs, dvv, mode="lines+markers")
            display(p)
            readline()
        end
    end
    close(f_out)
end

function test_WTS(finame::String, foname::String, InputDict::Dict, type::String, plot::Bool=true)
    valData = jldopen(finame)
    dvVlist = valData["info/dvVlist"]
    noiselist = valData["info/noiselist"]

    fmin    = InputDict["fmin"]
    fmax    = InputDict["fmax"]
    tmin    = InputDict["tmin"]
    tmax    = InputDict["tmax"]
    dt      = InputDict["dt"]
    dj      = InputDict["dj"]
    s0      = InputDict["s0"]
    J       = InputDict["J"]
    wvn     = InputDict["wvn"]
    ndv     = InputDict["ndv"]
    norm    = InputDict["norm"]
    dvmax   = InputDict["dvmax"]
    fbands  = InputDict["fbands"]
    mintime = InputDict["mintime"]
    winflag = InputDict["winflag"]
    nwindow = InputDict["nwindow"]

    twindow = [tmin, tmax]
    fwindow = [fmin, fmax]

    # output file
    f_out = jldopen(foname, "w")

    for dvV in dvVlist
        println("dv/v: $dvV")
        for noiselvl in noiselist
            signals = valData["$type/$dvV.$noiselvl"]
            ref = signals["ref"]
            cur = signals["cur"]

            npts = length(cur)
            tvec = (collect(0:npts-1) .* dt) .+ mintime

            # Wavelet trace stretching on discrete frequency bands
            (freqbands, dvv, cc, cdp, err) = wts_freqbands(cur, ref, tvec, twindow, fbands, dj, s0, J, dvmax=dvmax, ndv=ndv, wvn=wvn, normalize=norm, windowflag=winflag, nwindow=nwindow)
            f_out["bands/$dvV/$noiselvl"] = [freqbands, dvv, err, cc, cdp]

            # Wavelet trace stretching for each frequency
            (freqs, dvv, cc, cdp, err) = wts_allfreqs(cur, ref, tvec, twindow, fbands, dj, s0, J, dvmax=dvmax, ndv=ndv, wvn=wvn, normalize=norm, windowflag=winflag, nwindow=nwindow)
            f_out["full/$dvV/$noiselvl"] = [freqs, dvv, err, cc, cdp]
        end
    end
    close(f_out)
end

function test_WTDTW(finame::String, foname::String, InputDict::Dict, type::String; plot::Bool=true)
    valData = jldopen(finame)
    dvVlist = valData["info/dvVlist"]
    noiselist = valData["info/noiselist"]

    fmin    = InputDict["fmin"]
    fmax    = InputDict["fmax"]
    tmin    = InputDict["tmin"]
    tmax    = InputDict["tmax"]
    dt      = InputDict["dt"]
    dj      = InputDict["dj"]
    s0      = InputDict["s0"]
    J       = InputDict["J"]
    b       = InputDict["b"]
    wvn     = InputDict["wvn"]
    ndv     = InputDict["ndv"]
    dir     = InputDict["dir"]
    norm    = InputDict["norm"]
    dvmax   = InputDict["dvmax"]
    fbands  = InputDict["fbands"]
    maxLag  = InputDict["maxLag"]
    winflag = InputDict["winflag"]
    mintime = InputDict["mintime"]
    nwindow = InputDict["nwindow"]

    twindow = [tmin, tmax]
    fwindow = [fmin, fmax]

    # output file
    f_out = jldopen(foname, "w")

    for dvV in dvVlist
        println("dv/v: $dvV")
        for noiselvl in noiselist
            signals = valData["$type/$dvV.$noiselvl"]
            ref = signals["ref"]
            cur = signals["cur"]

            npts = length(cur)
            tvec = (collect(0:npts-1) .* dt) .+ mintime

            # Wavelet dynamic time warping on discrete frequency bands
            (freqbands, dvv, err) = wtdtw_freqbands(cur, ref, tvec, twindow, fbands, maxLag, b, dir, dj, s0, J, wvn=wvn, normalize=norm, windowflag=winflag, nwindow=nwindow)
            f_out["bands/$dvV/$noiselvl"] = [freqbands, dvv, err]

            # Wavelet dynamic time warping for each frequency
            (freqs, dvv, err) = wtdtw_allfreqs(cur, ref, tvec, twindow, fwindow, maxLag, b, dir, dj, s0, J, wvn=wvn, normalize=norm, windowflag=winflag, nwindow=nwindow)
            f_out["full/$dvV/$noiselvl"] = [freqs, dvv, err]
        end
    end
    close(f_out)
end

finame = "/Users/jared/SCECintern2019/data/validation/verificationData.jld2"
foname_WXS = "verificationData_WXS.jld2"
foname_WTS = "verificationData_WTS.jld2"
foname_DTW = "verificationData_DTW.jld2"

InputDict_WXS = Dict( "fmin"    => 0.0,
                      "fmax"    => 10.0,
                      "tmin"    => 0.0,
                      "tmax"    => 200.0,
                      "dt"      => 0.05,
                      "dj"      => 1/12,
                      "s0"      => -1,
                      "J"       => -1,
                      "wvn"     => "morlet",
                      "sig"     => false,
                      "unwrap"  => false,
                      "mintime" => 0,
                      "siglvl"  => 0.95,
                      "fbands"  => [0.1 0.2; 0.2 0.4; 0.4 0.8; 0.8 1.6; 1.6 3.2],
                      "winflag" => true,
                      "nwindow" => 1.5)

InputDict_WTS = Dict( "fmin"    => 0.0,
                      "fmax"    => 10.0,
                      "tmin"    => 0.0,
                      "tmax"    => 200.0,
                      "dt"      => 0.05,
                      "dj"      => 1/12,
                      "s0"      => -1,
                      "J"       => -1,
                      "wvn"     => "morlet",
                      "ndv"     => 100,
                      "dvmax"   => 0.1,
                      "mintime" => 0,
                      "norm"    => true,
                      "fbands"  => [0.1 0.2; 0.2 0.4; 0.4 0.8; 0.8 1.6; 1.6 3.2],
                      "winflag" => true,
                      "nwindow" => 1.5)

InputDict_DTW = Dict( "fmin"    => 0.0,
                      "fmax"    => 10.0,
                      "tmin"    => 0.0,
                      "tmax"    => 200.0,
                      "dt"      => 0.05,
                      "dj"      => 1/12,
                      "s0"      => -1,
                      "b"       => 1,
                      "J"       => -1,
                      "wvn"     => "morlet",
                      "dir"     => 0,
                      "ndv"     => 100,
                      "dvmax"   => 0.1,
                      "mintime" => 0,
                      "norm"    => true,
                      "maxLag"  => 150,
                      "fbands"  => [0.1 0.2; 0.2 0.4; 0.4 0.8; 0.8 1.6; 1.6 3.2],
                      "winflag" => true,
                      "nwindow" => 1.5)

for type in ["spectSynth", "rickerConv", "dampedSin", "realData"]
    test_WXS(finame, foname_WXS, InputDict_WXS, type)
    test_WTS(finame, foname_WTS, InputDict_WTS, type)
    test_WTDTW(finame, foname_DTW, InputDict_DTW, type)
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
###
# This script generates julia bindings for different Vulkan versions

using Clang.wrap_c
using Clang.cindex

const LLVM_VERSION = readchomp(`llvm-config --version`)
const LLVM_LIBDIR  = readchomp(`llvm-config --libdir`)
const LLVM_INCLUDE = joinpath(LLVM_LIBDIR, "clang", LLVM_VERSION, "include")

# Which variable contains the library later onr
header_library(x) = "libvulkan"

git = `git -C Vulkan-Docs`
# Setup Vulkan-Docs for getting the headers
if !isdir("Vulkan-Docs")
  run(`git clone https://github.com/KhronosGroup/Vulkan-Docs`)
else
  run(`$git fetch`)
end

const VK_INCLUDE = "Vulkan-Docs/src/vulkan"
const VK_HEADERS = map(x->joinpath(VK_INCLUDE, x), ["vk_platform.h", "vulkan.h"])
const VK_EXTENSIONS = [
  "VK_USE_PLATFORM_XCB_KHR",
  "VK_USE_PLATFORM_WAYLAND_KHR",
  "VK_USE_PLATFORM_MIR_KHR",
  "VK_USE_PLATFORM_WIN32_KHR",
  ]

# Map versions to revisions
const VK_VERSIONS = Dict(
  v"1.0.3" => "d204ac2",
  v"1.0.4" => "5a4c5e5",
  v"1.0.5" => "7380aee",
  v"1.0.6" => "13f85fa", # from here one VK_HEADER_VERSION is introduced
  v"1.0.7" => "8c3c9b4",
  v"1.0"   => "1.0", # Generate from master branch
)

# callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::ASCIIString, cursor_header::ASCIIString)
  return startswith(dirname(cursor_header), VK_INCLUDE)
end

# These rewriters are taken from CUDArt.jl
function rewriter(ex :: Expr)
  # Empty types get converted to Void
  # This is important for opaque handles
  if ex.head == :type
    a3 = ex.args[3]
    if isempty(a3.args)
      objname = ex.args[2]
      return :(typealias $objname Void)
    end
  end

  # Early exit for everything but functions
  ex.head == :function || return ex

  decl, body = ex.args[1], ex.args[2]
  # omit types from function prototypes
  for i = 2:length(decl.args)
    a = decl.args[i]
    # a can be a symbol (and thus already have no type information attached)
    if !(typeof(a) == Symbol) && a.head == :(::)
      decl.args[i] = a.args[1]
    end
  end

  return ex
end

rewriter(A::Array) = [rewriter(a) for a in A]

rewriter(arg) = arg

function generate_bindings(version = v"1.0")
  clang_includes = ASCIIString[]
  push!(clang_includes, LLVM_INCLUDE)
  push!(clang_includes, VK_INCLUDE)

  clang_extraargs = ASCIIString[]
  for extension in VK_EXTENSIONS
    push!(clang_extraargs, "-D")
    push!(clang_extraargs, extension)
  end

  git_commit = VK_VERSIONS[version]

  run(`$git checkout $git_commit`)

  const wc = wrap_c.init(;
                        headers = VK_HEADERS,
                        output_file = "api/vk_$(version).jl",
                        common_file = "api/vk_common_$(version).jl",
                        clang_includes = clang_includes,
                        clang_args = clang_extraargs,
                        header_wrapped = wrap_header,
                        header_library = header_library,
                        clang_diagnostics = false,
                        rewriter = rewriter)

  # wrap_structs, immutable_structs
  wc.options = wrap_c.InternalOptions(true, true)
  run(wc)
end

function generate_spirv()
  download("https://www.khronos.org/registry/spir-v/api/1.0/spirv.h", "spirv.h")

  clang_includes = ASCIIString[]
  push!(clang_includes, LLVM_INCLUDE)

  const wc = wrap_c.init(;
                        headers = ["./spirv.h"],
                        output_file = "api/spirv.jl",
                        common_file = "api/spirv_common_.jl",
                        clang_includes = clang_includes,
                        header_library = x->"libspirv",
                        clang_diagnostics = false,
                        rewriter = rewriter)

  # wrap_structs, immutable_structs
  wc.options = wrap_c.InternalOptions(true, true)
  run(wc)
end

generate_bindings()
generate_spirv()


text:
const DEBUG_ENABLED = Ref(false)
const DEBUG_CALLBACK = Ref{Function}()

@export struct GLDebugInfo <: Iterable
   type::String
   source::String
   message::String
   severity::String
end

function debug_message_callback(
   _source::GLenum,
   _type::GLenum,
   ::GLuint,
   _severity::GLenum,
   ::GLsizei,
   _msg::Ptr{GLchar},
   ::Ptr{Cvoid}
)::Cvoid
   message  = unsafe_string(_msg)
   source   = get(DebugSource, _source, "Unknown")
   type     = get(DebugType, _type, "Unknown")
   severity = get(DebugSeverity, _severity, "Unknown")

   if isassigned(DEBUG_CALLBACK)
      GLDebugInfo(
         type,
         source,
         message,
         severity,
      ) |> DEBUG_CALLBACK[]
   end

   return nothing
end

"""
      ondebug_msg(callback::Function)

Enables debugging of OpenGL functions and sets the OpenGL debug message `callback` function.

`callback` only accepts one parameter of type [`GLDebugInfo`](@ref),
that contains the information about the debug message.
"""
@export function ondebug_msg(callback::Function)
   @assert hasmethod(callback, Tuple{GLDebugInfo}) "
      The `ondebug_msg` function expects a `callback` that accepts
      a single parameter of type `GLDebugInfo` : `callback(info::GLDebugInfo)`
   "

   DEBUG_CALLBACK[] = callback

   if DEBUG_ENABLED[] == false
      DEBUG_ENABLED[] = true

      ptr_debug_callback = @cfunction(debug_message_callback, Cvoid, (
         GLenum,
         GLenum,
         GLuint,
         GLenum,
         GLsizei,
         Ptr{GLchar},
         Ptr{Cvoid}
      ))

      glEnable(GL_DEBUG_OUTPUT)
      glEnable(GL_DEBUG_OUTPUT_SYNCHRONOUS)

      GC.@preserve ptr_debug_callback begin
         glDebugMessageCallback(ptr_debug_callback, C_NULL)
         ids = Cuint(0)
         @c glDebugMessageControl(
            GL_DONT_CARE,
            GL_DONT_CARE,
            GL_DONT_CARE,
            0,
            &ids,
            GL_TRUE,
         )
      end
   end
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
# stuff in this file takes Julia types and converts them to Ctypes compatible
# to be used inside of

# Convert Data structures to C
function to_mtl_ref(arg::T) where T
    if !Base.datatype_pointerfree(T)
        error("Types should not contain pointers: $T")
    end
    if contains_different_layout(T)
        x = replace_different_layout(arg)
        return Base.RefValue(x), sizeof(x)
    end
    Base.RefValue(arg), sizeof(arg)
end

## Identify if there is a mtlbuffer in this structure...

# returns true if structure contains a device buffer, false otherwise
function _contains_mtlbuffer(::Type{T}) where T
    if T.hasfreetypevars
        throw("Only concrete types can be sent to GPU!")
    end
    for typ in T.types
        typ <: DeviceBuffer && return true
    end

    return false
end

@generated function contains_mtlbuffer(::Type{T}) where T
    :($(_contains_mtlbuffer(T)))
end

##
function _contains_different_layout(::Type{T}) where T
    sizeof(T) == 0 && return true
    nfields(T) == 0 && return false
    for fname in fieldnames(T)
        contains_different_layout(fieldtype(T, fname)) && return true
    end
    return false
end

contains_different_layout(::Type{NTuple{3, T}}) where {T <: Union{Float32, Float64, Int8, Int32,
                                                                  Int64, UInt8, UInt32, UInt64}} = true

"""
    contains_different_layout(T)
Empty types and NTuple{3, CLNumber} have different layouts and need to be replaced
(Where `CLNumber <: Union{Float32, Float64, Int8, Int32, Int64, UInt8, UInt32, UInt64}`)
TODO: Float16 + Int16 should also be in CLNumbers
"""
@generated function contains_different_layout(::Type{T}) where T
    :($(_contains_different_layout(T)))
end

function struct2tuple(x::T) where T
    ntuple(nfields(x)) do i
        getfield(x, i)
    end
end

replace_different_layout(red::NTuple{N, Any}, rest::Tuple{}) where N = red
function replace_different_layout(red::NTuple{N, Any}, rest) where N
    elem1 = first(rest)
    T = typeof(elem1)
    repl = if sizeof(T) == 0 && nfields(elem1) == 0
        Int32(0)
    elseif contains_different_layout(T)
        replace_different_layout(elem1)
    else
        elem1
    end
    replace_different_layout((red..., repl), Base.tail(rest))
end

# TODO UInt16/Float16?
# Handle different sizes of OpenCL Vec3, which doesn't agree with julia
function replace_different_layout(arg::NTuple{3, T}) where T <: Union{Float32, Float64, Int8, Int32, Int64, UInt8, UInt32, UInt64}
    pad = T(0)
    (arg..., pad)
end

Base.@pure datatype_align(x::T) where {T} = datatype_align(T)
Base.@pure function datatype_align(::Type{T}) where {T}
    # typedef struct {
    #     uint32_t nfields;
    #     uint32_t alignment : 9;
    #     uint32_t haspadding : 1;
    #     uint32_t npointers : 20;
    #     uint32_t fielddesc_type : 2;
    # } jl_datatype_layout_t;
    field = T.layout + sizeof(UInt32)
    unsafe_load(convert(Ptr{UInt16}, field)) & convert(Int16, 2^9-1)
end


text:
mutable struct FileCachedBoolVariate{S<:AbstractBoolVariate} <: AbstractBoolVariate
    cachepath::String
    basevariate::S
    lockobj::Threads.TatasLock

    FileCachedBoolVariate{S}(cachepath, basevariate, lockobj) where {S<:AbstractBoolVariate} = (res = new(cachepath, basevariate, lockobj); finalizer(res, x -> rm(x.cachepath; force = true)); res)
end

Base.length(var::FileCachedBoolVariate) = length(var.basevariate)

getname(var::FileCachedBoolVariate) = getname(var.basevariate)

function FileCachedBoolVariate(basevariate::S) where {S<:AbstractBoolVariate}
    FileCachedBoolVariate{S}("", basevariate, Threads.TatasLock())
end

function filecache(basevariate::S) where {S<:AbstractBoolVariate}
    FileCachedBoolVariate{S}("", basevariate, Threads.TatasLock())
end

function slice(boolvar::FileCachedBoolVariate{S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:Union{FileCachedBoolVariate, CachedBoolVariate, BoolVariate}}
    basevariate = boolvar.basevariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(basevariate, fromobs, toobs, slicelength)
end

function slice(boolvar::FileCachedBoolVariate{S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractBoolVariate}
    basevariate = boolvar.basevariate
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    lockobj = boolvar.lockobj
    lock(lockobj)
    try
        if !isfile(boolvar.cachepath)
            v = convert(Vector{Bool}, basevariate)
            tmpjuml = joinpath(tempdir(), JUMLDIR)
            mkpath(tmpjuml)
            outdatpath = joinpath(tmpjuml, "$(randstring(10)).dat")
            open(outdatpath, "w") do f
                write(f, v)
            end
            boolvar.cachepath = outdatpath
        end
    finally
        unlock(lockobj)
    end
    sizeT = sizeof(Bool)
    iostream = open(boolvar.cachepath, "r")
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{Bool}(undef, slicelength)
    map(Seq(Vector{Bool}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{Bool,1,Array{Bool,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(boolvariate::FileCachedBoolVariate, dataframe::AbstractDataFrame)
    FileCachedBoolVariate(map(boolvariate.basevariate, dataframe))
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# Make the parameter structure of a dense layer
dense_layer(n_in, n_out) = ComponentArray(W=Flux.glorot_uniform(n_out, n_in), b=zeros(Float32, n_out))

# Returns a function that runs the simulation for a given control parameters
function make_run_func(prob; kwargs...)
    return function (control)
        p = deepcopy(prob.p)
        p = (; p..., control)
        u0 = eltype(control).(prob.u0)
        new_prob = remake(prob; p=p, u0=u0)
        return solve(new_prob, Tsit5(); kwargs...)
    end
end

function make_objective(prob, out_fun; kwargs...)
    run_func = make_run_func(prob; kwargs...)
    return function (x)
        sol = run_func(x)
        return out_fun(sol), sol
    end
end

function make_callback(; plot_every=1, doplot=false)
    # plt = plot()
    all_loss = Float64[]
    i = 0
    return function (θ, loss, pred; doplot=doplot)
        display(loss)
        if doplot
            if i % plot_every == 0
                push!(all_loss, loss)
                display(plot((0:length(all_loss)-1)*plot_every, all_loss))
            end
        end
        i += 1
        return false
    end
end

# # Don't think we'll actually use this
# function run_neural(x, layers, act_fun, out_fun=identity)
#     out =  reduce(valkeys(layers); init=x) do y, key
#         @unpack W, b = layers[key]
#         return act_fun.(W*y .+ b)
#     end
#     return out_fun(out)
# end


text:
n = 17
ijs = [(3,1), (3,2), (4,1), (4,2), (4,3), (5,1), (5,3), (5,4), (8,7),
       (9,5), (9,6), (9,7), (9,8), (11,10), (13,10), (13,11), (13,12),
       (14,10), (14,11), (14,12), (14,13),
       (15, 1), (15,3), (15,4), (15,5), (15,7), (15,8), (15,9),
       (16,5), (16,6), (16,9), (16,12), (16,13), (16,14), (16,15),
       (17,10), (17,11), (17,12), (17,13), (17,14), (17,15), (17,16)]
append!(ijs, [(i,i) for i in 1:n])
I, J = Chordal.unzip(ijs)
A = sparse(I, J, ones(length(I)))
A = A + A'
L = qdldl(A; perm=nothing, logical=true).L

@testset "Test elimination tree" begin
    # VA survey, Figure 4.2
    par_true = vec([3 3 4 5 9 9 8 9 15 11 13 13 14 16 16 17 0])
    ch_true =  [Set{Int}() for _ in 1:n]
    ch_true[3] = Set([1, 2])
    ch_true[4] = Set([3])
    ch_true[5] = Set([4])
    ch_true[8] = Set([7])
    ch_true[9] = Set([5, 6, 8])
    ch_true[11] = Set([10])
    ch_true[13] = Set([11, 12])
    ch_true[14] = Set([13])
    ch_true[15] = Set([9])
    ch_true[16] = Set([15, 14])
    ch_true[17] = Set([16])

    et_par = CD.get_etree(L)
    @test all(et_par .== par_true)
    et_ch = CD.get_children_from_par(et_par)
    @test all([issetequal(et_ch[i], ch_true[i]) for i in 1:n])
    perm = CD.get_postordering(et_par, et_ch)
    iperm = invperm(perm)
    @test all([iperm[i] .< et_par[iperm[i]] for i in 1:17 if et_par[iperm[i]] ≠ 0])
end

@testset "Test Clique Tree" begin
    deg⁺_true = vec([4 2 3 2 3 2 3 2 2 4 3 4 3 2 2 1 0])

    et_par = CD.get_etree(L)
    et_ch = CD.get_children_from_par(et_par)
    deg⁺ = CD.get_higher_deg(L)
    @test all(deg⁺_true .== deg⁺)

    vreps, snd_par, snd_mem = CD.max_supernode_etree(L, et_par)
    # Test vrep (Thm 4.2)
    for v in vreps
        @test all([deg⁺[w] < deg⁺[v] + 1 for w in et_ch[v]])
    end
    # Test supernodes
    for (ind, v) in enumerate(vreps)
        nv = count(x->x==ind, snd_mem) - 1
        pv = v
        for k = 1:nv
            pv = et_par[pv]
            @test deg⁺[v] == deg⁺[pv] + k
            @test snd_mem[pv] == ind
        end
    end
    # Test clique tree (Thm 4.3)
    # Test root
    root_et = findfirst(x->x==0, et_par)
    root_ind = findfirst(x->x==0, snd_par)
    @test snd_mem[root_et] == root_ind
    for (ind, v) in enumerate(vreps)
        ind == root_ind && continue
        col_v = vcat([v], rowvals(L)[nzrange(L, v)])
        qv = vreps[snd_par[ind]]
        col_qv = vcat([qv], rowvals(L)[nzrange(L, qv)])
        col_minus_snd = filter(x->(snd_mem[x] != ind), col_v)
        # println(" --- v = $v --- ")
        # @show col_minus_snd
        # @show col_qv
        @test issubset(col_minus_snd, col_qv) && length(col_minus_snd) < length(col_qv)
    end

end


text:
using SmartShockFinder

datatyp = "NNRH"
n_dims = 2
η = 0.005                # learning rate
β = 0.002                # regularization paramater
number_epochs = 100
Sb = 500                 # batch size
L = 10                   # early stopping parameter
hidden_units_1 = 40
hidden_units_2 = 20
hidden_units_3 = 20
hidden_units_4 = 20
hidden_units_5 = 10

train_network2d(datatyp, n_dims, η, β, number_epochs, Sb, L,
                hidden_units_1, hidden_units_2, hidden_units_3,
                hidden_units_4, hidden_units_5)


text:
using Distributed
addprocs(2)
nprocs()
@everywhere using FRGRealTime
@everywhere using FRGRealTime.ThresholdFunctions
@everywhere using DifferentialEquations
@everywhere using HCubature
using Plots
@everywhere include("flow.jl")
include("solver.jl")

@everywhere function tchanelLPASolve(
    T,
    IRScale,
    UVScale,
    massflowfun,
    lambdaflowfun;
    u0,
    progress = true,
    atol = 1e-14,
    rtol = 1e-14,
    adaptive = true,
    dtmax = 0.01,
    kwargs...,
)
    function flow_LPA(du, u, p, k)
        #display(plot(x->-h2(x)[2],k,Λ))
        #println(k, " m2=", -u[2], " λ=", u[1])
        du[1] = massflowfun(k, u[1], p, u[2])
        du[2] = lambdaflowfun(k, u[1], p, u[2])
    end
    alg = Tsit5()
    kspan = (UVScale, IRScale)
    prob = ODEProblem(flow_LPA, u0, kspan, T)
    solve(
        prob,
        alg,
        progress = progress,
        atol = atol,
        rtol = rtol,
        adaptive = adaptive,
        dtmax = dtmax,
        kwargs...,
    )
end

@everywhere sol1 = tchanelLPASolve(
    150.0,
    1.0,
    800.0,
    massflow2,
    lambdaflow2;
    u0 = [-116585.3, 6 * 8.0],
)

@everywhere sol1 = tchanelLPASolve(
    145.0,
    1.0,
    1000.0,
    massflow1,
    lambdaflow1;
    u0 = [-146585.3, 5 * 7.0],
)

plot(k -> sol1(k)[1], 1.0, 800.0)

@time tsol1 = tchanelSolve(
    2000.0,
    145.0,
    1.0,
    1000.0,
    x -> sol1(x)[1],
    Imlambdaflow1,
    Relambdaflow1,
    0.1,
    u0 = [-0.0001, 5 * 7.0],
    dtmax = 0.5,
)

@time tsolparallel1 = tchanelSolveParallel(
    145.0,
    1.0,
    1000.0,
    x -> sol1(x)[1],
    Imlambdaflow1,
    Relambdaflow1,
    4000,
    0.2,
    u0 = [-0.0001, 5 * 7.0],
    dtmax = 0.05,
)

@time tsolTwoPointparallel1 = tchanelSolveTwoPointParallel(
    150.0,
    1.0,
    800.0,
    x -> sol1(x)[1],
    x -> sol1(x)[2],
    Imlambdaflow2,
    Relambdaflow2,
    100,
    0.05,
    u0 = [-0.00001, 6 * 8.0],
    dtmax = 0.04,
)

imvec=hcat((tsolTwoPointparallel1.u)...)'[:,1]
revec=hcat((tsolTwoPointparallel1.u)...)'[:,2]

lambda_omega = range(0.0, stop = 400.0, length = 100)

plot(lambda_omega,revec)

plot!(lambda_omega,abs.(revec),yaxis=:log)

plot(lambda_omega,-imvec./(imvec.^2+revec.^2),yaxis=:log)
tsolparallel1[3](40.0)

plot(k -> sol1(k)[1], 1.0, 800.0)

plot(k -> -tsol1(k)[1], 1, 1000.0)

plot(k -> tsol1(k)[2], 1.0, 800.0, yaxis = :log)


text:
"""
    push!!(collection, items...)

Push one or more `items` to collection.  Create a copy of `collection`
if it cannot be mutated or the element type does not match.

# Examples
```jldoctest
julia> using BangBang

julia> push!!((1, 2), 3)
(1, 2, 3)

julia> push!!([1, 2], 3)
3-element Array{Int64,1}:
 1
 2
 3

julia> push!!([1, 2], 3.0)
3-element Array{Float64,1}:
 1.0
 2.0
 3.0
```
"""
push!!(xs, i1, i2, items...) =
    foldl(push!!, items, init=push!!(push!!(xs, i1), i2))
push!!(xs, x) = may(push!, xs, x)

pure(::typeof(push!)) = NoBang.push
_asbb(::typeof(push!)) = push!!
possible(::typeof(push!), x, ::Any) = ismutable(x)
possible(::typeof(push!), ::C, ::S) where {C <: MaybeMutableContainer, S} =
    ismutable(C) && promote_type(eltype(C), S) <: eltype(C)

"""
    append!!(dest, src)

Append items in `src` to `dest`.  Mutate `dest` if possible.  See also
[`push!!`](@ref).

# Examples
```jldoctest
julia> using BangBang

julia> append!!((1, 2), (3, 4))
(1, 2, 3, 4)

julia> append!!([1, 2], (3, 4))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
"""
append!!(xs, ys) = may(append!, xs, ys)

pure(::typeof(append!)) = NoBang.append
_asbb(::typeof(append!)) = append!!
possible(::typeof(append!), x, ::Any) = ismutable(x)
possible(::typeof(append!), ::C, ys) where {C <: MaybeMutableContainer} =
    ismutable(C) && promote_type(eltype(C), eltype(ys)) <: eltype(C)

"""
    pushfirst!!(collection, items...)

# Examples
```jldoctest
julia> using BangBang

julia> pushfirst!!((1, 2), 3, 4)
(3, 4, 1, 2)

julia> pushfirst!!([1, 2], 3, 4)
4-element Array{Int64,1}:
 3
 4
 1
 2

julia> pushfirst!!([1, 2], 3, 4.0)
4-element Array{Float64,1}:
 3.0
 4.0
 1.0
 2.0
```
"""
pushfirst!!(xs, ys...) = may(pushfirst!, xs, ys...)

pure(::typeof(pushfirst!)) = NoBang.pushfirst
_asbb(::typeof(pushfirst!)) = pushfirst!!
possible(::typeof(pushfirst!), ::Tuple, ::Vararg) = false
possible(::typeof(pushfirst!), ::C, ys...) where {C <: AbstractVector} =
    ismutable(C) && promote_type(eltype(C), map(typeof, ys)...) <: eltype(C)

"""
    pop!!(sequence) -> (sequence′, value)
    pop!!(assoc, key) -> (assoc′, value)
    pop!!(assoc, key, default) -> (assoc′, value)

# Examples
```jldoctest
julia> using BangBang

julia> pop!!([0, 1])
([0], 1)

julia> pop!!((0, 1))
((0,), 1)

julia> pop!!(Dict(:a => 1), :a)
(Dict{Symbol,Int64}(), 1)

julia> pop!!((a=1,), :a)
(NamedTuple(), 1)
```
"""
pop!!(xs, args...) = may(_pop!, xs, args...)
_pop!(xs, args...) = xs, pop!(xs, args...)

pure(::typeof(_pop!)) = NoBang.pop
possible(::typeof(_pop!), ::C, ::Vararg) where C = ismutable(C)

"""
    delete!!(assoc, key) -> assoc′
"""
delete!!(xs, key) = may(delete!, xs, key)

pure(::typeof(delete!)) = NoBang.delete
_asbb(::typeof(delete!)) = delete!!
possible(::typeof(delete!), ::C, ::Any) where C = ismutable(C)

"""
    popfirst!!(sequence) -> (sequence′, value)

# Examples
```jldoctest
julia> using BangBang

julia> popfirst!!([0, 1])
([1], 0)

julia> popfirst!!((0, 1))
((1,), 0)
```
"""
popfirst!!(xs) = may(_popfirst!, xs)
_popfirst!(xs) = xs, popfirst!(xs)

pure(::typeof(_popfirst!)) = NoBang.popfirst
possible(::typeof(_popfirst!), ::C,) where C = ismutable(C)

"""
    empty!!(collection) -> collection′

# Examples
```jldoctest
julia> using BangBang

julia> empty!!((1, 2, 3))
()

julia> empty!!((a=1, b=2, c=3))
NamedTuple()

julia> xs = [1, 2, 3];

julia> empty!!(xs)
0-element Array{Int64,1}

julia> xs
0-element Array{Int64,1}
```
"""
empty!!(xs) = may(empty!, xs)

pure(::typeof(empty!)) = NoBang._empty
_asbb(::typeof(empty!)) = empty!!
possible(::typeof(empty!), ::C) where C = ismutable(C)

#=
"""
    splice!!(sequence, i, [replacement]) -> (sequence′, item)
    splice!!(sequence, range, [replacement]) -> (sequence′, items)
"""
=#

"""
    splice!!(sequence, i, [replacement]) -> (sequence′, item)
"""
splice!!(xs, args...) = may(_splice!, xs, args...)
_splice!(xs, args...) = xs, splice!(xs, args...)

pure(::typeof(_splice!)) = NoBang.splice
possible(::typeof(_splice!), ::C, ::Any) where C = ismutable(C)
possible(::typeof(_splice!), xs::C, i::Integer, ::S) where {C, S} =
    possible(_splice!, xs, i) && promote_type(eltype(C), S) <: eltype(C)
possible(::typeof(_splice!), xs::C, i::Any, ::S) where {C, S} =
    possible(_splice!, xs, i) && promote_type(eltype(C), eltype(S)) <: eltype(C)

"""
    setindex!!(collection, value, indices...) -> collection′
"""
Base.@propagate_inbounds setindex!!(xs, v, I...) = may(_setindex!, xs, v, I...)

Base.@propagate_inbounds _setindex!(xs, v, I...) = (setindex!(xs, v, I...); xs)

pure(::typeof(_setindex!)) = NoBang._setindex
possible(::typeof(_setindex!), ::Union{Tuple, NamedTuple}, ::Vararg) = false
possible(::typeof(_setindex!), ::C, ::T, ::Vararg) where {C <: AbstractArray, T} =
    ismutable(C) && promote_type(eltype(C), T) <: eltype(C)
possible(::typeof(_setindex!), ::C, ::V, ::K) where {C <: AbstractDict, V, K} =
    ismutable(C) &&
    promote_type(keytype(C), K) <: keytype(C) &&
    promote_type(valtype(C), V) <: valtype(C)

"""
    setproperty!!(value, name, x)

# Examples
```jldoctest
julia> using BangBang

julia> setproperty!!((a=1, b=2), :b, 3)
(a = 1, b = 3)

julia> struct Immutable
           a
           b
       end

julia> setproperty!!(Immutable(1, 2), :b, 3)
Immutable(1, 3)

julia> mutable struct Mutable
           a
           b
       end

julia> s = Mutable(1, 2);

julia> setproperty!!(s, :b, 3)
Mutable(1, 3)

julia> s
Mutable(1, 3)
```
"""
setproperty!!(value, name, x) = may(_setproperty!, value, name, x)

function _setproperty!(value, name, x)
    setproperty!(value, name, x)
    return value
end

pure(::typeof(_setproperty!)) = NoBang.setproperty
possible(::typeof(_setproperty!), x, ::Any, ::Any) = ismutablestruct(x)

"""
    materialize!!(dest, x)
"""
@inline materialize!!(dest, x) = may(materialize!, dest, x)
# TODO: maybe instantiate `x` and be aware of `x`'s style

pure(::typeof(materialize!)) = NoBang.materialize
possible(::typeof(materialize!), x, ::Any) = ismutable(x)


text:
# ------------------------------------------------------------------------------------------
# # 软件包（Packages）
#
# Julia 有超过 2000 个登记在案的包，这使得第三方的软件包成为 Julia 语言生态的重要组成部分。
# 尽管如此，软件包生态依然有很大的成长空间。
#
# 值得一提的是，调用其它语言的函数在 Julia 中也是一等公民，，提供了很棒的外部函数接口（FFI）。
# 我们可以很容易地通过 `PyCall` 或 `Rcall` 调用 Python 或 R。
# 这意味着你不必等到 Julia 的生态完全成熟就能迁移现有项目到 Julia，
# 而且迁移到 Julia 并不意味着放弃你在其他语言中惯用的包或者库！
#
# 查看全部可获得的包，请看
# https://pkg.julialang.org/
# 或者
# https://juliaobserver.com/
#
# 现在我们先学习一下该如何使用包。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 在安装好的 Julia 中第一次要使用一个包的时候，你需要用包管理器（package manager）明确添加（add）它：
# ------------------------------------------------------------------------------------------

using Pkg
Pkg.add("Example")

# ------------------------------------------------------------------------------------------
# 之后每次使用 Julia 时（比如打开一个新的 REPL 会话或者第一次打开一个 notebook），
# 你需要通过 `using` 关键字载入要用的包
# ------------------------------------------------------------------------------------------

using Example

# ------------------------------------------------------------------------------------------
# 在 `Example.jl` 的源码
# https://github.com/JuliaLang/Example.jl/blob/master/src/Example.jl
# 中我们可以看到声明了下面这个函数
#
# ```julia
# hello(who::String) = "Hello, $who"
# ```
#
# 载入了 `Example` 我们现在应该可以调用 `hello` 函数了
# ------------------------------------------------------------------------------------------

hello("it's me. I was wondering if after all these years you'd like to meet.")

# ------------------------------------------------------------------------------------------
# 现在让我们来玩 `Colors` 包吧
# ------------------------------------------------------------------------------------------

Pkg.add("Colors")

using Colors

# ------------------------------------------------------------------------------------------
# 我们来创建一个 100 个不同颜色的调色板
# ------------------------------------------------------------------------------------------

palette = distinguishable_colors(100)

# ------------------------------------------------------------------------------------------
# 然后通过 `rand` 函数来创建随机矩阵
# ------------------------------------------------------------------------------------------

rand(palette, 3, 3)

# ------------------------------------------------------------------------------------------
# 在下一个 notebook 中，我们将使用一个新的包来绘制数据集。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 练习
#
# #### 7.1
# 加载 `Primes` 包（源码在https://github.com/JuliaMath/Primes.jl ）
# ------------------------------------------------------------------------------------------

# 在此作答。做点什么以安装要求的包。

using Primes

@assert @isdefined Primes

# ------------------------------------------------------------------------------------------
# #### 7.2
# 使用 `primes` 函数获得所有小于 1,000,000 的质数，并把它存入变量 `primes_list`
# ------------------------------------------------------------------------------------------

# 在此作答
primes_list = [];

@assert primes_list == primes(1000000)

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
import LinearAlgebra.BLAS.@blasfunc

for (fname, elty) in ((:dgbmv_,:Float64),
                      (:sgbmv_,:Float32),
                      (:zgbmv_,:ComplexF64),
                      (:cgbmv_,:ComplexF32))
    @eval begin
             # SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER INCX,INCY,KL,KU,LDA,M,N
             #       CHARACTER TRANS
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function gbmv!(trans::Char, m::Int, n::Int, kl::Int, ku::Int, alpha::($elty),
                       A::Ptr{$elty}, st::Int,
                       x::Ptr{$elty}, incx::Int, beta::($elty), y::Ptr{$elty}, incy::Int)
            ccall((@blasfunc($fname), libblas), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}),
                 trans, m, n, kl,
                 ku, alpha, A, st,
                 x, incx, beta, y, incy)
            y
        end
    end
end

for (fname, elty) in ((:dsbmv_,:Float64),
                      (:ssbmv_,:Float32),
                      (:zsbmv_,:ComplexF64),
                      (:csbmv_,:ComplexF32))
    @eval begin
                # SUBROUTINE DSBMV(UPLO, N, K, ALPHA, A, LDA,
                # X, INCX, BETA, Y, INCY)
                # DOUBLE PRECISION ALPHA,BETA
                # INTEGER INCX,INCY,K,LDA,N
                # CHARACTER*1 UPLO
                # DOUBLE PRECISION A(LDA,*), X(*), Y(*)

        function sbmv!(uplo::Char, n::Int, k::Int, alpha::($elty),
                       A::Ptr{$elty}, lda::Int,
                       x::Ptr{$elty}, incx::Int, beta::($elty), y::Ptr{$elty}, incy::Int)
            ccall((@blasfunc($fname), libblas), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}),
                 uplo, n, k,
                 alpha, A, lda,
                 x, incx, beta,
                 y, incy)
            y
        end
    end
end

sbmv!(uplo::Char, k::Int, alpha::T,
      A::AbstractMatrix{T}, x::AbstractVector{T}, beta::T, y::AbstractVector{T}) where {T<:BlasFloat} =
    sbmv!(uplo,size(A,2),k,alpha,pointer(A),max(1,stride(A,2)),x,stride(x,1),beta,y,stride(y,1))

## Triangular *

function tbmv! end

for (fname, elty) in ((:dtbmv_,:Float64),
                        (:stbmv_,:Float32),
                        (:ztbmv_,:ComplexF64),
                        (:ctbmv_,:ComplexF32))
    @eval begin
                #       SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
                # *     .. Scalar Arguments ..
                #       INTEGER INCX,LDA,N
                #       CHARACTER DIAG,TRANS,UPLO
                # *     .. Array Arguments ..
                #       DOUBLE PRECISION A(LDA,*),X(*)
        function tbmv!(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar,
                       m::Int, k::Int, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            require_one_based_indexing(A, x)
            n = size(A,2)
            size(A,1) ≥ k+1 || throw(ArgumentError("triangular banded data missing"))
            n == m || throw(DimensionMismatch("matrix is not square: dimensions are $n, $m"))
            if n != length(x)
                throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, diag, m, k,
                 A, max(1,stride(A,2)), x, max(1,stride(x, 1)))
            x
        end
        function tbmv(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, m::Int, k::Int,
                      A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            tbmv!(uplo, trans, diag, m, k, A, copy(x))
        end
    end
end

## Triangular \

function tbsv end

for (fname, elty) in ((:dtbsv_,:Float64),
                        (:stbsv_,:Float32),
                        (:ztbsv_,:ComplexF64),
                        (:ctbsv_,:ComplexF32))
    @eval begin
                #       SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
                #       .. Scalar Arguments ..
                #       INTEGER INCX,LDA,N
                #       CHARACTER DIAG,TRANS,UPLO
                #       .. Array Arguments ..
                #       DOUBLE PRECISION A(LDA,*),X(*)
        function tbsv!(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar,
                        m::Int, k::Int, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            require_one_based_indexing(A, x)
            n = size(A,2)
            size(A,1) ≥ k+1 || throw(ArgumentError("triangular banded data missing"))
            n == m || throw(DimensionMismatch("matrix is not square: dimensions are $n, $m"))
            if n != length(x)
                throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, diag, m, k,
                 A, max(1,stride(A,2)), x, stride(x, 1))
            x
        end
        function tbsv(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar,
                        m::Int, k::Int, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            tbsv!(uplo, trans, diag, m, k, A, copy(x))
        end
    end
end


text:
"""
    degree_out(N::AbstractEcologicalNetwork)
"""
function degree_out(N::AbstractEcologicalNetwork)
  d_o = vec(sum(N.edges, dims=2))
  return Dict(zip(species(N; dims=1), d_o))
end

"""
    degree_in(N::AbstractEcologicalNetwork)
"""
function degree_in(N::AbstractEcologicalNetwork)
  d_i = vec(sum(N.edges, dims=1))
  return Dict(zip(species(N; dims=2), d_i))
end

"""
    degree(N::AbstractEcologicalNetwork; dims::Union{Nothing,Integer}=nothing)

Returns the degrees of nodes in a network; `dims` can be `1` for out degree, or
`2` for in degree.

#### References

Delmas, E., Besson, M., Brice, M.-H., Burkle, L.A., Dalla Riva, G.V., Fortin,
M.-J., Gravel, D., Guimarães, P.R., Hembry, D.H., Newman, E.A., Olesen, J.M.,
Pires, M.M., Yeakel, J.D., Poisot, T., 2018. Analysing ecological networks of
species interactions. Biological Reviews 112540.
https://doi.org/10.1111/brv.12433

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468

Williams, R.J., 2011. Biology, Methodology or Chance? The Degree Distributions
of Bipartite Ecological Networks. PLoS One 6, e17645.
https://doi.org/10.1371/journal.pone.0017645
"""
function degree(N::AbstractEcologicalNetwork; dims::Union{Nothing,Integer}=nothing)
  dims == 1 && return degree_out(N)
  dims == 2 && return degree_in(N)
  if isnothing(dims)
    if typeof(N) <: AbstractBipartiteNetwork
      return merge(degree_out(N), degree_in(N))
    else
      din = degree_in(N)
      dout = degree_out(N)
      for k in keys(din)
        din[k] += dout[k]
      end
      return din
    end
  end
  throw(ArgumentError("dims can only be 1 (out-degre) or 2 (in-degree) or `nothing` (both), you used $(dims)"))
end

function degree(N::QuantitativeNetwork; dims::Union{Nothing,Integer}=nothing)
  degree(convert(BinaryNetwork, N), dims=dims)
end

"""
**Variance in the outgoing degree**

    degree_out_var(N::ProbabilisticNetwork)
"""
function degree_out_var(N::ProbabilisticNetwork)
  d_o_v = vec(mapslices(_additive_bernoulli_variance, convert(Matrix, N.edges); dims=2))
  return Dict(zip(species(N; dims=1), d_o_v))
end

"""
**Variance in the ingoing degree**

    degree_in_var(N::ProbabilisticNetwork)
"""
function degree_in_var(N::ProbabilisticNetwork)
  d_i_v = vec(mapslices(_additive_bernoulli_variance, convert(Matrix, N.edges); dims=1))
  return Dict(zip(species(N; dims=2), d_i_v))
end

"""
    degree_var(N::ProbabilisticNetwork; dims::Union{Nothing,Integer}=nothing)

Variance in the degree of species in a probabilistic network.
"""
function degree_var(N::ProbabilisticNetwork; dims::Union{Nothing,Integer}=nothing)
  dims == 1 && return degree_out_var(N)
  dims == 2 && return degree_in_var(N)
  if isnothing(dims)
    if typeof(N) <: AbstractBipartiteNetwork
      return merge(degree_out_var(N), degree_in_var(N))
    else
      din = degree_in_var(N)
      dout = degree_out_var(N)
      for k in keys(din)
        din[k] += dout[k]
      end
      return din
    end
  end
  throw(ArgumentError("dims can only be 1 (out-degre) or 2 (in-degree) or `nothing` (both), you used $(dims)"))
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MPFI"
version = v"1.5.5"# <-- This is a lie, we're bumping from 1.5.4 to 1.5.5 to create a Julia v1.6+ release with experimental platforms

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://mirror.eu.oneandone.net/linux/distributions/gentoo/gentoo/distfiles/3a/mpfi-1.5.4.tgz", "3b3938595d720af17973deaf727cfc0dd41c8b16c20adc103a970f4a43ae3a56"),
]

# Bash recipe for building across all platforms
script = raw"""
    cd $WORKSPACE/srcdir/mpfi-*

    # The autotools scripts are not provided setup, so we must configure the configure script
    ./autogen.sh

    # Don't build the static library, documentation, or test set
    ./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-gmp=${prefix} --with-mpfr=${prefix} --disable-static
    make -j${nproc} SUBDIRS="src"
    make install SUBDIRS="src"

    # Remove the libtool files that are created
    rm -f ${prefix}/lib/*.la
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmpfi", :libmpfi)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.2.0"),
    Dependency("MPFR_jll", v"4.1.1"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
Copyright Marcus Greiff 2020

Full test suite for the aerial vehicle controller generation
"""

using LinearAlgebra, Test, AerialVehicleControl

recompile()        # Recompile the C-code
runAllTests = true # If defined, will not recompile the C-code

@testset "All tests" begin
    @testset "Math library tests" begin
        include("controller_math/test_math_functions.jl")
        include("controller_math/test_SO3_maps.jl")
        include("controller_math/test_SU2_maps.jl")
    end
    @testset "Attitude controllers" begin
        include("controller_implementations/test_FSF_attitude_utils.jl")
        include("controller_implementations/test_FSF_continuous_SO3.jl")
        include("controller_implementations/test_FSF_robust_SO3.jl")
        include("controller_implementations/test_FSF_continuous_SU2.jl")
        include("controller_implementations/test_FSF_discontinuous_SU2.jl")
        include("controller_implementations/test_FSF_robust_SU2.jl")
        include("controller_implementations/test_FOF_continuous_SO3.jl")
    end
    @testset "General utility functions" begin
        include("controller_utils/test_power_distribution.jl")
    end
    @testset "Matrix math functions" begin
        include("matrix_math/test_matrix_addition_inplace.jl")
        include("matrix_math/test_matrix_addition.jl")
        include("matrix_math/test_matrix_eigvals.jl")
        include("matrix_math/test_matrix_multiplication.jl")
        include("matrix_math/test_matrix_solve_posdef.jl")
        include("matrix_math/test_matrix_subtraction_inplace.jl")
        include("matrix_math/test_matrix_subtraction.jl")
        include("matrix_math/test_matrix_transposition.jl")
    end
end


text:
"""
    change_sense(objective_sense)

Change the objective sense of optimization.
Possible arguments are `MOI.MAX_SENSE` and `MOI.MIN_SENSE`.

If you want to change the objective and sense at the same time, use
[`change_objective`](@ref) instead to do both at once.
"""
change_sense(objective_sense) =
    (_, opt_model) -> set_objective_sense(opt_model, objective_sense)

"""
    change_optimizer(optimizer)

Change the JuMP optimizer used to run the optimization.

This may be used to try different approaches for reaching the optimum, and in
problems that may require different optimizers for different parts, such as the
[`parsimonious_flux_balance_analysis`](@ref).
"""
change_optimizer(optimizer) = (_, opt_model) -> set_optimizer(opt_model, optimizer)

"""
    change_optimizer_attribute(attribute_key, value)

Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer
to the JuMP documentation and the documentation of the specific optimizer for
usable keys and values.
"""
change_optimizer_attribute(attribute_key, value) =
    (_, opt_model) -> set_optimizer_attribute(opt_model, attribute_key, value)

"""
    silence

Modification that disable all output from the JuMP optimizer (shortcut for
`set_silent` from JuMP).
"""
const silence = (_, opt_model) -> set_silent(opt_model)


text:
module TestMCMC
using GaussianProcesses, Distributions
using Test, Random

Random.seed!(1)

@testset "MCMC" begin
    d, n = 1, 20
    ll = rand(d)
    X = 2π * rand(d, n)
    y = randn(n) .+ 0.5

    kern = RQ(1.0, 1.0, 1.0)

    # Just checks that it doesn't crash
    @testset "HMC" begin
        @testset "Without likelihood" begin
            gp = GP(X, y, MeanZero(), kern)
            set_priors!(gp.kernel, [Distributions.Normal(-1.0, 1.0) for i in 1:3])
            global hmc_chain = mcmc(gp)
        end

        @testset "With likelihood" begin
            lik = GaussLik(-1.0)
            gp = GP(X, y, MeanZero(), kern, lik)
            set_priors!(gp.kernel, [Distributions.Normal(-1.0, 1.0) for i in 1:3])
            mcmc(gp)
        end
    end

    @testset "ESS" begin
        gpess = GP(X, y, MeanZero(), kern)
        set_priors!(gpess.kernel, [Distributions.Normal(-1.0, 1.0) for i in 1:3])
        set_priors!(gpess.logNoise, [Distributions.Normal(-1.0, 1.0)])
        global ess_chain = ess(gpess)
    end

end
end


text:
using ADChain
using Test

@testset "ADChain.jl" begin
    # Write your own tests here.
end


text:
# ---
# title: 668. Kth Smallest Number in Multiplication Table
# id: problem668
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search
# link: <https://leetcode.com/problems/kth-smallest-number-in-multiplication-table/description/>
# hidden: true
# ---
#
# Nearly every one have used the [Multiplication
# Table](https://en.wikipedia.org/wiki/Multiplication_table). But could you find
# out the `k-th` smallest number quickly from the multiplication table?
#
# Given the height `m` and the length `n` of a `m * n` Multiplication Table, and
# a positive integer `k`, you need to return the `k-th` smallest number in this
# table.
#
# **Example 1:**
#
#
#
#     Input: m = 3, n = 3, k = 5
#     Output:
#     Explanation:
#     The Multiplication Table:
#     1	2	3
#     2	4	6
#     3	6	9
#
#     The 5-th smallest number is 3 (1, 2, 2, 3, 3).
#
#
# **Example 2:**
#
#
#
#     Input: m = 2, n = 3, k = 6
#     Output:
#     Explanation:
#     The Multiplication Table:
#     1	2	3
#     2	4	6
#
#     The 6-th smallest number is 6 (1, 2, 2, 3, 4, 6).
#
#
# **Note:**
#
#   1. The `m` and `n` will be in the range [1, 30000].
#   2. The `k` will be in the range [1, m * n]
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end