text:
using AnnealedIS
using AdvancedHMC
using Distributions
using DrWatson
using Dates

@with_kw struct SIRParameters
    tmax::Float32 = 15.0
    total_population::Int = 10_000
    i_0_true::Int = 10
    β_true::Float32 = 0.25
    γ::Float32 = 0.25
    neg_bin_dispersion_param::Float32 = 10
    seed::Int = 1234
end

function transition_kernel(prior_sample::T, i) where {T<:Real}
    if (i > 100) && (i <= 700)
        return Normal(0, 0.01)
    elseif (i > 700)
        return Normal(0, 0.001)
    else
        return Normal(0, 0.1)
    end
end

function transition_kernel(prior_sample::NamedTuple, i)
    return map(x -> transition_kernel(x, i), prior_sample)
end

abstract type MHTransitionConfig end

# TODO: Relax type constraint so that we can have other distributions than Normal.
@with_kw struct MHFixedTuple <: MHTransitionConfig
    proposal::NamedTuple{(:i₀, :β),Tuple{Normal{Float64},Normal{Float64}}} = (
        i₀ = Normal(0, 1), β = Normal(0, 0.1)
    )
end

function get_transition_kernel(transition_kernel_params::MHFixedTuple)
    function kernel(prior_sample, i)
        return transition_kernel_params.proposal
    end
    return kernel
end

@with_kw struct MHFunction <: MHTransitionConfig
    kernel_fn::Function
end

function get_transition_kernel(transition_kernel_params::MHFunction)
    return transition_kernel_params.kernel_fn
end

struct MHThreeLevel <: MHTransitionConfig end

function get_transition_kernel(transition_kernel_params::MHThreeLevel)
    return transition_kernel
end

abstract type AnISConfig end

@with_kw struct AnISMHConfig{T<:MHTransitionConfig,U<:RejectionSampler} <: AnISConfig
    transition_kernel_cfg::T
    rejection_type::U = SimpleRejection()
end

function get_anis_alg(
    anis_params::AnISMHConfig, betas
)
    return AnIS(
        betas,
        get_transition_kernel(anis_params.transition_kernel_cfg),
        anis_params.rejection_type
    )
end

struct ISConfig <: AnISConfig end

get_anis_alg(anis_params::ISConfig, betas) = IS()

@with_kw struct AnISHMCConfig{AD,T<:AdvancedHMC.AbstractProposal,U<:RejectionSampler} <: AnISConfig
    proposal::T
    num_steps::Int = 10
    rejection_type::U = SimpleRejection()
end

AnISHMCConfig(args...; kwargs...) = AnISHMCConfig{Turing.Core.ForwardDiffAD{1}}(args...; kwargs...)
function AnISHMCConfig{AD}(;
    proposal::T,
    num_steps::Int=10,
    rejection_type::U=SimpleRejection()
) where {AD,T<:AdvancedHMC.AbstractProposal,U<:RejectionSampler}
    return AnISHMCConfig{AD,T,U}(proposal, num_steps, rejection_type)
end

function get_anis_alg(
    anis_params::AnISHMCConfig{AD,T,U},
    betas
) where {AD,T<:AdvancedHMC.AbstractProposal,U<:RejectionSampler}
    anis_alg = AnISHMC{AD}(
        betas,
        anis_params.proposal,
        anis_params.num_steps,
        anis_params.rejection_type
    )
end

@with_kw struct AnISHMCGranularConfig{
    AD,
    T<:AdvancedHMC.AbstractProposal,
    M<:AdvancedHMC.AbstractMetric,
    U<:RejectionSampler
} <: AnISConfig
    proposals::Array{T}
    metrics::Array{M}
    num_steps::Int = 10
    rejection_type::U = SimpleRejection()
end

AnISHMCGranularConfig(args...; kwargs...) = AnISHMCGranularConfig{Turing.Core.ForwardDiffAD{1}}(args...; kwargs...)
function AnISHMCGranularConfig{AD}(;
    proposals::Array{T,1},
    metrics::Array{M,1}=repeat(
        [AdvancedHMC.UnitEuclideanMetric(1)], length(proposals)
    ),
    num_steps::Int=10,
    rejection_type::U=SimpleRejection()
) where {
    AD,
    T<:AdvancedHMC.AbstractProposal,
    M<:AdvancedHMC.AbstractMetric,
    U<:RejectionSampler
}
    return AnISHMCGranularConfig{AD,T,M,U}(
        proposals, metrics , num_steps, rejection_type)
end

function get_anis_alg(
    anis_params::AnISHMCGranularConfig{AD,T,M,U},
    betas
) where {
    AD,
    T<:AdvancedHMC.AbstractProposal,
    M<:AdvancedHMC.AbstractMetric,
    U<:RejectionSampler
}
    num_betas = length(betas)
    num_metrics = length(anis_params.metrics)
    num_proposals = length(anis_params.proposals)

    @assert num_metrics == num_proposals
    @assert num_metrics == (num_betas - 1)

    anis_alg = AnISHMC{AD}(
        betas,
        anis_params.proposals,
        anis_params.metrics,
        anis_params.num_steps,
        anis_params.rejection_type
    )
end

@with_kw struct SIRExperimentConfig{T<:AnISConfig}
    experiment_name::String = "test"
    num_samples::Int = 10
    num_betas::Int = 101
    anis_params::T = AnISMHConfig(MHFixedTuple(), SimpleRejection())
    seed::Int = 1234
    data_filename::String = "neg_binomial_data.jld2"
    cost_fn_name::Symbol = :cost_fn_default
end

function DrWatson.default_prefix(e::SIRExperimentConfig)
    datestring = Dates.format(Dates.now(), "yyyymmdd_HHMMSS")
    return "$(datestring)_$(e.experiment_name)_$(typeof(e.anis_params))"
end

DrWatson.allignore(::SIRExperimentConfig) = (:data_filename, :experiment_name)

@with_kw struct ISGroundTruthConfig
    num_samples::Int = Int(1e6)
    qi::Distributions.Distribution = truncated(Normal(10, 10), 0, 1_000)
    qb::Distributions.Distribution = truncated(Normal(1.5, 1.5), 0, Inf)
    cost_fn_name::Symbol = :cost_fn_default
    data_filename::String = "neg_binomial_data.jld2"
    seed::Int = 1234
end


text:
module edit

end

ext: jl

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using Printf
using GLMakie
using Oceananigans
using Oceananigans.Units

include("cyclic_interpolate_utils.jl")

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/"

dh = DataDep("near_global_lat_lon",
    "Forcing data for global latitude longitude simulation",
    [path * "bathymetry_lat_lon_128x60_FP32.bin",
     path * "sea_surface_temperature_25_128x60x12.jld2",
     path * "tau_x_128x60x12.jld2",
     path * "tau_y_128x60x12.jld2"]
)

DataDeps.register(dh)

datadep"near_global_lat_lon"

# 2.8125 degree resolution
Nx = 128
Ny = 60
Nz = 18
reference_density = 1035

#####
##### Load forcing files roughly from CORE2 paper
##### (Probably https://data1.gfdl.noaa.gov/nomads/forms/core/COREv2.html)
#####

filename = [:sea_surface_temperature_25_128x60x12, :tau_x_128x60x12, :tau_y_128x60x12]

for name in filename
    datadep_path = @datadep_str "near_global_lat_lon/" * string(name) * ".jld2"
    file = Symbol(:file_, name)
    @eval $file = jldopen($datadep_path)
end

bathymetry_data = Array{Float32}(undef, Nx*Ny)
bathymetry_path = @datadep_str "near_global_lat_lon/bathymetry_lat_lon_128x60_FP32.bin"
read!(bathymetry_path, bathymetry_data)

bathymetry_data = bswap.(bathymetry_data) |> Array{Float64}
bathymetry_data = reshape(bathymetry_data, Nx, Ny)

τˣ = zeros(Nx, Ny, Nmonths)
τʸ = zeros(Nx, Ny, Nmonths)
T★ = zeros(Nx, Ny, Nmonths)

for month in 1:Nmonths
    τˣ[:, :, month] = file_tau_x_128x60x12["tau_x/$month"] ./ reference_density
    τʸ[:, :, month] = file_tau_y_128x60x12["tau_y/$month"] ./ reference_density
    T★[:, :, month] = file_sea_surface_temperature_25_128x60x12["sst25/$month"]
end

times = 0:1days:24*30days
discrete_times = 0:30days:11*30days

#=
i = 80
j = 12 # Southern ocean?
Nt = length(times)
τ_ij = zeros(Nt)

for (m, t) in enumerate(times)
    n₁ = current_time_index(t)
    n₂ = next_time_index(t)
    τ_ij[m] = cyclic_interpolate(τˣ[i, j, n₁], τˣ[i, j, n₂], t)
end

fig = Figure(resolution = (1200, 600))

ax = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Target SST at (i, j) = ($i, $j)")
scatter!(ax, discrete_times ./ day, τˣ[i, j, :], markersize=20, marker=:utriangle, color=:red)
lines!(ax, times ./ day, τ_ij)

display(fig)
=#

#####
##### Analyze forcing
#####

fig = Figure(resolution = (1200, 800))

t = Node(0.0)

T★_continuous(t) = cyclic_interpolate(T★, t)
τˣ_continuous(t) = cyclic_interpolate(τˣ, t)
τʸ_continuous(t) = cyclic_interpolate(τʸ, t)

T★_t = @lift T★_continuous($t)
τˣ_t = @lift τˣ_continuous($t)
τʸ_t = @lift τʸ_continuous($t)

max_τ = maximum(abs, τˣ)
max_T = maximum(T★)
min_T = minimum(T★)

fig = Figure(resolution = (1200, 600))

ax_τˣ = Axis(fig[1, 1], title="East-west wind stress (m² s⁻²)")
hm_τˣ = heatmap!(ax_τˣ, τˣ_t, colorrange=(-max_τ, max_τ), colormap=:balance)
cb_τˣ = Colorbar(fig[1, 2], hm_τˣ)

ax_τʸ = Axis(fig[2, 1], title="North-south wind stress (m² s⁻²)")
hm_τʸ = heatmap!(ax_τʸ, τʸ_t, colorrange=(-max_τ, max_τ), colormap=:balance)
cb_τʸ = Colorbar(fig[2, 2], hm_τʸ)

ax_T = Axis(fig[3, 1], title="Sea surface temperature (ᵒC)")
hm_T = heatmap!(ax_T, T★_t, colorrange=(min_T, max_T), colormap=:thermal)
cb_T = Colorbar(fig[3, 2], hm_T)

title_str = @lift "Earth day = " * prettytime($t)
ax_t = fig[0, :] = Label(fig, title_str)

display(fig)

record(fig, "annual_cycle_forcing.mp4", times, framerate=8) do tn
    t[] = tn
end


text:
module PiCollisions

using Lynx
using Luxor

@info "Running 'PiCollisions'"

Lynx.init("PiCollisions", 400, 400)

function update(dt)
    # draw here
    background("#111")
end

run!(update, await=true)

end # PiCollisions

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
"""
volume(M::AbstractManifold)

Compute the area/volume of the manifold `M`.
"""
function volume end

Manifolds.base_manifold(μ::AbstractMeasure) = base_manifold(basemeasure(μ))
Manifolds.manifold_dimension(μ::AbstractMeasure) = manifold_dimension(base_manifold(μ))
MeasureTheory.testvalue(μ::PrimitiveMeasure) = μ.point

# The Hausdorff measure generalizes the notion of volume or area to a manifold.
# For manifolds embedded in real space, it can be intuitively understoof as the addition
# of an infinitesimal dimension and the limit of the volume in the embedded space as the
# infinitesimal dimension goes to 0.
# A more useful description comes from the area formula. Given an m-
# dimensional manifold M and a Lipschitz map f: ℝᵐ → M ⊆ ℝⁿ that embeds M in ℝⁿ, then the Hausdorff
# measure dμ(x) on M is is the pushforward of the Lebesgue measure through f, given as
# dμ(x) = det(G(x))^(-1/2) dλ(f⁻¹(x)), where G(x) is the matrix representation of the metric tensor
# in M, given by G(x) = J(x)ᵀ J(x), where J is the Jacobian of f.
# This is the un-normalized Hausdorff measure, which integrates to the area/volume of M.
struct Hausdorff{M<:AbstractManifold,P,A} <: PrimitiveMeasure
    manifold::M
    point::P  # for random sampling
    atlas::A  # optional, when the Hausdorff measure is constructed from a Lebesgue measure
end
Hausdorff(M::AbstractManifold, point) = Hausdorff(M, point, nothing)

function Base.show(io::IO, ::MIME"text/plain", μ::Hausdorff)
    return print(io, "Hausdorff", (μ.manifold, μ.point, μ.atlas))
end

# TODO: how to specify density of Hausdorff wrt Lebesgue in MeasureTheory?
function MeasureTheory.basemeasure(
    μ::Hausdorff{<:AbstractManifold{ℝ},P,<:AbstractAtlas}
) where {P}
    return Lebesgue(MeasureTheory.ℝ)^manifold_dimension(μ.manifold)
end

Manifolds.base_manifold(μ::Hausdorff) = μ.manifold

MeasureTheory.sampletype(::Hausdorff{M,P}) where {M,P} = P

MeasureTheory.logdensity(::Hausdorff, x) = zero(eltype(x))

LinearAlgebra.normalize(μ::Hausdorff) = inv(volume(μ.manifold)) * μ

function Base.rand(rng::AbstractRNG, T::Type, μ::WeightedMeasure{S,<:Hausdorff}) where {S}
    return rand(rng, T, μ.base)
end

# The Haar measure is a left- or right-invariant measure on a group manifold.
# That is, given p ∈ G and q ∈ H ⊆ G, where τ(H, p) is the set of all q τ-translated
# (left- or right-) by p, then μ(τ(H, p)) = μ(H). This is the τ-invariant Haar measure
# on G.
struct Haar{M<:AbstractManifold,P,D<:ActionDirection} <: PrimitiveMeasure
    group::M
    point::P
    direction::D
end
Haar(G::AbstractManifold, point) = Haar(G, point, LeftAction())

const LeftHaar{M,P} = Haar{M,P,LeftAction}
LeftHaar(G, point) = Haar(G, point, LeftAction())

const RightHaar{M,P} = Haar{M,P,RightAction}
RightHaar(G, point) = Haar(G, point, RightAction())

function Base.show(io::IO, ::MIME"text/plain", μ::Haar)
    return print(io, "Haar", (μ.group, μ.point, μ.direction))
end

Manifolds.base_manifold(μ::Haar) = μ.group
Manifolds.base_group(μ::Haar) = μ.group
Manifolds.direction(μ::Haar) = μ.direction

MeasureTheory.sampletype(::Haar{M,P}) where {M,P} = P

MeasureTheory.logdensity(::Haar, x) = zero(eltype(x))


text:
# Based on Quantum Monte Carlo Methods
# https://www.amazon.com/Quantum-Monte-Carlo-Methods-Algorithms/dp/1107006422
# Chapter 3.4

################################################################################
### For generic methods
################################################################################

_reliable_level(B::FullBinner) = max(1, fld(length(B), 32))

_eachlevel(B::FullBinner) = 1:_reliable_level(B)

################################################################################
### Statistics
################################################################################

"""
    mean(B::LogBinner[, lvl])

Calculates the mean for a given level in the Binning Analysis.
"""
mean(B::FullBinner, binsize = 1) = @views mean(B.x[1:fld(end, binsize)*binsize])

"""
    var(B::LogBinner[, lvl])

Calculates the variance of a given level in the Binning Analysis.
"""
var(B::FullBinner, binsize = _reliable_level(B)) = _var(B.x, binsize)

"""
    varN(B::LogBinner[, lvl])

Calculates the variance/N of a given level in the Binning Analysis.
"""
function varN(B::FullBinner, binsize = _reliable_level(B))
    return _var(B.x, binsize) / fld(length(B.x), binsize)
end

# This uses Welford following the implementation in log binner
function _var(A::Vector{T}, binsize) where {T <: Number}
    δ = m1 = m2 = zero(T)
    N = 0
    for i in 1 : binsize : length(A) - binsize + 1
        blockmean = @views mean(A[i : i+binsize-1])
        δ = blockmean - m1
        N += 1
        m1 += δ / N
        m2 += _prod(δ, blockmean - m1)
    end

    return __var(m2, N)
end

function _var(A::Vector{T}, binsize) where {T <: AbstractArray}
    x = A[1]
    δ  = zeros(eltype(x), size(x))
    m1 = zeros(eltype(x), size(x))
    m2 = zeros(eltype(x), size(x))
    blockmean = zeros(eltype(x), size(x))
    N = 0
    for i in 1 : binsize : length(A) - binsize + 1
        blockmean .= zero(eltype(x))
        for j in i : i+binsize-1
            blockmean .+= A[j]
        end
        blockmean ./= binsize
        @. δ = blockmean - m1
        N += 1
        @. m1 += δ / N
        @. m2 += _prod(δ, blockmean - m1)
    end

    return __var.(m2, N)
end

__var(m2::Real, N) = m2 / (N - 1)
__var(m2::Complex, N) = (real(m2) + imag(m2)) / (N - 1)

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
##
# The default processor core.

type ThriftHandler{I,O}
    name::AbstractString
    fn::Function
    intyp::Type{I}
    outtyp::Type{O}
end

type ThriftProcessor
    handlers::Dict{AbstractString, ThriftHandler}
    use_spawn::Bool
    extends::ThriftProcessor
    ThriftProcessor() = (o=new(); o.use_spawn=false; o.handlers=Dict{AbstractString, ThriftHandler}(); o)
end

handle(p::ThriftProcessor, handler::ThriftHandler) = (p.handlers[handler.name] = handler; nothing)
extend(p::ThriftProcessor, extends::ThriftProcessor) = (setfield!(p, :extends, extends); nothing)
distribute(p::ThriftProcessor, use_spawn::Bool=true) = (setfield!(p, :use_spawn, use_spawn); nothing)

function _reply(outp::TProtocol, name::AbstractString, seqid::Int32, mtyp::Int32, m::Any)
    @logmsg("_reply $name:$seqid $m")
    writeMessageBegin(outp, name, mtyp, seqid)
    write(outp, m)
    writeMessageEnd(outp)
    flush(outp.t)
    nothing
end

_exception(extyp::Int32, exmsg::AbstractString, outp::TProtocol, name::AbstractString, seqid::Int32) = _reply(outp, name, seqid, MessageType.EXCEPTION, TApplicationException(extyp, exmsg))

function process(p::ThriftProcessor, inp::TProtocol, outp::TProtocol)
    @logmsg("process begin")
    (name, typ, seqid) = readMessageBegin(inp)

    haskey(p.handlers, name) && (return _process(p, inp, outp, name, typ, seqid))

    isdefined(p, :extends) && (return _process(p.extends, inp, outp, name, typ, seqid))

    skip(inp, TSTRUCT)
    readMessageEnd(inp)
    _exception(ApplicationExceptionType.UNKNOWN_METHOD, "Unknown function $name", outp, name, seqid)
end

function _process(p::ThriftProcessor, inp::TProtocol, outp::TProtocol, name::AbstractString, typ::Int32, seqid::Int32)
    handler = p.handlers[name]
    @logmsg("_process: reading instruct of type $(handler.intyp)")
    instruct = read(inp, handler.intyp)
    readMessageEnd(inp)
    @logmsg("_process: calling handler function")
    if p.use_spawn
        outstruct = fetch(@spawn handler.fn(instruct))
    else
        outstruct = handler.fn(instruct)
    end
    @logmsg("_process: out of handler function. return val: $outstruct")
    if !isa(outstruct, handler.outtyp)
        _exception(ApplicationExceptionType.MISSING_RESULT, "Invalid return type. Expected $(handler.outtyp). Got $(typeof(outstruct))", outp, name, seqid)
        return
    end
    isa(outstruct, Void) || _reply(outp, name, seqid, MessageType.REPLY, outstruct)
end


text:
@timedtestset "test A_canonical" for ix in 1:10
    psi = bimps(rand, 9, 4)
    T = mpo_triangular_AF_ising_alternative()

    _, psi_a = A_canonical(T, psi)

    # test psi_a.A is right canonical
    A_tmp = permute(psi_a.A, (1, ), (2, 3))
    @test A_tmp * A_tmp' ≈ id(ℂ^9)
    # fidelity
    @test abs(ln_fidelity(psi_a.A, psi.A)) < 1e-14
    # test the fixed point equation
    @tensor lhs[-1, -2; -3] := psi_a.B[1, 3, 4] * T[-2, 5, 2, 3] * psi_a.A[-1, 2, 1] * psi_a.A'[4, -3, 5]
    @test lhs ≈ (lhs[1] / psi_a.B[1]) * psi_a.B
end

@timedtestset "test B_canonical" for ix in 1:10
    psi = bimps(rand, 9, 4)
    T = mpo_triangular_AF_ising_alternative()

    _, psi_b = B_canonical(T, psi)

    # test psi_b.B is left canonical
    @test psi_b.B' * psi_b.B ≈ id(ℂ^9)
    # fidelity
    @test abs(ln_fidelity(psi_b.B, psi.B)) < 1e-14
    # test the fixed point equation
    @tensor lhs[-1, -2; -3] := psi_b.A[4, 2, 1] * T[5, -2, 2, 3] * psi_b.B[1, 3, -3] * psi_b.B'[-1, 4, 5]
    @test lhs ≈ (lhs[1] / psi_b.A[1]) * psi_b.A
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct FileInput <: Check
        f1::String
        f2::String
        colnames::Array{Symbol,1}
    end

Information for files to compare.

# Arguments
- `f1::String`: First file to compare
- `f2::String`: Second file to compare
- `colnames::Array{Symbol,1}`: column names shared by the two files
"""
mutable struct FileInput <: Check
    "First file to compare"
    f1::String
    "Second file to compare"
    f2::String
    "column names shared by the two files"
    colnames::Array{Symbol,1}
end

function FileInput(; f1, f2, colnames, )
    FileInput(f1, f2, colnames, )
end

"""Get [`FileInput`](@ref) `f1`."""
get_f1(value::FileInput) = value.f1
"""Get [`FileInput`](@ref) `f2`."""
get_f2(value::FileInput) = value.f2
"""Get [`FileInput`](@ref) `colnames`."""
get_colnames(value::FileInput) = value.colnames

"""Set [`FileInput`](@ref) `f1`."""
set_f1!(value::FileInput, val) = value.f1 = val
"""Set [`FileInput`](@ref) `f2`."""
set_f2!(value::FileInput, val) = value.f2 = val
"""Set [`FileInput`](@ref) `colnames`."""
set_colnames!(value::FileInput, val) = value.colnames = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# Define convolutional layer:
mutable struct Conv; w; b; f; end
(c::Conv)(x) = c.f.(pool(conv4(c.w, x) .+ c.b))
Conv(w1,w2,cx,cy,f=relu; atype=_atype) = Conv(
    param(w1,w2,cx,cy; atype=atype), param0(1,1,cy,1; atype=atype), f)

# Define dense layer:
mutable struct Dense; w; b; f; end
(d::Dense)(x) = d.f.(d.w * mat(x) .+ d.b)
Dense(i::Int,o::Int,f=relu; atype=_atype) = Dense(
    param(o,i; atype=atype), param0(o; atype=atype), f)

# Define a chain of layers and a loss function:
mutable struct Chain; layers; end
(c::Chain)(x) = (for l in c.layers; x = l(x); end; x)
(c::Chain)(x,y) = nll(c(x),y)

"""
    Lenet(; atype)

Returns a randomly initialized Lenet model. Used for demo and testing.
Take a look at for more detailed information. atype is KnetArray{Float32}
if GPU is detected, otherwise Array{Float64}.
"""
Lenet(; atype=_atype) = Chain(
    (Conv(5,5,1,20; atype=atype),
     Conv(5,5,20,50; atype=atype),
     Dense(800,500; atype=atype),
     Dense(500,10,identity; atype=atype)))


text:
# https://www.math.wisc.edu/~anderson/papers/AndPostleap.pdf
# Incorporating postleap checks in tau-leaping
# J. Chem. Phys. 128, 054103 (2008); https://doi.org/10.1063/1.2819665
function cpp_bridge(cpp,W,W0,Wh,q,h,u,p,t,rng)
  rand.(rng,Distributions.Binomial.(Int.(Wh),float.(q)))
end
function cpp_bridge!(rand_vec,cpp,W,W0,Wh,q,h,u,p,t,rng)
  rand_vec .= rand.(rng,Distributions.Binomial.(Int.(Wh),float.(q)))
end

mutable struct CompoundPoissonProcess{R,CR}
  rate::R
  currate::CR
  computerates::Bool
  function CompoundPoissonProcess(rate,t0,W0;computerates=true,kwargs...)
    cpp = new{typeof(rate),typeof(W0)}(rate,W0,computerates)
    NoiseProcess{false}(t0,W0,nothing,cpp,(W,W0,Wh,q,h,u,p,t,rng)->cpp_bridge(cpp,W,W0,Wh,q,h,u,p,t,rng);continuous=false,cache=cpp,kwargs...)
  end
end
function (P::CompoundPoissonProcess)(W,dt,u,p,t,rng)
  P.computerates && (P.currate = P.rate(u,p,t))
  PoissonRandom.pois_rand.(rng,dt.*P.currate)
end

struct CompoundPoissonProcess!{R,CR}
  rate::R
  currate::CR
  computerates::Bool
  function CompoundPoissonProcess!(rate,t0,W0;computerates=true,kwargs...)
    cpp = new{typeof(rate),typeof(W0)}(rate,copy(W0),computerates)
    NoiseProcess{true}(t0,W0,nothing,cpp,(rand_vec,W,W0,Wh,q,h,u,p,t,rng)->cpp_bridge!(rand_vec,cpp,W,W0,Wh,q,h,u,p,t,rng);continuous=false,cache=cpp,kwargs...)
  end
end
function (P::CompoundPoissonProcess!)(rand_vec,W,dt,u,p,t,rng)
  P.computerates && P.rate(P.currate,u,p,t)
  @. rand_vec = PoissonRandom.pois_rand(rng,dt*P.currate)
end


text:
using ForwardDiff

"""
dict = make_dict(argstrings, x, [starting_dict=Dict()] )

Given a list of strings, and a list of values, makes a dictionary of Symbols to values, with the Symbols
corresponding to each of the strings.  Mostly used to pass arguments as a keyword-value set into a function.
If one of the elements of argstrings is *not* a string, but is instead a 2-long list, the first element of that
list should be a string, and the second element of that list should be a positive integer. This will be
interpreted as "don't take only one value, take this number of values and this parameter will be a vector"

# PARAMS:

* argstrings     A list of strings. Each element may also be a two-long list of a string, positive integer, e.g., ["this" 3]

* x              A vector of numeric values. Its length must be such that all the strings in argstrings
                 can take their corresponding element(s), sequentially, from x

* starting_dict  An optional initial dictionary to work with.  Any key in this starting dictionary matching an argstring
                 will be replaced by the new value. Keys not matched will remain.

# RETURNS:

dict             The symbol dictionary.

# EXAMPLES:

>> make_dict(["this", "that", ["there", 2]], [10, 20, 3, 4])

Dict{Any,Any} with 3 entries:
  :this  => 10
  :that  => 20
  :there => [3,4]

>> make_dict(["doo", "gaa"], [10, 20], Dict(:blob=>100, :gaa=>-44))

Dict{Symbol,Int64} with 3 entries:
  :gaa  => 20
  :blob => 100
  :doo  => 10

"""
function make_dict(args, x, starting_dict=Dict())
    # For error diagnostics, check that the length of the param vector specified in args matches the length of x
    nargs = 0
    for i in [1:length(args);]
        if typeof(args[i])==String # if the entry in args is a string, then there's one corresponding scalar entry in x0
            nargs += 1
        else
            nargs += args[i][2]    # otherwise, the entry in args should be a  [varnamestring, nvals] vector,
            # indicating that the next nvals entries in x0 are all a single vector, belonging to variable
            # with name varnamestring.
        end
    end
    if nargs != length(x)
        error("Oy! args and x must indicate the same total number of variables!")
    end

    # ---- done error-checking, now main function

    kwargs = starting_dict;
    i = 1; j=1
    while i<=length(args)
        if typeof(args[i])==String
            kwargs = merge(kwargs, Dict(Symbol(args[i]) => x[j]))
        else
            if length(args[i]) == 2
                extra = args[i][2]-1
                kwargs = merge(kwargs, Dict(Symbol(args[i][1]) => x[j:(j+extra)]))
                j = j+extra
            else
                error("Each element of the args vector must be either a string, or a 2-long vector, first element a string, second integer")
            end
        end
        i = i+1; j=j+1
    end
    return kwargs
end

# Here we're going to define a closure over x so that when this code runs, it sets the local variable
# x to report ForwardDiff's verison number; then we export the function FDversion, that simply returns
# x.  When we call FDversion(), it simply returns the value of x, stored locally inside the let block.
# So it is extremely fast.
let x
    try
        x = Pkg.installed("ForwardDiff").major + 0.1*Pkg.installed("ForwardDiff").minor
    catch
        error("Is ForwardDiff really installed???")
    end

    global FDversion

    @doc """
    v = FDversion()

    Return the installed version of ForwardDiff as a floating point, major.minor
    """ function FDversion()
        return x
    end
end

if FDversion() < 0.6
    # --------------------------------------------------------------
    #
    #               FOR FORWARDDIFF < 0.6   (Julia 0.5.2)
    #
    # --------------------------------------------------------------

    # """
    # We define functions to convert Duals, the variable types used by ForwardDiff,
    # to Floats. This is useful if we want to print out the value of a variable
    # (since print doesn't know how to Duals). Note that after being converted to a Float, no
    # differentiation by ForwardDiff can happen!  e.g. after
    #     x = convert(Float64, y)
    # ForwardDiff can still differentiate y, but it can't differentiate x
    # """

    import Base.convert
    convert(::Type{Float64}, x::ForwardDiff.Dual) = Float64(x.value)
    function convert(::Array{Float64}, x::Array{ForwardDiff.Dual})
        y = zeros(size(x));
        for i in 1:prod(size(x))
            y[i] = convert(Float64, x[i])
        end
        return y
    end

    """
    function M = ForwardDiffZeros(m, n; nderivs=0, difforder=0)

    Use instead of zeros(). Creates a matrix of zeros, of size m rows by n columns, with elements appropriate for
    differentiation by ForwardDiff. If nderivs==0 or difforder==0 then the elements will be regular
    Float64, not ForwardDiff types.

    PARAMETERS:
    ===========

    m        Integer, number of rows

    n        Integer, number of columns

    OPTIONAL PARAMETERS:
    ====================

    nderivs=0       The number of variables that we'll be differentiating with respect to. In other
                    words, this number is equal to the length of the gradient. If this is left as zero (the default) then
                    the data type will be regular Float64

    difforder=0     The order of the derivative we will want to take.  Zero means nothing, stick with
                    regular Float64, 1 means gradient, 2 means hessian

    RETURNS:
    ========

    An m-by-n matrix of zeros that can be used with Forward Diff.

    """
    function ForwardDiffZeros(m, n; nderivs=0, difforder=0)
        if nderivs == 0 || difforder == 0
            return zeros(m, n)
        elseif difforder == 1
            return zeros(ForwardDiff.Dual{nderivs, Float64}, m , n)
        elseif difforder == 2
            return zeros(ForwardDiff.Dual{nderivs, ForwardDiff.Dual{nderivs, Float64}}, m, n)
        else
            error("Don't know how to do that order of derivatives!", nderivs)
        end
    end

else

    # --------------------------------------------------------------
    #
    #         FOR FORWARDDIFF >= 0.6   (Julia 0.6 and onwards)
    #
    # --------------------------------------------------------------

    @doc """
    e = get_eltype(vars)

    vars should be a tuple of variables. If any of them is a ForwardDiff Dual, this function returns
    the typeof of that one (the first one encountered); otherwise it returns Float64.
    """ function get_eltype(vars)
        # print("vars is "); print(vars); print("\n")
        # print("typeof(vars) is "); print(typeof(vars)); print("\n")
        if ! (typeof(vars)<:Tuple)
            error("vars must be a Tuple (of variables)")
        end
        for v in vars
            # print("About to check v = "); print(v); print("\n")
            # If it's an array, check its elements
            if typeof(v)<:Array && length(v)> 0 && typeof(v[1])<:ForwardDiff.Dual
                return typeof(v[1])
            end
            # If it's an array of tuples, check each one
            if typeof(v)<:Array && length(v)> 0 && typeof(v[1])<:Tuple
                for tup in v; if typeof(tup[2])<:ForwardDiff.Dual; return typeof(tup[2]); end; end
            end
            # If it's a Pair, check its value
            if typeof(v)<:Pair && typeof(v[2])<:ForwardDiff.Dual
                return typeof(v[2])
            end
            # If it's a tuple, try to turn it into a Dict
            if typeof(v)<:Tuple
                try; myv = Dict(v); catch; error("Sorry don''t know how to deal with that kind of Tuple"); end
            else
                myv = v
            end
            # If it's a Dict, check all the key contents
            if typeof(myv)<:Dict
                for k in keys(myv)
                    # print("about to check key "); print(k); print(" with value "); print(myv[k]); print("\n")
                    if get_eltype(Tuple(myv[k]))<:ForwardDiff.Dual
                        return get_eltype(Tuple(myv[k]))
                    end
                end
            end
            if typeof(myv)<:ForwardDiff.Dual
                return typeof(myv)
            end
        end
        return Float64
    end

end

# We're deprecating automatic converts from Dual to Float64;
# instead, let's use an explicit call to a get_value function.
# That means that any accidental casts into Float64s will cause errors, alerting us of
# problems rather than letting the code run but producing derivatives that are zero

@doc """
v = get_value(x)

If you're going to @print something that might be a ForwardDiff Dual, use this function. It'll
return the value of the number if it is a Dual and just the number if it was not a Dual, suitable for
printing, e.g.

    @printf("%g\n", get_value(x))

will work regardless of whether x is a ForwardDiff Dual, a Float64, or an Int64
""" function get_value(x)
    if typeof(x)<:Array && length(x)>0
        if typeof(x[1])<:ForwardDiff.Dual
            y = zeros(size(x));
            if typeof(x[1].value)<:ForwardDiff.Dual  # nested, taking 2nd derivative
                for i=1:length(y); y[i] = x[i].value.value; end;
            else # not nested, 1st derivative
                for i=1:length(y); y[i] = x[i].value; end;
            end
            return y
        elseif typeof(x[1])<:Int64 || typeof(x[1])<:Float64
            return x
        else
            error(@sprintf("Don't know how to get the value of type %s", typeof(x[1])))
        end
    elseif typeof(x)<:Array
        return zeros(size(x))
    elseif typeof(x)<:ForwardDiff.Dual
        if typeof(x.value)<:ForwardDiff.Dual; return x.value.value; else; return x.value; end;
    elseif typeof(x)<:Int64 || typeof(x)<:Float64
        return x
    else
        error(@sprintf("Don't know how to get the value of type %s", typeof(x)))
    end
end

if FDversion() < 0.6
    # --------------------------------------------------------------
    #
    #               FOR FORWARDDIFF < 0.6   (Julia 0.5.2)
    #
    # --------------------------------------------------------------

    using DiffBase
    @doc """
    function value, gradient, hessian = vgh(func, x0)

    Wrapper for ForwardDiff.hessian!() that computes and returns all three of a function's value, gradient, and hessian.

    EXAMPLE:
    ========

    function tester(x::Vector)

        return sum(x.*x)
    end

    value, grad, hess = vgh(tester, [10, 3.1])
    """ function vgh(func, x0, do_hess = true)
        if do_hess
            out = DiffBase.HessianResult(x0)
            ForwardDiff.hessian!(out, func, x0)
            value    = DiffBase.value(out)
            gradient = DiffBase.gradient(out)
            hessian  = DiffBase.hessian(out)
        else
            out = DiffBase.GradientResult(x0)
            ForwardDiff.gradient!(out, func, x0)
            value    = DiffBase.value(out)
            gradient = DiffBase.gradient(out)
            hessian  = []
        end
        return value, gradient, hessian
    end

    @doc """
    function value, gradient, hessian = keyword_vgh(func, args, x0)

    Wrapper for vgh() that computes and returns all three of a function's value, gradient, and hessian, but now
    uses make_dict() to apply it to a function that only takes keyword-value pairs.

    *Note that func MUST also take the keyword parameters nderivs and difforder*. If you declare any vectors or
    matrices inside func() (or inside any function inside func()), use ForwardDiffZeros with these two parameters,
    do NOT use zeros(). Your gradients will come out as zero is you use zeros().

    # PARAMETERS

    * func    A function that takes keyword-value pairs only, including nderivs and difforder.  I.e., it must be a function declared as `function func(; nderivs=0, difforder=0, other_kw_value_pairs)` or as `function func(; nderivs=0, difforder=0, other_kw_value_pairs_dict...)`
    * args    A list of strings indicating names of variables to work with
    * x0      A vector with the value of the variables indicates in args.  **See make_dict() for how to pass both scalars and vectors as variables**

    # IMPORTANT JULIA BUG

    If you modify func, it is possible that keyword_vgh() will still work on the previously defined version. AACK!
    That's horrible! Alice Yoon's tip on the workaround: instead of func(), use (;params...) -> func(; params...) and then
    everything will be fine. Perhaps this bug will be fixed in Julia 0.6

    # EXAMPLE:

    function tester(;a=10, b=20, c=30, nderivs=0, difforder=0)
        M = ForwardDiffZeros(3, 3; nderivs=nderivs, difforder=difforder)
        M[1,1] = a^2*10
        M[2,2] = b*20
        M[3,3] = a*sqrt(c)*30.1
        return trace(M)
    end

    value, grad, hess = keyword_vgh(tester, ["a", "c"], [10, 3.1])

    value, grad, hess = keyword_vgh((;params...) -> tester(;params...), ["a", "c"], [10, 3.1])

    """ function keyword_vgh(func, args, x0, do_hess = true)

        value, gradient, hessian = vgh(x -> func(;nderivs=length(x), difforder=2, make_dict(args, x)...), x0, do_hess)

        return value, gradient, hessian
    end

else
    # --------------------------------------------------------------
    #
    #         FOR FORWARDDIFF >= 0.6   (Julia 0.6 and onwards)
    #
    # --------------------------------------------------------------

    using DiffResults

    @doc """
    function value, gradient, hessian = vgh(func, x0)

    Wrapper for ForwardDiff.hessian!() that computes and returns all three of a function's value, gradient, and hessian.

    EXAMPLE:
    ========

    function tester(x::Vector)

        return sum(x.*x)
    end

    value, grad, hess = vgh(tester, [10, 3.1])
    """ function vgh(func, x0, do_hess=true)
        if do_hess
            out = DiffResults.HessianResult(x0)
            out = ForwardDiff.hessian!(out, func, x0)
            value    = DiffResults.value(out)
            gradient = DiffResults.gradient(out)
            hessian  = DiffResults.hessian(out)
        else
            out = DiffResults.GradientResult(x0)
            out = ForwardDiff.gradient!(out, func, x0)
            value    = DiffResults.value(out)
            gradient = DiffResults.gradient(out)
            hessian  = []
        end
        return value, gradient, hessian

    end

    @doc """
    function value, gradient, hessian = keyword_vgh(func, args, x0)

    Wrapper for vgh() that computes and returns all three of a function's value, gradient, and hessian, but now
    uses make_dict() to apply it to a function that only takes keyword-value pairs.

    *Note that if you declare any vectors or matrices inside func() (or inside any function inside func()),
    you will need to make sure they are ForwardDiff Duals if you want to differentiate w.r.t. their contents.*
    The function get_eltype()  can help you with this.  For example, if you have three variables, a, b, and c,
    and you don't know in advance which one you will differentiate w.r.t., you could declare a matrix using

        new_matrix = zeros(get_eltype((a,b,c)), 2, 3)

    and that will make sure that it is the right type if you later want derivatives of contents of that
    matrix.

    # PARAMETERS

    * func    A function that takes keyword-value pairs only, including nderivs and difforder.  I.e., it must be a function declared as `function func(; nderivs=0, difforder=0, other_kw_value_pairs)` or as `function func(; nderivs=0, difforder=0, other_kw_value_pairs_dict...)`
    * args    A list of strings indicating names of variables to work with
    * x0      A vector with the value of the variables indicates in args.  **See make_dict() for how to pass both scalars and vectors as variables**

    # IMPORTANT JULIA BUG

    If you modify func, it is possible that keyword_vgh() will still work on the previously defined version. AACK!
    That's horrible! Alice Yoon's tip on the workaround: instead of func(), use (;params...) -> func(; params...) and then
    everything will be fine. Not sure yet whether this bug is fixed in Julia 0.6

    # EXAMPLE:

    function tester(;a=10, b=20, c=30)
        M = zeros(get_eltype((a,b,c)), 3, 3)
        M[1,1] = a^2*10
        M[2,2] = b*20
        M[3,3] = a*sqrt(c)*30.1
        return trace(M)
    end

    value, grad, hess = keyword_vgh(tester, ["a", "c"], [10, 3.1])

    value, grad, hess = keyword_vgh((;params...) -> tester(;params...), ["a", "c"], [10, 3.1])

    """ function keyword_vgh(func, args, x0, do_hess=true)

        value, gradient, hessian = vgh(x -> func(;make_dict(args, x)...), x0, do_hess)

        return value, gradient, hessian
    end

end


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 15, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function solve(::Puzzle{2020, 15, 1}, input)
    nothing
end

function solve(::Puzzle{2020, 15, 2}, input)
    nothing
end

@testset "Day 15" begin
    TEST_INPUT_1 = """
"""

    @test_skip(solve(Puzzle(15, 1), TEST_INPUT_1) == nothing)
    @test_skip(solve(Puzzle(15, 1), input(Puzzle(15, 1))) == nothing)

    TEST_INPUT_2 = """
"""

    @test_skip(solve(Puzzle(15, 2), TEST_INPUT_2) == nothing)
    @test_skip(solve(Puzzle(15, 2), input(Puzzle(15, 2))) == nothing)
end


text:
module SumBassSolver
using DataStructures, Ipopt, JuMP, Roots, SigmoidalProgramming, Optim

export Bass, Bass_derivative, Bass_inverse, Bass_derivative_inverse, Bass_inflection
export sumBassSolver
export sumBassHomogeneousSolver, sumBassDiffStartSolver
export sumSigmoidSolverSB

include("Bass.jl")
include("UdellBoyd.jl")
include("SrivastavaBullo.jl")
include("LocalOpt.jl")
include("HomogeneousSolver.jl")
include("DiffstartSolver.jl")
include("SrivastavaBulloGen.jl")

@doc raw"""
    sumBassSolver(p,q,m,c0,alpha,F,T,n,Pi;alg="UB", localopt=false)

This function maximizes the sum Bass functions and returns the maximizer
and maximum. `alg` is one of UB (Udell and Boyd's algorithm) or
SB (Srivastava and Bullo's algorithm). `localopt` indicates making a local search
after the algorithms terminate.
"""
function sumBassSolver(p,q,m,c0,alpha,F,T,Pi::Real;alg = "SB",localopt=false)
    n = length(p) # Number of sites
    if sum(F) > Pi
        throw(DomainError(Pi,"the problem is unfeasible because the capacity is too small."))
    end
    if n==0
        throw(DomainError(p,"inputs need to be of size at least 1"))
    end
    if length(q) != n || length(m) != n || length(c0) != n || length(alpha) != n || length(F) != n || length(T) != n
        throw(DomainError("p,q,m,c0,alpha,F,T","inputs must be vectors of the same length."))
    end
    if sum(c0 .> m) > 0
        throw(DomainError("c0,m","existing clients needs to be smaller than the market potential."))
    end

    if alg == "UB"
        maximizer,maximum = sumBassUB(p,q,m,c0,alpha,F,T,Pi)
    elseif alg == "SB"
        maximizer,maximum = sumBassSB(p,q,m,c0,alpha,F,T,Pi)
    else
        throw(DomainError(alg,"alg is one of UB (Udell and Boyd's algorithm) or SB (Srivastava and Bullo's algorithm)."))
    end

    if localopt
        maximizer,maximum = sumBassLocal(p,q,m,c0,alpha,F,T,Pi,maximizer)
    end

    return maximizer,maximum
end

end

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
# This supports nonbroadcasting math on NamedDimArrays

################################################
# Matrix product

valid_matmul_dims(a::Tuple{Symbol}, b::Tuple{Vararg{Symbol}}) = true
function valid_matmul_dims(a::Tuple{Symbol,Symbol}, b::Tuple{Vararg{Symbol}})
    a_dim = a[end]
    b_dim = b[1]

    return a_dim === b_dim || a_dim === :_ || b_dim === :_
end

matmul_names((a1, a2)::Tuple{Symbol,Symbol}, (b,)::Tuple{Symbol}) = (a1,)
matmul_names((a1, a2)::Tuple{Symbol,Symbol}, (b1, b2)::Tuple{Symbol,Symbol}) = (a1, b2)
matmul_names((a1,)::Tuple{Symbol}, (b1, b2)::Tuple{Symbol,Symbol}) = (a1, b2)

function throw_matrix_dim_error(a, b)
    msg = "Cannot take matrix product of arrays with different inner dimension names. $a vs $b"
    return throw(DimensionMismatch(msg))
end

function matrix_prod_names(a, b)
    # 0 Allocations. See `@btime (()-> matrix_prod_names((:foo, :bar),(:bar,)))()`
    valid_matmul_dims(a, b) || throw_matrix_dim_error(a, b)
    res = matmul_names(a, b)
    return compile_time_return_hack(res)
end

for (NA, NB) in ((1, 2), (2, 1), (2, 2))  #Vector * Vector, is not allowed
    @eval function Base.:*(
        a::NamedDimsArray{A,T,$NA}, b::NamedDimsArray{B,S,$NB},
    ) where {A,B,T,S}
        L = matrix_prod_names(A, B)
        data = *(parent(a), parent(b))
        return NamedDimsArray{L}(data)
    end
end

# vector^T * vector
function Base.:*(
    a::NamedDimsArray{A,T,2,<:CoVector}, b::NamedDimsArray{B,S,1},
) where {A,B,T,S}
    valid_matmul_dims(A, B) || throw_matrix_dim_error(last(A), first(B))
    return *(parent(a), parent(b))
end

function Base.:*(a::NamedDimsArray{L,T,2,<:CoVector}, b::AbstractVector) where {L,T}
    return *(parent(a), b)
end

# Using `CoVector` results in Method ambiguities; have to define more specific methods.
for A in (Adjoint{<:Number,<:AbstractVector}, Transpose{<:Real,<:AbstractVector{<:Real}})
    @eval function Base.:*(a::$A, b::NamedDimsArray{L,T,1,<:AbstractVector{T}}) where {L,T}
        return *(a, parent(b))
    end
end

"""
    @declare_matmul(MatrixT, VectorT=nothing)

This macro helps define matrix multiplication for the types
with 2D type parameterization `MatrixT` and 1D `VectorT`.
It defines the various overloads for `Base.:*` that are required.
It should be used at the top level of a module.
"""
macro declare_matmul(MatrixT, VectorT=nothing)
    dim_combos = VectorT === nothing ? ((2, 2),) : ((1, 2), (2, 1), (2, 2))
    codes = map(dim_combos) do (NA, NB)
        TA_named = :(NamedDims.NamedDimsArray{<:Any,<:Any,$NA})
        TB_named = :(NamedDims.NamedDimsArray{<:Any,<:Any,$NB})
        TA_other = (VectorT, MatrixT)[NA]
        TB_other = (VectorT, MatrixT)[NB]

        quote
            function Base.:*(a::$TA_named, b::$TB_other)
                return *(a, NamedDims.NamedDimsArray{dimnames(b)}(b))
            end
            function Base.:*(a::$TA_other, b::$TB_named)
                return *(NamedDims.NamedDimsArray{dimnames(a)}(a), b)
            end
        end
    end
    return esc(Expr(:block, codes...))
end

# The following two methods can be defined by using
# @declare_matmul(Diagonal, AbstractVector)
# but that overwrites existing *(1D NDA, Vector) methods
# should improve the macro above to deal with this case
function Base.:*(a::Diagonal, b::NamedDimsArray{<:Any,<:Any,1})
    return *(NamedDimsArray{dimnames(a)}(a), b)
end

function Base.:*(a::NamedDimsArray{<:Any,<:Any,1}, b::Diagonal)
    return *(a, NamedDimsArray{dimnames(b)}(b))
end

@declare_matmul(AbstractMatrix, AbstractVector)
@declare_matmul(
    Adjoint{<:Any,<:AbstractMatrix{T1}} where {T1}, Adjoint{<:Any,<:AbstractVector}
)
@declare_matmul(Diagonal,)

################################################
# Others

function Base.inv(nda::NamedDimsArray{L,T,2}) where {L,T}
    data = inv(parent(nda))
    names = reverse(L)
    return NamedDimsArray{names}(data)
end

# Statistics
for fun in (:cor, :cov)
    @eval function Statistics.$fun(a::NamedDimsArray{L,T,2}; dims=1, kwargs...) where {L,T}
        numerical_dims = dim(a, dims)
        data = Statistics.$fun(parent(a); dims=numerical_dims, kwargs...)
        names = symmetric_names(L, numerical_dims)
        return NamedDimsArray{names}(data)
    end
end

# CovarianceEstimation
for E in (:LinearShrinkage, :SimpleCovariance, :AnalyticalNonlinearShrinkage)
    @eval function CovarianceEstimation.cov(
        estimator::$E, a::NamedDimsArray{L,T,2}; dims=1, kwargs...
    ) where {L,T}
        numerical_dims = dim(a, dims)
        # cov returns a Symmetric matrix which needs to be rewrapped in a NamedDimsArray
        data = cov(estimator, parent(a); dims=numerical_dims, kwargs...)
        names = symmetric_names(L, numerical_dims)
        return NamedDimsArray{names}(data)
    end
end

function symmetric_names(L::Tuple{Symbol,Symbol}, dims::Integer)
    # 0 Allocations. See `@btime (()-> symmetric_names((:foo, :bar), 1))()`
    names = if dims == 1
        (L[2], L[2])
    elseif dims == 2
        (L[1], L[1])
    else
        (:_, :_)
    end
    return compile_time_return_hack(names)
end


text:
module test_loltools_matchv4

using Test
using LOLTools.MatchV4

@test MatchV4.match_by_tournament_code isa Function
@test MatchV4.match_by_id isa Function

include("helpers.jl")

struct MatchController <: ApplicationController
    conn::Conn
end

teams = [(teamId=0, win=0, firstBlood=0, firstTower=0, firstInhibitor=0, firstBaron=0, firstDragon=0, firstRiftHerald=0, towerKills=0, inhibitorKills=0, baronKills=0, dragonKills=0, vilemawKills=0, riftHeraldKills=0, dominionVictoryScore=0, bans=0)]
match_tup = (seasonId=0, queueId=0, gameId=0, participantIdentities=[(participantId=1, player=(platformId = "KR", accountId = "", summonerName = "", summonerId = "", currentPlatformId = "KR", currentAccountId = "", matchHistoryUri = "/v1/stats/player_history/KR/", profileIcon = 0))], gameVersion=0, platformId=0, gameMode=0, mapId=0, gameType=0, teams=teams, participants=[(participantId = 1, teamId = 100, championId = 516, spell1Id = 4, spell2Id = 12, stats =(), timeline = (participantId = 7,))], gameDuration=0, gameCreation=0)

function match_by_tournament_code(c::MatchController)
    # @info :c c.params.tournamentCode
    render(JSON, match_tup)
end

function match_by_id(c::MatchController)
    # @info :c c.params.matchId
    render(JSON, match_tup)
end

routes() do
    get("/lol/match/v4/matches/by-tournament-code/:tournamentCode/ids", MatchController, match_by_tournament_code)
    get("/lol/match/v4/matches/:matchId", MatchController, match_by_id)
end

api_key = ""
platform = "na1"
MatchV4.match_by_tournament_code(api_key, platform, "tournamentCode"; action=mock_action)
MatchV4.match_by_id(api_key, platform, 0; action=mock_action)

end # module test_loltools_matchv4


text:
function rbc()
    @parameters α β ρ δ σ Ω_1
    @variables k z c q i
    @make_markov k z c q i
    x = [k, z]
    y = [c, q, i]
    p = [α, β, ρ]
    p_f = [δ, σ, Ω_1]

    H = [1 / c - (β / c_p) * (α * exp(z_p) * k_p^(α - 1) + (1 - δ)),
         c + k_p - (1 - δ) * k - q,
         q - exp(z) * k^α,
         z_p - ρ * z,
         i - (k_p - (1 - δ) * k)]

    x̄ = [k_ss ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z_ss ~ 0]
    ȳ = [c_ss ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
         q_ss ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)),
         i_ss ~ δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1))]

    x̄_iv = [z_ss ~ 0, k_ss ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1))]
    ȳ_iv = [q_ss ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)),
            c_ss ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                    δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
            i_ss ~ δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1))]

    n_ϵ = 1
    n_z = 2
    n_x = length(x)
    n_y = length(y)
    n_p = length(p)
    Γ = reshape([σ], n_ϵ, n_ϵ)
    η = reshape([0; -1], n_x, n_ϵ) # η is n_x * n_ϵ matrix

    Q = zeros(n_z, n_y + n_x) # The order is [y, x]
    Q[1, 1] = 1.0 # c
    Q[2, 3] = 1.0 # i

    Ω = [Ω_1, Ω_1]

    return H, (; x, y, x̄, ȳ, Γ, η, p_f, p, x̄_iv, ȳ_iv, Ω, Q), "rbc"
end


text:
using LightGraphs
using LinearAlgebra
using BALUtils
using Rotations
using PyCall
using StaticArrays

# get camera locations
pairwise_poses = Dict()
es = []
for line in readlines(ARGS[1])
    s = split(line, " ")
    i, j = parse.(Int, s[1:2])
    push!(es, LightGraphs.SimpleEdge(i, j))
    # rotation is row major
    rot = transpose(reshape(parse.(Float64, s[3:11]), 3, 3))
    trans = parse.(Float64, s[12:14])
    pairwise_poses[(i,j)] = (rot, trans)
    pairwise_poses[(j,i)] = (transpose(rot), -(rot' * trans))
end
graph = SimpleGraphFromIterator(es)

# compute global rotations
mst = kruskal_mst(graph) |> SimpleDiGraphFromIterator
_, root = findmax(degree(mst))

poses = Dict()
poses[root] = (Matrix(I, 3, 3), zeros(3))
to_visit = [root]
while length(to_visit) > 0
    v = pop!(to_visit)
    rot_base, trans_base = poses[v]
    for n in neighbors(mst, v)

        rot, trans = pairwise_poses[v, n]
        poses[n] = (rot * rot_base, rot * trans_base + trans)

        push!(to_visit, n)
    end
end

# get point locations and image features
coords = map(x -> Vector(), 1:nv(graph))
camera_intrins = Dict()
open(ARGS[2]) do f
    while !eof(f)
        header = readline(f)
        index = parse(Int, match(r"index = ([0-9.]+)", header).captures[1])
        focal = parse(Float64, match(r"focal = ([0-9.]+)", header).captures[1])
        if haskey(poses, index)
            camera_intrins[index] = (focal, 0, 0)
        end

        px = parse(Float64, match(r"px = ([0-9.]+)", header).captures[1])
        py = parse(Float64, match(r"py = ([0-9.]+)", header).captures[1])

        keys = parse(Int, match(r"keys = ([0-9.]+)", header).captures[1])
        for i in 1:keys
            l = split(readline(f), " ")
            if haskey(poses, index)
                push!(coords[index], (parse(Int, l[1])
                                     , parse(Float64, l[2]) - px/2
                                     , parse(Float64, l[3]) - py/2
                                     ))
            end
        end
    end
end

# mapping from features to camera
features_map = Dict()
for (i, obs) in enumerate(coords)
    for o in obs
        x = o[1]
        y = (i, o[2], o[3])
        if haskey(features_map, x)
            push!(features_map[x], y)
        else
            features_map[x] = [y]
        end
    end
end
features_map = filter(x -> length(x[2]) >= 3, features_map)

renumbering = Dict(zip(1:length(poses), keys(poses)))
renumbering_inv = Dict(zip(keys(poses), 1:length(poses)))
cameras = map(1:length(poses)) do i
    rot, trans = poses[renumbering[i]]
    rv = RodriguesVec(rot)
    intrin = camera_intrins[renumbering[i]]
    Camera(SVector{3,Float64}(rv.sx, rv.sy, rv.sz), SVector{3,Float64}(trans), SVector{3,Float64}(intrin))
end

function intrinsic_matrix(intrin)
    f = intrin[1]
    [f 0 0
     0 f 0
     0 0 1]
end

function rot_matrix(x)
    RotMatrix(RodriguesVec(x[1], x[2], x[3]))
end

# triangulate points
projs = map(c -> intrinsic_matrix(c.intrinsics) * hcat(rot_matrix(c.rotation), c.pose), cameras)
opencv = pyimport("cv2")
points = map(collect(keys(features_map))) do feat
    obs = features_map[feat]
    points = map(x -> reshape([x[2], x[3]], 2, 1), obs)
    proj_mats = map(x -> projs[renumbering_inv[x[1]]], obs)
    points4d = opencv.triangulatePoints(proj_mats[1], proj_mats[2], points[1], points[2])
    points4d[1:3]/points4d[4]
end

points_map = Dict(zip(keys(features_map), 1:length(features_map)))
observations = map(1:length(cameras)) do i
    lobs = []
    for obs in coords[renumbering[i]]
        if haskey(points_map, obs[1])
            push!(lobs, (points_map[obs[1]], obs[2],  obs[3]))
        end
    end
    lobs
end

ba = BA(cameras, observations, points)
writebal(ARGS[3], ba)


text:
# calculates colsums of a large matrix stored in HDFS
# julia -p 3
# include("test_colsum.jl")

using Elly
using HadoopBlocks

const INP = "hdfs://tan@localhost:9000/colsuminp.csv"

if myid() == 1
    function gendata(M, N)
        open(HDFSFile(INP), "w") do f
            for m in 1:M
                write(f, join(map(string, rand(N)), ","), "\n")
            end
        end
    end
    gendata(10^5, 5)
end

# split script here and uncomment below lines to run data generation and computation separately
#using Elly
#using HadoopBlocks
#const INP = "hdfs://tan@localhost:9000/colsuminp.csv"

function findrow(r::HdfsBlockReader, iter_status)
    rec = HadoopBlocks.find_rec(r, iter_status, Vector)
    #HadoopBlocks.logmsg("findrow found rec:$rec")
    rec
end

function maprow(rec)
    #HadoopBlocks.logmsg("maprow called with rec:$rec")
    [tuple([parse(Float64, x) for x in rec]...)]
end

function collectrow(results, rec)
    #HadoopBlocks.logmsg("collectrow called with results:$results rec:$rec")
    isempty(rec) && (return results)
    (results == nothing) && (return rec)
    tuple([results[x]+rec[x] for x in 1:length(results)]...)
end

function reducerow(reduced, results...)
    #HadoopBlocks.logmsg("reducerow called with reduced:$reduced results:$results")
    (nothing == reduced) && (reduced = zeros(Float64, length(results[1])))
    for res in results
        (nothing == res) && continue
        #HadoopBlocks.logmsg("reducerow res:$res")
        #HadoopBlocks.logmsg("reducerow res:$([res...])")
        for x in 1:length(res)
            reduced[x] += res[x]
        end
    end
    #HadoopBlocks.logmsg("reducerow returning reduced:$reduced")
    reduced
end

function wait_results(j_mon)
    loopstatus = true
    while(loopstatus)
        sleep(5)
        jstatus,jstatusinfo = status(j_mon,true)
        ((jstatus == "error") || (jstatus == "complete")) && (loopstatus = false)
        (jstatus == "running") && println("$(j_mon): $(jstatusinfo)% complete...")
    end
    wait(j_mon)
    println("time taken (total time, wait time, run time): $(times(j_mon))")
    println("")
end

j = dmapreduce(MRHdfsFileInput([INP], findrow), maprow, collectrow, reducerow)
wait_results(j)
println(results(j))

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
######################
# 1: The Julia type for ToricDivisors
######################

struct ToricDivisor
           polymake_divisor::Polymake.BigObject
end
export ToricDivisor

function pm_tdivisor(td::ToricDivisor)
    return td.polymake_divisor
end

######################
# 2: Generic constructors
######################

@doc Markdown.doc"""
    ToricDivisor(coeffs::AbstractVector, v::AbstractNormalToricVariety )

Construct the torus invariant divisor on the normal toric variety `v` as linear combination of the torus invariant prime divisors of `v`. The coefficients of thi linear combination are passed as list of integers as first argument.

# Examples
```jldoctest
julia> show( ToricDivisor( [1,1,2], projective_space( 2 ) ) )
A torus invariant divisor on a normal toric variety
```
"""
function ToricDivisor( coeffs::AbstractVector, v::AbstractNormalToricVariety )
    if length(coeffs) != pm_ntv(v).N_RAYS
        throw(ArgumentError("Number of coefficients needs to match number of prime divisors!"))
    end
    ptd = Polymake.fulton.TDivisor(COEFFICIENTS=coeffs)
    pmntv = pm_ntv(v)
    Polymake.add(pmntv, "DIVISOR", ptd)
    return ToricDivisor(ptd)
end
export ToricDivisor

######################
# 3: Properties
######################

@doc Markdown.doc"""
    is_cartier( d::ToricDivisor )

Checks if the divisor `d` is Cartier.

# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> iscartier(td)
true
```
"""
iscartier(td::ToricDivisor) = pm_tdivisor(td).CARTIER::Bool
export iscartier

@doc Markdown.doc"""
    isprincipal(td::ToricDivisor)

Determine whether the toric divisor `td` is principal.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isprincipal(td)
false
```
"""
isprincipal(td::ToricDivisor) = pm_tdivisor(td).PRINCIPAL::Bool
export isprincipal

@doc Markdown.doc"""
    isbasepoint_free(td::ToricDivisor)

Determine whether the toric divisor `td` is basepoint free.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isbasepoint_free(td)
true
```
"""
isbasepoint_free(td::ToricDivisor) = pm_tdivisor(td).BASEPOINT_FREE::Bool
export isbasepoint_free

@doc Markdown.doc"""
    iseffective(td::ToricDivisor)

Determine whether the toric divisor `td` is effective.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> iseffective(td)
true
```
"""
iseffective(td::ToricDivisor) = pm_tdivisor(td).EFFECTIVE::Bool
export iseffective

@doc Markdown.doc"""
    isintegral(td::ToricDivisor)

Determine whether the toric divisor `td` is integral.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isintegral(td)
true
```
"""
isintegral(td::ToricDivisor) = pm_tdivisor(td).INTEGRAL::Bool
export isintegral

@doc Markdown.doc"""
    isample(td::ToricDivisor)

Determine whether the toric divisor `td` is ample.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isample(td)
false
```
"""
isample(td::ToricDivisor) = pm_tdivisor(td).AMPLE::Bool
export isample

@doc Markdown.doc"""
    isvery_ample(td::ToricDivisor)

Determine whether the toric divisor `td` is very ample.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isvery_ample(td)
false
```
"""
isvery_ample(td::ToricDivisor) = pm_tdivisor(td).VERY_AMPLE::Bool
export isvery_ample

@doc Markdown.doc"""
    isnef(td::ToricDivisor)

Determine whether the toric divisor `td` is nef.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isnef(td)
true
```
"""
isnef(td::ToricDivisor) = pm_tdivisor(td).NEF::Bool
export isnef

@doc Markdown.doc"""
    isq_cartier(td::ToricDivisor)

Determine whether the toric divisor `td` is Q-Cartier.
# Examples
```jldoctest
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isq_cartier(td)
true
```
"""
isq_cartier(td::ToricDivisor) = pm_tdivisor(td).Q_CARTIER::Bool
export isq_cartier

@doc Markdown.doc"""
    polyhedron_of_divisor(td::ToricDivisor)

Construct the polyhedron $P_D$ of a torus invariant divisor $D:=td$ as in 4.3.2
of [CLS11](@cite). The lattice points of this polyhedron correspond to the
global sections of the divisor.

# Examples
The polyhedron of the divisor with all coefficients equal to zero is a point,
if the ambient variety is complete. Changing the coefficients corresponds to
moving hyperplanes. One direction moves the hyperplane away from the origin,
the other moves it across. In the latter case there are no global sections
anymore and the polyhedron becomes empty.
```
julia> H = hirzebruch_surface(4)
A normal toric variety corresponding to a polyhedral fan in ambient dimension 2

julia> td0 = ToricDivisor([0,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isfeasible(polyhedron_of_divisor(td0))
true

julia> dim(polyhedron_of_divisor(td0))
0

julia> td1 = ToricDivisor([1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isfeasible(polyhedron_of_divisor(td1))
true

julia> td2 = ToricDivisor([-1,0,0,0], H)
A torus invariant divisor on a normal toric variety

julia> isfeasible(polyhedron_of_divisor(td2))
false
```
"""
function polyhedron_of_divisor(td::ToricDivisor)
    pmtd = pm_tdivisor(td)
    return Polyhedron(pmtd.SECTION_POLYTOPE)
end

###############################################################################
###############################################################################
### Display
###############################################################################
###############################################################################
function Base.show(io::IO, td::ToricDivisor)
    print(io, "A torus invariant divisor on a normal toric variety")
end


text:
"""
    _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation

Parse `SBML.GeneProductAssociation` structure to the simpler GeneAssociation.
The input must be (implicitly) in a positive DNF.
"""
function _parse_grr(gpa::SBML.GeneProductAssociation)::GeneAssociation
    parse_ref(x) =
        typeof(x) == SBML.GPARef ? [x.gene_product] :
        begin
            @_models_log @warn "Could not parse a part of gene association, ignoring: $x"
            String[]
        end
    parse_and(x) =
        typeof(x) == SBML.GPAAnd ? vcat([parse_and(i) for i in x.terms]...) : parse_ref(x)
    parse_or(x) =
        typeof(x) == SBML.GPAOr ? vcat([parse_or(i) for i in x.terms]...) : [parse_and(x)]
    return parse_or(gpa)
end

"""
    _unparse_grr(
        ::Type{SBML.GeneProductAssociation},
        x::GeneAssociation,
    )::SBML.GeneAssociation

Convert a GeneAssociation to the corresponding `SBML.jl` structure.
"""
function _unparse_grr(
    ::Type{SBML.GeneProductAssociation},
    x::GeneAssociation,
)::SBML.GeneAssociation
    SBML.GPAOr([SBML.GPAAnd([SBML.GPARef(j) for j in i]) for i in x])
end

"""
    _parse_grr(s::String)::GeneAssociation

Parse a DNF gene association rule in format `(YIL010W and YLR043C) or (YIL010W
and YGR209C)` to `GeneAssociation. Also accepts `OR`, `|`, `||`, `AND`, `&`,
and `&&`.

# Example
```
julia> _parse_grr("(YIL010W and YLR043C) or (YIL010W and YGR209C)")
2-element Array{Array{String,1},1}:
 ["YIL010W", "YLR043C"]
 ["YIL010W", "YGR209C"]
```
"""
function _parse_grr(s::String)::GeneAssociation
    # first get the gene id list in string format
    gene_reaction_rules = Vector{Vector{String}}()
    or_genes = split(s, r"\s?(or|OR|(\|\|)|\|)\s?") # separate or terms
    for or_gene in or_genes
        and_genes = split(replace(or_gene, r"\(|\)" => ""), r"\s?(and|AND|(\&\&)|\&)\s?")
        push!(gene_reaction_rules, and_genes)
    end
    return gene_reaction_rules
end

"""
    unparse_grr(grr::Vector{Vector{Gene}}

Converts a nested string gene reaction array  back into a gene reaction rule
string.

# Example
```
julia> _unparse_grr(String, [["YIL010W", "YLR043C"], ["YIL010W", "YGR209C"]])
"(YIL010W and YLR043C) or (YIL010W and YGR209C)"
```
"""
function _unparse_grr(::Type{String}, grr::GeneAssociation)::String
    grr_strings = String[]
    for gr in grr
        push!(grr_strings, "(" * join([g for g in gr], " and ") * ")")
    end
    grr_string = join(grr_strings, " or ")
    return grr_string
end


text:
# Dual infeasible test problems
# Choose A, P and q in a way such that the problem becomes unbounded below
# here the last element of x appears in the cost function with a negative sign and is unconstrained

rng = Random.MersenneTwister(555)
nn = 1

@testset "Dual infeasible QP problems - Testset 1" begin
  for iii = 1:nn

  # choose size of problem
  n = rand(rng,5:50)
  m = 2*n
  A = sprandn(m,n,0.7)*50
  A[:,end] .= 0

  P = spzeros(n,n)

  q = randn(n,1)*50
  q[end] = -1
  q = vec(q)
  xtrue = randn(n,1)*50
  strue = rand(m,1)*50
  b = A*xtrue+strue
  b = vec(b)

  constraint = COSMO.Constraint(-A,b,COSMO.Nonnegatives)

  settings = COSMO.Settings(max_iter=10000,eps_abs = 1e-5,eps_rel=1e-5)
  model = COSMO.Model()
  assemble!(model,P,q,[constraint], settings = settings)
  res = COSMO.optimize!(model);

  @test res.status == :Dual_infeasible
  end
end


text:
# bind console api for julia

# for linux
const STDIN_FILENO = 0
const ICANON = 0x02
const ECHO = 0x08
const TCSANOW = 0

# for windows
const STD_OUTPUT_HANDLE = reinterpret(Culong, Clong(-11))
const HANDLE = Ptr{Cvoid}

mutable struct termios
    c_iflag::Cuint		#  input mode flags
    c_oflag::Cuint		#  output mode flags
    c_cflag::Cuint		#  control mode flags
    c_lflag::Cuint		#  local mode flags
    c_line::Cuchar		#  line discipline
    c_cc::NTuple{32,Cuchar}	#  control characters
    c_ispeed::Cint		#  input speed
    c_ospeed::Cint		#  output speed
end

mutable struct COORD
    x::Cshort
    y::Cshort
end

termios() = termios(0, 0, 0, 0, 0, ntuple(_ -> Cuchar(0), 32), 0, 0)

function getch_linux()::Cint
    oldt = termios()
    ccall(:tcgetattr, Cint, (Cint, Ref{termios}), STDIN_FILENO, Ref(oldt))
    newt = oldt
    newt.c_lflag &= ~(ICANON | ECHO)
    ccall(:tcsetattr, Cint, (Cint, Cint, Ref{termios}), STDIN_FILENO, TCSANOW, Ref(newt))
    ch = ccall(:getchar, Cint, ())
    ccall(:tcsetattr, Cint, (Cint, Cint, Ref{termios}), STDIN_FILENO, TCSANOW, Ref(oldt))
    return ch
end

getch()::Cint = begin
    @static if Sys.iswindows()
        return ccall(:_getch, Cint, ())
    else
        return getch_linux()
    end
end

function clean_console()
    @static if Sys.iswindows()
        ccall(:system, Cint, (Cstring,), "cls")
    else
        print("\033[2J")
    end
end

function set_cursor_position(x::Int, y::Int)
    @static if Sys.iswindows()
        handle = ccall(:GetStdHandle, HANDLE, (UInt32,), STD_OUTPUT_HANDLE)
        pos = COORD(x % Cshort, y % Cshort)
        ccall(:SetConsoleCursorPosition, Cint, (HANDLE, COORD), handle, pos)
    else
        @printf("\033[%d;%dH", y + 1, x + 1)
    end
end


text:
include("../utility.jl")

path = Pkg.dir("MKL") * "/benchmark/b1/"

df = DataFrames.readtable(path*"b1_data.csv")
N = union(df[:N])
T = union(df[:Type])
imp = union(df[:Implementation])
gen = union(df[:Generic])

# this is an option to get normalized data relatively to MKL
# for i in 1:nrow(df)
#     imkl = 3*div(i-1,3) + 2 # MKL row
#     df[i,:CPU_mean] = df[i,:CPU_mean] / df[imkl,:CPU_mean]
#     println(df[i,:CPU_mean])
# end

for n in N
    for t in T

        println(n)
        println(t)

        filter_set = (df[:Type] .== string(t)) & (df[:N] .== n)
        sdf = df[filter_set , :]
        println(sdf)

        set_default_plot_size(20cm, 10cm)
        p = Gadfly.plot(
                    sdf, xgroup="Generic", x="Implementation", y="CPU_mean",
                    color="Implementation",
                    Scale.color_discrete_manual(Tdarkgray,Tblue,Tlightgray),
                    Scale.y_continuous(minvalue=0, maxvalue=20),
                    Geom.subplot_grid(Geom.bar, Guide.xticks(orientation=:vertical, label=false)),
                    Guide.xlabel(""),
                    Guide.ylabel(""),
                    Guide.title("CPU time on n-dim vectors (n=1e"*string(convert(Int32,log10(n)))),
                    Theme(bar_spacing=0.1cm)
                    )

        draw(PNG(path*"output/b1_$t"*"_$n.png", 25cm, 12cm), p)
        # draw(SVG(path*"output/b1_$t"*"_$n.svg", 25cm, 12cm), p)
        p
    end
end

for n in N

    println(n)

    filter_set = (df[:N] .== n)
    sdf = df[filter_set , :]
    println(sdf)

    p = Gadfly.plot(
                layer(
                sdf, xgroup="Generic", ygroup="Type", x="Implementation", y="CPU_mean",
                color="Implementation",
                Geom.subplot_grid(Geom.bar, Guide.xticks(orientation=:vertical, label=false))
                ),
                Scale.color_discrete_manual(Tdarkgray,Tblue,Tlightgray),
                Scale.y_continuous(minvalue=0, maxvalue=20),
                Guide.xlabel(""),
                Guide.ylabel(""),
                Guide.title("CPU time on n-dim vectors (n=1e"*string(convert(Int32,log10(n)))*")"),
                Theme(bar_spacing=0.1cm)
                )

    # draw(PNG(path*"output/b1_$n.png", 25cm, 12cm), p)
    draw(SVG(path*"output/b1_$n.svg", 29cm, 20cm), p)
    draw(PDF(path*"output/b1_$n.pdf", 29cm, 20cm), p)
    p
end

# POST TREATMENT FOR .CSV EXPORT (using DatFrames)
# begin
#     hrow = union(df[:Generic]) # list of generic functions
#     hcol = ("BLAS","MKL","VJulia") #union(df[:Implementation])
#     resdf = DataFrame()
#
#     col = ASCIIString[]
#     for i=1:length(hcol)
#         push!(col,string(hcol[i]," (mean)"))
#         push!(col,string(hcol[i]," (conf)"))
#     end
#     resdf[:Generic] = col
#
#     for imp in hrow
#         resdf[symbol(imp)] = [0.0 for i=1:2*length(hcol)]
#     end
#
#     subdf = df[ df[:Type] .== "Float32", :]
#
#     for i=1:nrow(subdf)
#         icol = 1 + findfirst(hrow,subdf[i,1])
#         irow = 1 + 2*(findfirst(hcol,subdf[i,2])-1)
#         resdf[irow,icol] = subdf[i,5]
#         resdf[irow+1,icol] = subdf[i,6]
#     end
#     # println(hcol)
#     # println(hrow)
#     println(df)
#     println(subdf)
#
#     println(resdf)
#     # writetable("output.csv", resdf, separator='\t')
#     writetable(path * "b1_cpu_$n.csv", resdf, separator='\t')
#     # writetable(path * "output.csv", df, separator='\t')
# end

# Float32 / Float64 PERFORMANCE COMPARISON
# begin
#     hrow = union(df[:Generic]) # list of generic functions
#     hcol = ("BLAS","MKL","VJulia") #union(df[:Implementation])
#     resdf = DataFrame()
#
#     col = ASCIIString[]
#     for i=1:length(hcol)
#         push!(col,string(hcol[i],""))
#     end
#     resdf[:Generic] = col
#
#     for imp in hrow
#         resdf[symbol(imp)] = [0.0 for i=1:1*length(hcol)]
#     end
#
#     subdf_32 = df[ df[:Type] .== "Float32", :]
#     subdf_64 = df[ df[:Type] .== "Float64", :]
#
#     for i=1:nrow(subdf_64)
#         subdf_64[i,5] = subdf_64[i,5] / subdf_32[i,5]
#         println(subdf_64[i,5])
#     end
#
#     for i=1:nrow(subdf_64)
#         icol = 1 + findfirst(hrow,subdf_64[i,1])
#         irow = 1 + 1*(findfirst(hcol,subdf_64[i,2])-1)
#         resdf[irow,icol] = subdf_64[i,5]
#     end
#
#     # println(subdf_64)
#     # println(resdf)
#     writetable(path * "b1_Float32_vs_Float64_$n.csv", resdf, separator='\t')
#
#     # println(resdf)
#     # writetable("output.csv", resdf, separator='\t')
#
#     # writetable(path * "output.csv", df, separator='\t')
# end

# set_default_plot_size(10cm, 10cm)
# p = plot(df, ygroup="Type", x="Method", y="CPU_mean", color="Implementation", Geom.subplot_grid(Geom.bar, Guide.xticks(orientation=:vertical)),Theme(bar_spacing=0.1cm))


text:
using BenchmarkTools
using Yao
using Yao.Intrinsics
using Yao.Boost

v = randn(ComplexF64, 1 << 16)
bench = BenchmarkGroup()
xg = xgate(16, 2)
#bg["indices_with"] = @benchmarkable indices_with(2, bss)
bg = bench["Basic Gate"] = BenchmarkGroup()
bg["X"] = @benchmarkable xgate(16, 2)
bg["Y"] = @benchmarkable ygate(16, 2)
bg["Z"] = @benchmarkable zgate(16, 2)
bg["X*v"] = @benchmarkable $xg * v

############# controlled gates ##################
bg = bench["CZ Gate"] = BenchmarkGroup()
bg["CZ"] = @benchmarkable czgate(ComplexF64, 16, 3, 7)
bg["Diag-CZ"] = @benchmarkable controlled_U1(16, PAULI_Z, [3], 7)
bg["General-CZ"] = @benchmarkable general_controlled_gates(16, [P1], [7], [PAULI_Z], [3])

bg = bench["CX Gate"] = BenchmarkGroup()
bg["CNOT"] = @benchmarkable cxgate(ComplexF64, 16, 7, 3)
bg["PM-CNOT"] = @benchmarkable controlled_U1(16, PAULI_X, [3], 7)
bg["General-CNOT"] = @benchmarkable general_controlled_gates(16, [P1], [7], [PAULI_X], [3])

bg = bench["CY Gate"] = BenchmarkGroup()
bg["CY"] = @benchmarkable cygate(ComplexF64, 16, 7, 3)
bg["PM-CY"] = @benchmarkable controlled_U1(16, PAULI_Y, [3], 7)
bg["General-CY"] = @benchmarkable general_controlled_gates(16, [P1], [7], [PAULI_Y], [3])

##### TOFFOLI gate
bg = bench["Toffoli Gate"] = BenchmarkGroup()
bg["Toffoli"] = @benchmarkable toffoligate(16, 2, 3, 1)
bg["General-Toffoli"] = @benchmarkable general_controlled_gates(
    16,
    [P1, P1],
    [2, 3],
    [PAULI_X],
    [1],
)

showall(run(bench, verbose = true))


text:
module FlightGNC
using LinearAlgebra

using Reexport
@reexport using FSimBase
const FSBase = FSimBase
import FSimBase: State, Params, Dynamics!, Dynamics
import FSimBase: Command
using UnPack, ComponentArrays
using Plots  # required for fig_print

# exports -------------------------------------------------------------
## algorithms
export BPNG, BPNG_cmd
export Bias_zero, Bias_IACG_StationaryTarget, Bias_IACG_StationaryTarget_2D

## environments
export PointMassMissile, PursuerEvadorMissile

## utils
export fig_print

# includes -------------------------------------------------------------
## algorithms
include("algorithms/BPNG.jl")

## environments
include("environments/missiles.jl")

## utils
include("utils/fig_print.jl")
end


text:
module LanguageServer
import URIParser
using JSON, REPL, CSTParser, DocumentFormat, SymbolServer, StaticLint
using CSTParser: EXPR, Tokenize.Tokens, typof, kindof, parentof, valof
using StaticLint: refof, scopeof, bindingof
using UUIDs
using Base.Docs, Markdown
import JSONRPC
using JSONRPC: Outbound, @dict_readable

export LanguageServerInstance, runserver

include("exception_types.jl")
include("uri2.jl")
include("protocol/protocol.jl")
include("extensions/extensions.jl")
include("document.jl")
include("languageserverinstance.jl")
include("runserver.jl")
include("staticlint.jl")

include("requests/misc.jl")
include("requests/textdocument.jl")
include("requests/features.jl")
include("requests/hover.jl")
include("requests/completions.jl")
include("requests/workspace.jl")
include("requests/actions.jl")
include("requests/init.jl")
include("utilities.jl")

end


text:
export AutoBcast

"""
    AutoBcast{T,N} <: IWrapper{T}

A vectorized variable.
"""
struct AutoBcast{T,N} <: IWrapper{T} x::Vector{T} end
AutoBcast(x::Vector{T}) where {T} = AutoBcast{T, length(x)}(x)
AutoBcast(x::AutoBcast{T,N}) where {T,N} = x # to avoid ambiguity error
AutoBcast{T,N}(x::AutoBcast{T,N}) where {T,N} = x
value(x::AutoBcast) = x.x
NiLangCore.chfield(x::AutoBcast, ::typeof(value), xval) = chfield(x, Val(:x), xval)
Base.zero(x::AutoBcast) = AutoBcast(zero(x.x))
Base.zero(::Type{AutoBcast{T,N}}) where {T,N} = AutoBcast{T,N}(zeros(T, N))
Base.length(ab::AutoBcast{T,N}) where {T, N} = N

for F1 in [:(Base.:-), :INC, :FLIP, :DEC]
    @eval function $F1(a!::AutoBcast)
        @instr @invcheckoff @inbounds for i=1:length(a!)
            $F1(a!.x[i])
        end
        a!
    end
end

for F2 in [:SWAP, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    F2 != :SWAP && @eval function $F2(a::AutoBcast, b::Real)
        @instr @invcheckoff @inbounds for i=1:length(a)
            $F2(a.x[i], b)
        end
        a, b
    end
    @eval function $F2(a::AutoBcast, b::AutoBcast)
        @instr @invcheckoff @inbounds for i=1:length(a)
            $F2(a.x[i], b.x[i])
        end
        a, b
    end
end

for F3 in [:ROT, :IROT, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    if !(F3 in [:ROT, :IROT])
        @eval function $F3(a::AutoBcast, b::Real, c::Real)
            @instr @invcheckoff @inbounds for i=1:length(a)
                $F3(a.x[i], b, c)
            end
            a, b, c
        end
        @eval function $F3(a::AutoBcast, b::Real, c::AutoBcast)
            @instr @invcheckoff for i=1:length(a)
                $F3(a.x[i], b, c.x[i])
            end
            a, b, c
        end
    end
    @eval function $F3(a::AutoBcast, b::AutoBcast, c::Real)
        @instr @invcheckoff @inbounds for i=1:length(a)
            $F3(a.x[i], b.x[i], c)
        end
        a, b, c
    end
    @eval function $F3(a::AutoBcast, b::AutoBcast, c::AutoBcast)
        @instr @invcheckoff @inbounds for i=1:length(a)
            $F3(a.x[i], b.x[i], c.x[i])
        end
        a, b, c
    end
end

(f::PlusEq{typeof(identity)})(x::T, y::T) where T<:AutoBcast = invoke(f, Tuple{T,T} where T, x, y)
(f::MinusEq{typeof(identity)})(x::T, y::T) where T<:AutoBcast = invoke(f, Tuple{T,T} where T, x, y)


text:
module ClickHouse

include("Net.jl")
include("Query.jl")

export ClickHouseSock
export Block
export select
export select_callback
export select_channel
export select_df
export insert
export execute
export connect
export ping
export ClickHouseServerException

end # module


text:
using Random, Test, HDF5

@testset "hyperslab" begin
    N = 10
    v = [randstring(rand(5:10)) for i in 1:N, j in 1:N]

    fn = tempname()
    h5open(fn, "w") do f
        f["data"] = v
    end

    h5open(fn, "r") do f
        dset = f["data"]
        indices = (1, 1)
        @test dset[indices...] == v[indices...]

        indices = (1:10, 1)
        @test dset[indices...] == v[indices...]

        indices = (1, 1:10)
        @test dset[indices...] == v[indices...]

        indices = (1:2:10, 1:3:10)
        @test dset[indices...] == v[indices...]
    end
end


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

JLLWrappers.@generate_main_file_header("unixODBC")
JLLWrappers.@generate_main_file("unixODBC", UUID("1841a5aa-d9e2-579c-8226-32ed2af93ab1"))
end  # module unixODBC_jll


text:
# Functions for time-domain modeling and inversion using devito
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January, 2017
# Updated, December 2020, Mathias Louboutin, mlouboutin3@gatech.edu

#############################################################################
# Containers
include("Types/ModelStructure.jl")    # model container
include("Types/InfoStructure.jl") # basic information required by all operators
include("Types/GeometryStructure.jl") # source or receiver setup, recording time and sampling
include("Types/OptionsStructure.jl")

#############################################################################
# Abstract vectors
include("Types/judiWavefield.jl") # dense RHS (wavefield)
include("Types/judiWeights.jl")    # Extended source weight vector
include("Types/judiVector.jl")    # Julia data container
include("Types/judiComposites.jl")    # A composite type to work with hcat/vcat of judi types
include("Types/typeutils.jl") # dense RHS (wavefield)

#############################################################################
# Utils
include("Utils/auxiliaryFunctions.jl")

#############################################################################
# PDE solvers
include("Modeling/python_interface.jl")  # forward/adjoint linear/nonlinear modeling
include("Modeling/time_modeling_serial.jl")  # forward/adjoint linear/nonlinear modeling
include("Modeling/extended_source_interface_serial.jl")  # forward/adjoint linear/nonlinear modeling w/ extended source
include("Modeling/fwi_objective_serial.jl")  # FWI objective function value and gradient
include("Modeling/lsrtm_objective_serial.jl")  # LSRTM objective function value and gradient
include("Modeling/twri_objective_serial.jl")  # TWRI objective function value and gradient
include("Modeling/distributed.jl") # Modeling functions utilities

#############################################################################
# Linear operators
include("LinearOperators/judiRHS.jl")   # sparse RHS (point source(s))
include("LinearOperators/judiExtendedSource.jl")   # sparse RHS (point source(s))
include("LinearOperators/judiModeling.jl")  # nonlinear modeling operator F (no projection operators)
include("LinearOperators/judiProjection.jl")    # source/receiver projection operator
include("LinearOperators/judiLRWF.jl")   # low rank wavefield (point source(s))
include("LinearOperators/judiPDEfull.jl")   # modeling operator with source and receiver projection: P*F*P'
include("LinearOperators/judiPDEextended.jl")   # modeling operator for extended sources
include("LinearOperators/judiPDE.jl")   # modeling operator with lhs projection only: P*F
include("LinearOperators/judiJacobian.jl")  # linearized modeling operator J
include("LinearOperators/judiJacobianExtendedSource.jl")  # Jacobian of extended source modeling operator

#############################################################################
# Preconditioners and optimization
include("Utils/seismic_preconditioners.jl")

#############################################################################
# Utility types
const SourceTypes = Union{judiVector, Tuple{judiWeights, judiLRWF}}
# PDE types are callable w.r.t non-linear parameters, i.e F(model) or F(;m=m, epsilon=new_epsilon)
const pde_types = Union{judiModeling, judiModelingAdjoint, judiPDEfull, judiPDE, judiPDEadjoint,
                        judiJacobian, judiJacobianExQ, judiPDEextended}

if VERSION>v"1.2"
  function (F::pde_types)(m::Model)
    Fl = deepcopy(F)
    Fl.model.n = m.n
    Fl.model.d = m.d
    Fl.model.o = m.o
    for (k, v) in m.params
      Fl.model.params[k] = v
    end
    Fl
  end

  function (F::pde_types)(;kwargs...)
    Fl = deepcopy(F)
    for (k, v) in kwargs
      k in keys(Fl.model.params) && Fl.model.params[k] .= v
    end
    Fl
  end
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
@inline function pow_by_square(_v, e::IntegerTypesHW)
  v = float(_v)
  if e < 0
    v = y = inv(v)
    e = -e
  elseif e == 0
    return one(v)
  else
    y = v
  end
  tz = trailing_zeros(e)
  e >>= (tz + one(tz))
  while tz > zero(tz)
    y = Base.FastMath.mul_fast(y, y)
    tz -= one(tz)
  end
  x = y
  while e ≠ zero(e)
    y = Base.FastMath.mul_fast(y, y)
    tz = trailing_zeros(e)
    e >>= (tz + one(tz))
    while tz > zero(tz)
      y = Base.FastMath.mul_fast(y, y)
      tz -= one(tz)
    end
    x = Base.FastMath.mul_fast(x, y)
  end
  return x
end
@generated function pow_by_square(_v, ::StaticInt{E}) where {E}
  e = E
  q = Expr(:block, Expr(:meta, :inline), :(v = float(_v)))
  xdefined = false
  if e < 0
    push!(q.args, :(v = y = inv(v)))
    e = -e
  else
    push!(q.args, :(y = v))
  end
  mf = Base.FastMath.mul_fast
  while e ≠ zero(e)
    xdefined && push!(q.args, :(y = $mf(y, y)))
    tz = trailing_zeros(e)
    e >>= (tz + one(tz))
    while tz > zero(tz)
      push!(q.args, :(y = $mf(y, y)))
      tz -= one(tz)
    end
    if xdefined
      push!(q.args, :(x = $mf(x, y)))
    else
      xdefined = true
      push!(q.args, :(x = y))
    end
  end
  push!(q.args, :x)
  return q
end
# 5 = 101 = 2^2 + 2^0 # x^4 * x^1
# x^5 = x^4 * x

@inline Base.:^(v::AbstractSIMD{W,T}, i::IntegerTypesHW) where {W,T} = pow_by_square(v, i)
@inline Base.:^(
  v::AbstractSIMD{W,T},
  i::IntegerTypesHW,
) where {W,T<:Union{Float32,Float64}} = pow_by_square(v, i)
@inline Base.:^(v::AbstractSIMD, ::StaticInt{N}) where {N} =
  pow_by_square(v, StaticInt{N}())
@inline Base.FastMath.pow_fast(v::AbstractSIMD, ::StaticInt{N}) where {N} =
  pow_by_square(v, StaticInt{N}())
@inline Base.FastMath.pow_fast(v::AbstractSIMD{W,T}, i::IntegerTypesHW) where {W,T} =
  pow_by_square(v, i)
@inline Base.FastMath.pow_fast(
  v::AbstractSIMD{W,T},
  i::IntegerTypesHW,
) where {W,T<:Union{Float32,Float64}} = pow_by_square(v, i)
@inline Base.FastMath.pow_fast(v::AbstractSIMD, x::FloatingTypes) =
  exp2(Base.FastMath.log2_fast(v) * x)
@inline Base.FastMath.pow_fast(v::FloatingTypes, x::AbstractSIMD) =
  exp2(Base.FastMath.log2_fast(v) * x)
@inline Base.FastMath.pow_fast(v::AbstractSIMD, x::AbstractSIMD) =
  exp2(Base.FastMath.log2_fast(v) * x)
@inline Base.literal_pow(::typeof(^), x::AbstractSIMD, ::Val{N}) where {N} =
  pow_by_square(x, StaticInt(N))
# @inline relu(x) = (y = zero(x); ifelse(x > y, x, y))
@inline relu(x) = (y = zero(x); ifelse(x < y, y, x))

Base.sign(v::AbstractSIMD) = ifelse(v > 0, one(v), -one(v))

@inline Base.fld(x::AbstractSIMD, y::AbstractSIMD) = div(promote_div(x, y)..., RoundDown)
@inline Base.fld(x::AbstractSIMD, y::Real) = div(promote_div(x, y)..., RoundDown)
@inline Base.fld(x::Real, y::AbstractSIMD) = div(promote_div(x, y)..., RoundDown)

@inline function Base.div(
  x::AbstractSIMD{W,T},
  y::AbstractSIMD{W,T},
  ::RoundingMode{:Down},
) where {W,T<:IntegerTypes}
  d = div(x, y)
  d - (signbit(x ⊻ y) & (d * y != x))
end

@inline Base.mod(x::AbstractSIMD{W,T}, y::AbstractSIMD{W,T}) where {W,T<:IntegerTypes} =
  ifelse(y == -1, zero(x), x - fld(x, y) * y)

@inline Base.mod(x::AbstractSIMD{W,T}, y::AbstractSIMD{W,T}) where {W,T<:Unsigned} =
  rem(x, y)

@inline function Base.mod(
  x::AbstractSIMD{W,T1},
  y::T2,
) where {W,T1<:SignedHW,T2<:UnsignedHW}
  _x, _y = promote_div(x, y)
  unsigned(mod(_x, _y))
end
@inline function Base.mod(
  x::AbstractSIMD{W,T1},
  y::T2,
) where {W,T1<:UnsignedHW,T2<:SignedHW}
  _x, _y = promote_div(x, y)
  signed(mod(_x, _y))
end
@inline function Base.mod(
  x::AbstractSIMD{W,T1},
  y::AbstractSIMD{W,T2},
) where {W,T1<:SignedHW,T2<:UnsignedHW}
  _x, _y = promote_div(x, y)
  unsigned(mod(_x, _y))
end
@inline function Base.mod(
  x::AbstractSIMD{W,T1},
  y::AbstractSIMD{W,T2},
) where {W,T1<:UnsignedHW,T2<:SignedHW}
  _x, _y = promote_div(x, y)
  signed(mod(_x, _y))
end

@inline Base.mod(i::AbstractSIMD{<:Any,<:IntegerTypes}, r::AbstractUnitRange{<:IntegerTypes}) =
  mod(i - first(r), length(r)) + first(r)

@inline Base.mod(x::AbstractSIMD, y::NativeTypes) = mod(promote_div(x, y)...)
@inline Base.mod(x::NativeTypes, y::AbstractSIMD) = mod(promote_div(x, y)...)

# avoid ambiguity with clamp(::Missing, lo, hi) in Base.Math at math.jl:1258
# but who knows what would happen if you called it
for (X, L, H) in Iterators.product(fill([:Any, :Missing, :AbstractSIMD], 3)...)
  any(==(:AbstractSIMD), (X, L, H)) || continue
  @eval @inline function Base.clamp(x::$X, lo::$L, hi::$H)
    x_, lo_, hi_ = promote(x, lo, hi)
    ifelse(x_ > hi_, hi_, ifelse(x_ < lo_, lo_, x_))
  end
end

@inline Base.FastMath.hypot_fast(x::AbstractSIMD, y::AbstractSIMD) =
  sqrt(Base.FastMath.add_fast(Base.FastMath.mul_fast(x, x), Base.FastMath.mul_fast(y, y)))

@inline Base.clamp(x::AbstractSIMD{<:Any,<:IntegerTypes}, r::AbstractUnitRange{<:IntegerTypes}) =
  clamp(x, first(r), last(r))

@inline function Base.gcd(
  a::AbstractSIMDVector{W,I},
  b::AbstractSIMDVector{W,I},
) where {W,I<:Base.HWReal}
  aiszero = a == zero(a)
  biszero = b == zero(b)
  absa = abs(a)
  absb = abs(b)
  za = trailing_zeros(a)
  zb = ifelse(biszero, zero(b), trailing_zeros(b))
  k = min(za, zb)
  u = unsigned(ifelse(biszero, zero(a), abs(a >> za)))
  v = unsigned(ifelse(aiszero, zero(b), abs(b >> zb)))
  ne = u ≠ v
  while vany(ne)
    ulev = (u > v) & ne
    t = u
    u = ifelse(ulev, v, u)
    v = ifelse(ulev, t, v)
    d = v - u
    v = ifelse(ne, d >> trailing_zeros(d), v)
    ne = u ≠ v
  end
  ifelse(aiszero, absb, ifelse(biszero, absa, (u << k) % I))
end
@inline Base.gcd(a::VecUnroll, b::Real) = VecUnroll(fmap(gcd, data(a), b))
@inline Base.gcd(a::Real, b::VecUnroll) = VecUnroll(fmap(gcd, a, data(b)))
@inline Base.gcd(a::VecUnroll, b::VecUnroll) = VecUnroll(fmap(gcd, data(a), data(b)))
@inline function Base.lcm(a::AbstractSIMD, b::AbstractSIMD)
  z = zero(a)
  isz = (a == z) | (b == z)
  ifelse(isz, z, (b ÷ gcd(b, a)) * a)
end
@inline Base.lcm(a::AbstractSIMD, b::Real) = ((c, d) = promote(a, b); lcm(c, d))
@inline Base.lcm(a::Real, b::AbstractSIMD) = ((c, d) = promote(a, b); lcm(c, d))


text:
@test modeorthonorm("Hermite",100,1e-10) #test orthonormality
@test modeorthonorm("Hermite",100,1e-10,3.) #test scaling with frequency
@test modeorthonorm("Hermite",370,1e-8) #test maximum cutoff