text:
module StatsAPI

include("statisticalmodel.jl")
include("regressionmodel.jl")

"""
    params(model)

Return all parameters of a model.
"""
function params end

function params! end

#    pairwise(f, x[, y])
#
# Return a matrix holding the result of applying `f` to all possible pairs
# of entries in iterators `x` and `y`, and return it.
#
# This generic function is owned by StatsBase.jl, which is the sole provider
# of the default definition.
function pairwise end

#    pairwise!(f, dest::AbstractMatrix, x[, y])
#
# Store in matrix `dest` the result of applying `f` to all possible pairs
# of entries in iterators `x` and `y`, and return it.
#
# This generic function is owned by StatsBase.jl, which is the sole provider
# of the default definition.
function pairwise! end

end # module

max_stars_repo_name: JuliaStats/StatsAPI

max_issues_repo_name: JuliaStats/StatsAPI

max_forks_repo_name: JuliaStats/StatsAPI


text:
using ConvexTests, COSMO

@info "Starting COSMO tests"

do_tests("COSMO", () -> COSMO.Optimizer(decompose = false, max_iter = 40000, verbose = false, eps_abs = 5e-7, eps_rel = 5e-7); exclude = [r"mip", r"dual"], description = "Tests run with `decompose = false`, `max_iter = 40000`, `eps_abs = 5e-7`, and `eps_rel = 5e-7`.")

ext: jl


text:
module FrostHelperEasedCameraZoomTrigger

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/EasedCameraZoomTrigger" EasedCameraZoomTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16, easing::String="Linear", targetZoom::Number=2.0, revertOnLeave::Bool=false, easeDuration::Number=1.0, focusOnPlayer::Bool=false, revertMode::String="RevertToNoZoom", disableInPhotosensitiveMode::Bool=false)

const placements = Ahorn.PlacementDict(
    "Eased Camera Zoom (Frost Helper)" => Ahorn.EntityPlacement(
        EasedCameraZoomTrigger,
        "rectangle"
    )
)

Ahorn.nodeLimits(entity::EasedCameraZoomTrigger) = 0, 1

Ahorn.editingOptions(entity::EasedCameraZoomTrigger) = Dict{String, Any}(
    "easing" => easings,
    "revertMode" => revertModes
)

const easings = String[
    "BackIn",
    "BackInOut",
    "BackOut",
    "BounceIn",
    "BounceInOut",
    "BounceOut",
    "CubeIn",
    "CubeInOut",
    "CubeOut",
    "ElasticIn",
    "ElasticInOut",
    "ElasticOut",
    "ExpoIn",
    "ExpoInOut",
    "ExpoOut",
    "Linear",
    "QuadIn",
    "QuadInOut",
    "QuadOut",
    "SineIn",
    "SineInOut",
    "SineOut"
]

const revertModes = String[
    "RevertToNoZoom",
    "RevertToPreviousZoom",
]

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using PyCall: python

run(`$python -m pip install --upgrade pip`)
run(`$python -m pip install jqdatasdk`)


text:
function forward_nolog!(α, c, obs_density::AbstractMatrix, hmm::HMM)
    T, S = size(obs_density, 1), nb_states(hmm)
    π0, P = initial_distribution(hmm), transition_matrix(hmm)
    for i = 1:S
        α[1, i] = π0[i] * obs_density[1, i]
    end
    c[1] = sum(@view α[1, :])  # scaling
    for i = 1:S
        α[1, i] /= c[1]
    end
    for t = 1:T-1
        for j = 1:S
            α[t+1, j] = sum(α[t, i] * P[i, j] for i = 1:S) * obs_density[t+1, j]
        end
        c[t+1] = sum(@view α[t+1, :])  # scaling
        for j = 1:S
            α[t+1, j] /= c[t+1]
        end
    end
    for t = 1:T
        if all_zeros(@view α[t, :])
            throw(OverflowError("Probabilities are too small in forward step."))
        end
    end
    return nothing
end

function forward_log!(logα, obs_logdensity::AbstractMatrix, hmm::HMM)
    T, S = size(obs_logdensity, 1), nb_states(hmm)
    logπ0, logP = log.(initial_distribution(hmm)), log.(transition_matrix(hmm))
    for i = 1:S
        logα[1, i] = logπ0[i] + obs_logdensity[1, i]
    end
    for t = 1:T-1
        for j = 1:S
            logα[t+1, j] = logsumexp(logα[t, :] + logP[:, j]) + obs_logdensity[t+1, j]  # TODO
        end
    end
    for t = 1:T
        if all_nan(@view logα[t, :])
            throw(OverflowError("Log probabilities are too small in forward step."))
        end
    end
    return nothing
end

function backward_nolog!(β, c, obs_density::AbstractMatrix, hmm::HMM)
    T, S = size(obs_density, 1), nb_states(hmm)
    P = transition_matrix(hmm)

    for i = 1:S
        β[T, i] = 0.0
    end
    for t = T-1:-1:1
        for i = 1:S
            β[t, i] = sum(P[i, j] * obs_density[t+1, j] * β[t+1, j] for j = 1:S) / c[t]
        end
    end
    for t = 1:T
        if all_zeros(@view β[t, :])
            throw(OverflowError("Log probabilities are too small in backward step."))
        end
    end
    return nothing
end

function backward_log!(logβ, obs_logdensity::AbstractMatrix, hmm::HMM)
    T, S = size(obs_logdensity, 1), nb_states(hmm)
    logP = log.(transition_matrix(hmm))

    for i = 1:S
        logβ[T, i] = 0.0
    end
    for t = T-1:-1:1
        for i = 1:S
            logβ[t, i] = logsumexp(logP[i, :] .+ obs_logdensity[t+1, :] .+ logβ[t+1, :])  # TODO
        end
    end
    for t = 1:T
        if all_nan(@view logβ[t, :])
            throw(OverflowError("Log probabilities are too small in backward step."))
        end
    end
    return nothing
end

function forward_backward_nolog!(α, β, c, γ, ξ, obs_density::AbstractMatrix, hmm::HMM)
    T, S = size(obs_density, 1), nb_states(hmm)

    forward_nolog!(α, c, obs_density, hmm)
    backward_nolog!(β, c, obs_density, hmm)

    for t = 1:T
        for i = 1:S
            γ[t, i] = α[t, i] * β[t, i]
        end
        sumγ = sum(@view γ[t, :])
        for i = 1:S
            γ[t, i] /= sumγ
        end
    end

    for t = 1:T-1
        for i = 1:S, j = 1:S
            ξ[t, i, j] = α[t, i] * P[i, j] * obs_density[t+1, j] * β[t+1, j]
        end
        sumξ = sum(@view ξ[t, :, :])
        for i = 1:S, j = 1:S
            ξ[t, i, j] /= sumξ
        end
    end

    logL = sum(log.(c))  # TODO: @view

    return log(L)
end

function forward_backward_log!(
    logα,
    logβ,
    logγ,
    logξ,
    obs_logdensity::AbstractMatrix,
    hmm::HMM,
)
    T, S = size(obs_logdensity, 1), nb_states(hmm.transitions)
    logP = log.(transition_matrix(hmm))

    forward_log!(logα, obs_logdensity, hmm)
    backward_log!(logβ, obs_logdensity, hmm)

    for t = 1:T
        for i = 1:S
            logγ[t, i] = logα[t, i] + logβ[t, i]
        end
        logsumγ = logsumexp(@view logγ[t, :])
        for i = 1:S
            logγ[t, i] -= logsumγ
        end
    end

    for t = 1:T-1
        for i = 1:S, j = 1:S
            logξ[t, i, j] = logα[t, i] + logP[i, j] + obs_logdensity[t+1, j] + logβ[t+1, j]
        end
        logsumξ = logsumexp(@view logξ[t, :, :])
        for i = 1:S, j = 1:S
            logξ[t, i, j] -= logsumξ
        end
    end

    logL = logsumexp(logα[T, :])

    return logL
end

## Likelihood of observations

function update_obs_density!(
    obs_density::AbstractMatrix,
    hmm::HMM,
    observations::AbstractVector,
)
    T, S = length(observations), nb_states(hmm)
    for t = 1:T
        for s = 1:S
            obs_density[t, s] = density(emission(hmm, s), observations[t])
        end
    end
    for t = 1:T
        if all_zeros(@view obs_logdensity[t, :])
            throw(OverflowError("Densities are too small for observations."))
        end
    end
end

function update_obs_logdensity!(
    obs_logdensity::AbstractMatrix,
    hmm::HMM,
    observations::AbstractVector,
)
    T, S = length(observations), nb_states(hmm)
    for t = 1:T
        for s = 1:S
            obs_logdensity[t, s] = logdensity(emission(hmm, s), observations[t])
        end
    end
    for t = 1:T
        if all_minus_inf(@view obs_logdensity[t, :])
            throw(OverflowError("Log-densities are too small for observations."))
        end
    end
end

function baum_welch_nolog!(
    α,
    β,
    c,
    γ,
    ξ,
    obs_density::AbstractMatrix,
    hmm::HMM{Tr,Em},
    observations::AbstractVector;
    iterations,
) where {Tr,Em}
    logL_evolution = Float64[]
    for _ = 1:iterations
        update_obs_density!(obs_density, hmm, observations)
        logL = forward_backward_nolog!(α, β, c, γ, ξ, obs_density, hmm)
        push!(logL_evolution, logL)
        new_transitions = fit(Tr, γ = γ, ξ = ξ)
        new_emissions = [fit(Em, observations, γ[:, s]) for s = 1:nb_states(hmm)]  # TODO: @view
        hmm = HMM(new_transitions, new_emissions)
    end
    return hmm, logL_evolution
end

function baum_welch_log!(
    logα,
    logβ,
    logγ,
    logξ,
    obs_logdensity::AbstractMatrix,
    hmm::HMM{Tr,Em},
    observations::AbstractVector;
    iterations,
) where {Tr,Em}
    logL_evolution = Float64[]
    for _ = 1:iterations
        update_obs_logdensity!(obs_logdensity, hmm, observations)
        logL = forward_backward_log!(logα, logβ, logγ, logξ, obs_logdensity, hmm)
        push!(logL_evolution, logL)
        new_transitions = fit(Tr, exp.(logγ), exp.(logξ))
        new_emissions = [fit(Em, observations, exp.(logγ[:, s])) for s = 1:nb_states(hmm)]  # TODO: @view
        hmm = HMM(new_transitions, new_emissions)
    end
    return hmm, logL_evolution
end

## Non mutating version

function baum_welch(hmm::HMM, observations::AbstractVector; iterations, log = true)
    T, S = length(observations), nb_states(hmm)
    if log
        logα = Matrix{Float64}(undef, T, S)
        logβ = Matrix{Float64}(undef, T, S)
        logγ = Matrix{Float64}(undef, T, S)
        logξ = Array{Float64,3}(undef, T - 1, S, S)
        obs_logdensity = Matrix{Float64}(undef, T, S)
        hmm_est, logL_evolution = baum_welch_log!(
            logα,
            logβ,
            logγ,
            logξ,
            obs_logdensity,
            hmm,
            observations;
            iterations = iterations,
        )
    else
        α = Matrix{Float64}(undef, T, S)
        β = Matrix{Float64}(undef, T, S)
        c = Vector{Float64}(undef, T)
        γ = Matrix{Float64}(undef, T, S)
        ξ = Array{Float64,3}(undef, T - 1, S, S)
        obs_density = Matrix{Float64}(undef, T, S)
        hmm_est, logL_evolution = baum_welch_nolog!(
            α,
            β,
            c,
            γ,
            ξ,
            obs_density,
            hmm,
            observations;
            iterations = iterations,
        )
    end
    return hmm_est, logL_evolution
end


text:
# Statistical classification
struct BoxPlot <: AbstractStatisticalMapClassificator
    k::Int
    h::Float64
    BoxPlot(h::Float64 = 1.5) = return new(6, h);
end

struct StdMean <: AbstractStatisticalMapClassificator
    k::Int
    StdMean() = return new(6);
end

struct Percentiles <: AbstractStatisticalMapClassificator
    k::Int
    p::Union{Vector{Int}, Vector{Float64}}
    Percentiles(p::Union{Vector{Int}, Vector{Float64}} = [1, 10, 50, 90, 99]) = return new(length(p) + 1, p);
end

# Bounds
function mapclassifybounds(mcr::BoxPlot, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    h = mcr.h
    q = quantile(x, [0.25, 0.5, 0.75])
    q25 = q[1]
    q50 = q[2]
    q75 = q[3]
    IQR = q75 - q25

    lbound = [-Inf,
              q25 - h * IQR,
              q25,
              q50,
              q75,
              q75 + h * IQR]
    ubound = [q25 - h * IQR,
              q25,
              q50,
              q75,
              q75 + h * IQR,
              Inf]

    lbound, ubound
end

function mapclassifybounds(::StdMean, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    m = mean(x)
    sd = std(x, corrected = true)

    lbounds = [-Inf, -2, -1, 0, 1, 2] * sd .+ m
    ubounds = [-2, -1, 0, 1, 2, Inf]  * sd .+ m

    lbounds, ubounds
end

function mapclassifybounds(mcr::Percentiles, x::Vector{T} where T<:Real)::Tuple{Vector{Float64}, Vector{Float64}}
    p = percentile(x, mcr.p)

    lbound = vcat(minimum(x), p)
    ubound = vcat(p, maximum(x))

    lbound, ubound
end


text:
__precompile__()

module DiffEqSensitivity

using DiffEqBase, Compat, ForwardDiff, Flux, DiffEqDiffTools, Statistics
using DiffEqCallbacks, QuadGK, RecursiveArrayTools, LinearAlgebra
using DataFrames, GLM, RecursiveArrayTools

abstract type SensitivityFunction end

include("derivative_wrappers.jl")
include("local_sensitivity.jl")
include("adjoint_sensitivity.jl")
include("morris_sensitivity.jl")
include("sobol_sensitivity.jl")
include("regression_sensitivity.jl")
include("DGSM.jl")

export extract_local_sensitivities

export ODELocalSensitivityFunction, ODELocalSensitivityProblem, SensitivityFunction,
       ODEAdjointSensitivityProblem, ODEAdjointProblem, AdjointSensitivityIntegrand,
       adjoint_sensitivities, adjoint_sensitivities_u0,
       morris_sensitivity, MorrisSensitivity, sobol_sensitivity,
       SensitivityAlg, regression_sensitivity, DGSM

end # module


text:
struct HydroDispatchRunOfRiverCascade <: PSI.AbstractHydroDispatchFormulation end

# Constraint name
const CASCADING_FLOW = "cascading_flow"

function flow_balance_external_input_param_cascade(
    optimization_container::PSI.OptimizationContainer,
    inflow_data::Vector{PSI.DeviceTimeSeriesConstraintInfo},
    upstream::Vector{
        Vector{NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}},
    },
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol},
    param_reference::PSI.UpdateRef,
)
    time_steps = PSI.model_time_steps(optimization_container)
    resolution = PSI.model_resolution(optimization_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / 60.0

    name_index = [PSI.get_component_name(d) for d in inflow_data]

    varspill = PSI.get_variable(optimization_container, var_names[1])
    varout = PSI.get_variable(optimization_container, var_names[2])

    container_inflow = PSI.add_param_container!(
        optimization_container,
        param_reference,
        name_index,
        time_steps,
    )
    param_inflow = PSI.get_parameter_array(container_inflow)
    multiplier_inflow = PSI.get_multiplier_array(container_inflow)

    flow_constraint =
        PSI.add_cons_container!(optimization_container, cons_name, name_index, time_steps)

    for (ix, d) in enumerate(inflow_data)
        name = PSI.get_component_name(d)
        upstream_devices = upstream[ix]

        multiplier_inflow[name, 1] = d.multiplier
        param_inflow[name, 1] =
            PJ.add_parameter(optimization_container.JuMPmodel, d.timeseries[1])

        exp =
            multiplier_inflow[name, 1] * param_inflow[name, 1] - varspill[name, 1] -
            varout[name, 1] #+ initial_conditions[ix].value
        #= spillage isn't an initial condition because it's not stored in the struct, so we can't make the cascading flow constrainits for the first periiod.
            for j in upstream_devices
                JuMP.add_to_expression!(exp, varspill[IS.get_name(j), 1])
                JuMP.add_to_expression!(exp, varout[IS.get_name(j), 1])
            end
        =#
        flow_constraint[name, 1] =
            JuMP.@constraint(optimization_container.JuMPmodel, exp == 0.0)

        for t in time_steps[2:end]
            param_inflow[name, t] =
                PJ.add_parameter(optimization_container.JuMPmodel, d.timeseries[t])

            exp = d.multiplier * param_inflow[name, t] - varspill[name, t] - varout[name, t]
            for j in upstream_devices
                if t - j.lag >= 1
                    JuMP.add_to_expression!(
                        exp,
                        varspill[IS.get_name(j.unit), t - j.lag],
                        j.multiplier,
                    )
                    JuMP.add_to_expression!(
                        exp,
                        varout[IS.get_name(j.unit), t - j.lag],
                        j.multiplier,
                    )
                end
            end
            flow_constraint[name, t] =
                JuMP.@constraint(optimization_container.JuMPmodel, exp == 0.0)
        end
    end

    return
end

function flow_balance_external_input_cascade(
    optimization_container::PSI.OptimizationContainer,
    inflow_data::Vector{PSI.DeviceTimeSeriesConstraintInfo},
    upstream::Vector{
        Vector{NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}},
    },
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol},
)
    time_steps = PSI.model_time_steps(optimization_container)
    resolution = PSI.model_resolution(optimization_container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / 60.0

    name_index = [PSI.get_component_name(d) for d in inflow_data]

    varspill = PSI.get_variable(optimization_container, var_names[1])
    varout = PSI.get_variable(optimization_container, var_names[2])

    flow_constraint =
        PSI.add_cons_container!(optimization_container, cons_name, name_index, time_steps)

    for (ix, d) in enumerate(inflow_data)
        name = PSI.get_component_name(d)
        upstream_devices = upstream[ix]

        exp = d.multiplier * d.timeseries[1] - varspill[name, 1] - varout[name, 1] #+ initial_conditions[ix].value
        #= spillage isn't an initial condition because it's not stored in the struct, so we can't make the cascading flow constrainits for the first periiod.
            for j in upstream_devices
                JuMP.add_to_expression!(exp, varspill[IS.get_name(j), 1])
                JuMP.add_to_expression!(exp, varout[IS.get_name(j), 1])
            end
        =#
        flow_constraint[name, 1] =
            JuMP.@constraint(optimization_container.JuMPmodel, exp == 0.0)

        for t in time_steps[2:end]
            exp = d.multiplier * d.timeseries[t] - varspill[name, t] - varout[name, t]

            for j in upstream_devices
                if t - j.lag >= 1
                    JuMP.add_to_expression!(
                        exp,
                        varspill[IS.get_name(j.unit), t - j.lag],
                        j.multiplier,
                    )
                    JuMP.add_to_expression!(
                        exp,
                        varout[IS.get_name(j.unit), t - j.lag],
                        j.multiplier,
                    )
                end
            end
            flow_constraint[name, t] =
                JuMP.@constraint(optimization_container.JuMPmodel, exp == 0.0)
        end
    end
    return
end

function flow_balance_cascade_constraint!(
    optimization_container::PSI.OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{H},
    model::PSI.DeviceModel{H, HydroDispatchRunOfRiverCascade},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, PSI.AbstractAffectFeedForward},
) where {H <: HydroCascade}
    parameters = PSI.model_has_parameters(optimization_container)
    use_forecast_data = PSI.model_uses_forecasts(optimization_container)

    inflow_forecast_label = "max_active_power"
    constraint_infos_inflow =
        Vector{PSI.DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    upstream_data = Vector{
        Vector{NamedTuple{(:unit, :lag, :multiplier), Tuple{PSY.HydroGen, Int64, Float64}}},
    }(
        undef,
        length(devices),
    )

    for (ix, d) in enumerate(devices)
        ts_vector_inflow =
            PSI.get_time_series(optimization_container, d, inflow_forecast_label)
        constraint_info_inflow = PSI.DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_max_active_power(x),
            ts_vector_inflow,
        )
        PSI.add_device_services!(constraint_info_inflow.range, d, model)
        constraint_infos_inflow[ix] = constraint_info_inflow
        upstream_data[ix] = get_upstream(d)
    end

    if parameters
        flow_balance_external_input_param_cascade(
            optimization_container,
            constraint_infos_inflow,
            upstream_data,
            PSI.make_constraint_name(CASCADING_FLOW, H),
            (
                PSI.make_variable_name(PSI.SPILLAGE, H),
                PSI.make_variable_name(PSI.ACTIVE_POWER, H),
            ),
            PSI.UpdateRef{H}(PSI.INFLOW, inflow_forecast_label),
        )
    else
        flow_balance_external_input_cascade(
            optimization_container,
            constraint_infos_inflow,
            upstream_data,
            PSI.make_constraint_name(CASCADING_FLOW, H),
            (
                PSI.make_variable_name(PSI.SPILLAGE, H),
                PSI.make_variable_name(PSI.ACTIVE_POWER, H),
            ),
        )
    end
    return
end


text:
import Adapt

using Oceananigans.Utils: prettysummary

"""
    struct DiscreteForcing{P, F}

Wrapper for "discrete form" forcing functions with optional `parameters`.
"""
struct DiscreteForcing{P, F}
    func :: F
    parameters :: P
end

"""
    DiscreteForcing(func; parameters=nothing)

Construct a "discrete form" forcing function with optional parameters.
The forcing function is applied at grid point `i, j, k`.

When `parameters` are not specified, `func` must be callable with the signature

```
func(i, j, k, grid, clock, model_fields)
```

where `grid` is `model.grid`, `clock.time` is the current simulation time and
`clock.iteration` is the current model iteration, and `model_fields` is a
`NamedTuple` with `u, v, w` and the fields in `model.tracers`.

*Note* that the index `end` does *not* access the final physical grid point of
a model field in any direction. The final grid point must be explicitly specified, as
in `model_fields.u[i, j, grid.Nz]`.

When `parameters` _is_ specified, `func` must be callable with the signature.

```
func(i, j, k, grid, clock, model_fields, parameters)
```

Above, `parameters` is, in principle, arbitrary. Note, however, that GPU compilation
can place constraints on `typeof(parameters)`.
"""
DiscreteForcing(func; parameters=nothing) = DiscreteForcing(func, parameters)

@inline function (forcing::DiscreteForcing{P, F})(i, j, k, grid, clock, model_fields) where {P, F<:Function}
    parameters = forcing.parameters
    return forcing.func(i, j, k, grid, clock, model_fields, parameters)
end

@inline (forcing::DiscreteForcing{<:Nothing, F})(i, j, k, grid, clock, model_fields) where F<:Function =
    forcing.func(i, j, k, grid, clock, model_fields)

"""Show the innards of a `DiscreteForcing` in the REPL."""
Base.show(io::IO, forcing::DiscreteForcing{P}) where P =
    print(io, "DiscreteForcing{$P}", '\n',
        "├── func: $(prettysummary(forcing.func))", '\n',
        "└── parameters: $(forcing.parameters)")

Adapt.adapt_structure(to, forcing::DiscreteForcing) =
    DiscreteForcing(Adapt.adapt(to, forcing.func),
                    Adapt.adapt(to, forcing.parameters))


text:
module Met_Plots_Models

# package code goes here
using Met_Basic_Types
using NetCDF
using PyCall
using PyPlot
export PyPlot
@pyimport mpl_toolkits.basemap as basemap
@pyimport prettyplotlib as ppl
@pyimport brewer2mpl

include("colormaps/cmaps.jl")

export zero_to_one_ratio_cmap;
export purple_cmap, rev_purple_cmap
export sst_only_positive_cmap;
export sst_only_negative_cmap, sst_cmap, diff_cmap, pp_cmap;

include("pairs/all.jl")

#= export mapplot_cross_compare =#

function flush_figures(fpath :: AbstractString)
  println("$fpath.png")
  savefig("$fpath.png")
  savefig("$fpath.pdf")
  savefig("$fpath.svg")
end
export flush_figures

end # module


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module BernoulliNumbers
using Nemo, NumberTheory, PrimesIterator, AndreNumbers, ClausenNumbers, Products

export ModuleBernoulliNumbers
export BernoulliInt, BernoulliIntList, Bernoulli, BernoulliList
export V195441, V065619, V281586, V281588, V027641, L065619

"""
We are primarily concerned with the integer Bernoulli numbers.

Cf. ``A000182 (m=2), A293951 (m=3), A273352 (m=4), A318258 (m=5).``

```
[1] [0, 1,    0,       0,             0,                  0]
[2] [0, 1,   -2,      16,          -272,               7936]
[3] [0, 1,   -9,     477,        -74601,           25740261]
[4] [0, 1,  -34,   11056,     -14873104,        56814228736]
[5] [0, 1, -125,  249250,   -2886735625,    122209131374375]
[6] [0, 1, -461, 5699149, -574688719793, 272692888959243481]
```

The rational Bernoulli numbers are defined here with ``B(1) = 1/2``. Why this is preferred over ``B(1) = -1/2`` is explained in the [Bernoulli Manifesto](http://luschny.de/math/zeta/The-Bernoulli-Manifesto.html).

* BernoulliInt, BernoulliIntList, Bernoulli, BernoulliList
* V195441, V065619, V281586, V281588, V027641, L065619
"""
const ModuleBernoulliNumbers = ""

"""
Return the generalized integer Bernoulli numbers ``b_{m}(n) = n × ``André``(m, n-1)``.
"""
BernoulliInt(m::Int, n::Int) = n == 0 ? ZZ(0) : n * André(m, n - 1)

"""
Return the number of down-up permutations w on ``[n+1]`` such that ``w_2 = 1``.
"""
V065619(n::Int) = BernoulliInt(2, n)

"""
Return the generalized integer Bernoulli numbers ``b_{3}(n) = n × ``André``_{3}(n-1)``.
"""
V281586(n::Int) = BernoulliInt(3, n)

"""
Return the generalized integer Bernoulli numbers ``b_{4}(n) = n × `` André``_{4}(n-1)``.
"""
V281588(n::Int) = BernoulliInt(4, n)

"""
Return a list of length `len` of the integer Bernoulli numbers ``b_{m}(n)`` using Seidel's boustrophedon algorithm.
"""
function BernoulliIntList(m::Int, len::Int)
    len ≤ 0 && return fmpz[]
    R = zeros(ZZ, len)
    len == 1 && return R
    R[2] = 1
    len == 2 && return R
    A = zeros(ZZ, len)
    A[1] = 1; A[2] = 1

    for n in 1:len - 2
        if n % m ≠ 0
            for i in n:-1:1 A[i] += A[i + 1] end
            C = A[1]
        else
            C = 0
            for i in 1:(n + 2) A[i], C = C, A[i]; C = A[i] - C end
        end
        R[n + 2] = C
    end
    R
end

"""
Computes a list of length `len` of the integer Bernoulli numbers ``b_{2}(n)`` using Seidel's boustrophedon algorithm.
"""
function L065619(len::Int)
    len ≤ 0  && return fmpz[]
    R  = zeros(ZZ, len)
    len == 1 && return R[0]
    len == 2 && (R[0, 1] = 1; return R)

    A = Dict{Int,fmpz}(-1 => 1, 0 => 0)
    k = 0; e = 1
    for i in 0:len - 1
        Am = 0; A[k + e] = 0; e = -e
        for j in 0:i
            Am += A[k]; A[k] = Am; k += e
        end
        j = e < 0 ? div(-i, 2) : div(i, 2)
        R[i+1] = A[j]
    end
    R
end

"""
Return the rational Bernoulli number ``B_n``  (cf. A027641/A027642).
"""
function Bernoulli(n::Int)
    isOdd(n) && (n == 1 ? (return fmpq(1, 2)) : (return fmpq(0, 1)))
    n == 0 && return fmpq(1, 1)
    denom = ^(ZZ(4), n) - ^(ZZ(2), n)
    fmpq(BernoulliInt(2, n), denom)
end

"""
Return a list of the first `len` Bernoulli numbers ``B_n`` (cf. A027641/A027642).
"""
function BernoulliList(len::Int)
    if len ≤ 0 return fmpq[] end
    R = zeros(QQ, len)
    R[1] = fmpq(1, 1); len == 1 && return R
    R[2] = fmpq(1, 2); len == 2 && return R

    A = Dict{Int,fmpz}(0 => 1, -2 => 0, -1 => 1, 1 => 0)
    a = fmpz(12); b = fmpz(240)
    k = e = 1

    for i in 2:len - 1
        Am = 0; A[k + e] = 0; e = -e
        for j in 0:i
            Am += A[k]; A[k] = Am; k += e
        end
        if e > 0
            R[i + 1] = fmpq(0, 1)
        else
            d = i >> 1
            R[i + 1] = isEven(d) ? fmpq(-A[-d], a) : fmpq(A[-d], a)
            a, b = b, b << 4 + b << 2 - a << 6
        end
    end
    R
end

"""
Return the numerator of the Bernoulli number ``B_n``.
"""
function V027641(n::Int)
    isOdd(n) && (n == 1 ? (return ZZ(-1)) : return ZZ(0))
    n == 0 && return ZZ(1)
    denom = ^(ZZ(4), n) - ^(ZZ(2), n)
    Nemo.numerator(BernoulliInt(2, n) // denom)
end

# We could also define the denominator of the Bernoulli number as
#    V027642(n::Int) = denominator(Bernoulli(n))
# however V027642 is implemented more efficiently via the
# Clausen numbers (see the module "Clausen").

"""
Return denominator(Bernoulli ``_{n+1}(x) - `` Bernoulli ``_{n+1})``.
"""
function V195441(n::Int)
    n < 4 && return ZZ([1, 1, 2, 1][n + 1])
    P = Primes(2, div(n + 2, 2 + n % 2))
    ∏([p for p in P if p ≤ sum(digits(n + 1, base=Int(p)))])
end

#START-TEST-########################################################

using Test, SeqTests

function test()

    @testset "BernoulliNum" begin

        @test BernoulliInt(2, 10) == ZZ(79360)
        @test BernoulliInt(3, 30) == ZZ(-7708110416280010548302670)
        @test BernoulliInt(4, 40) == ZZ(-44494882577309421077208834962882560)

        @test isa(BernoulliInt(3, 30), Nemo.fmpz)
        @test isa(BernoulliIntList(2, 20)[end], Nemo.fmpz)
        @test isa(Bernoulli(0), Nemo.fmpq)

        @test isa(V027641(10), Nemo.fmpz)
        @test isa(V195441(10), Nemo.fmpz)

        for m in 1:20
            @test BernoulliInt(m, 200) == BernoulliIntList(m, 200 + 1)[end]
        end

        for n in 0:20
            @test denominator(Rational(Bernoulli(2 * n))) == ClausenNumber(n)
        end

        t = ZZ(8622529719094842064796322984685715031642180319435676189471082876882585178647210)
        @test V195441(10000) == t

        # A065619 in the OEIS has an arbitrary offset of 1.
        l = fmpz[0, 1, 2, 3, 8, 25, 96, 427, 2176]
        @test all(L065619(9) .== l)

        if is_oeis_installed()
            V = [V195441, V281586, V281588, V027641]
            for v in V SeqTest(v, 'V') end
        end
    end
end

function demo()
    println("\nGeneralized integer Bernoulli numbers BernoulliInt(m, n):")
    for m in 1:5
        println([BernoulliInt(m, n) for n in 0:10])
    end
    println("\nGeneralized integer Bernoulli numbers BernoulliInt(3, n):")
    for len in 1:6
        println(BernoulliIntList(3, len))
    end
    println("\nClassical rational Bernoulli numbers:")
    for n in 0:20
        print(fmpq(Nemo.numerator(Bernoulli(n)), Nemo.denominator(Bernoulli(n))), ", ")
    end
    println()
end

"""
Since the BernoulliInts are cached via the André numbers the benchmarks below
depend on the cache history.

for m in 1:10, n in 0:499 BernoulliInt(m,n) end
    0.983624 seconds (1.67 M allocations: 36.160 MiB, 9.71% gc time)
for m in 1:10 BernoulliIntList(m,500) end
    0.504717 seconds (1.26 M allocations: 19.258 MiB)
for n in 0:1000 Bernoulli(n) end
    1.104452 seconds (946.74 k allocations: 20.252 MiB)
BernoulliList(1000)
    0.789146 seconds (508.75 k allocations: 7.891 MiB, 11.06% gc time)
for n in 0:1000 V027641(n) end
    0.005921 seconds (7.74 k allocations: 144.313 KiB)
for n in 0:10000 V195441(n) end
    0.881609 seconds (3.49 M allocations: 330.094 MiB, 26.84% gc time)
"""
function perf()
    GC.gc()
    @time (for m in 1:10, n in 0:499 BernoulliInt(m, n) end)
    @time (for m in 1:10 BernoulliIntList(m, 500) end)

    GC.gc()
    @time (for n in 0:1000 Bernoulli(n) end)
    @time BernoulliList(1000)
    @time (for n in 0:1000 V027641(n) end)

    GC.gc()
    @time (for n in 0:10000 V195441(n) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

# [0] 0, 1,  2, 3,  4,  5,    6,    7,     8,     9,     10
# [1] 0, 1, -2, 3, -4,  5,   -6,    7,    -8,     9,    -10
# [2] 0, 1, 2, -3, -8,  25,  96, -427, -2176, 12465,  79360
# [3] 0, 1, 2,  3, -4, -15, -54,  133,   792,  4293, -15130
# [4] 0, 1, 2,  3,  4,  -5, -24,  -98,  -272,   621,   4960
# [5] 0, 1, 2,  3,  4,   5,  -6,  -35,  -160,  -495,  -1250


text:
using DifferentialEquations, Test

f = function (du,u,p,t)
  du[1] = 2 - 2u[1]
  du[2] = u[1] - 4u[2]
end
u0 = zeros(2)
prob = SteadyStateProblem(f,u0)

sol = solve(prob)

@test typeof(sol.alg)<: SSRootfind

max_stars_repo_name: ChrisRackauckas/DiffEq


text:
## random number generation

using Random
import RandomNumbers

# global state

# shared memory with the actual seed, per warp, loaded lazily or overridden by calling `seed!`
@eval @inline function global_random_keys()
    ptr = Base.llvmcall(
        $("""@global_random_keys = weak addrspace($(AS.Shared)) global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Shared))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Shared))* @global_random_keys, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Shared))* %ptr to i8 addrspace($(AS.Shared))*
                 ret i8 addrspace($(AS.Shared))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Shared}, Tuple{})
    CuDeviceArray((32,), ptr)
end

# shared memory with per-warp counters, incremented when generating numbers
@eval @inline function global_random_counters()
    ptr = Base.llvmcall(
        $("""@global_random_counters = weak addrspace($(AS.Shared)) global [32 x i32] zeroinitializer, align 32
             define i8 addrspace($(AS.Shared))* @entry() #0 {
                 %ptr = getelementptr inbounds [32 x i32], [32 x i32] addrspace($(AS.Shared))* @global_random_counters, i64 0, i64 0
                 %untyped_ptr = bitcast i32 addrspace($(AS.Shared))* %ptr to i8 addrspace($(AS.Shared))*
                 ret i8 addrspace($(AS.Shared))* %untyped_ptr
             }
             attributes #0 = { alwaysinline }
          """, "entry"), LLVMPtr{UInt32, AS.Shared}, Tuple{})
    CuDeviceArray((32,), ptr)
end

@device_override Random.make_seed() = clock(UInt32)

# generators

using Random123: philox2x_round, philox2x_bumpkey

# GPU-compatible/optimized version of the generator from Random123.jl
struct Philox2x32{R} <: RandomNumbers.AbstractRNG{UInt64}
    @inline function Philox2x32{R}() where R
        rng = new{R}()
        if rng.key == 0
            # initialize the key. this happens when first accessing the (0-initialized)
            # shared memory key from each block. if we ever want to make the device seed
            # controlable from the host, this would be the place to read a global seed.
            #
            # note however that it is undefined how shared memory persists across e.g.
            # launches, so we may not be able to rely on the zero initalization then.
            rng.key = Random.make_seed()
        end
        return rng
    end
end

# default to 7 rounds; enough to pass SmallCrush
@inline Philox2x32() = Philox2x32{7}()

@inline function Base.getproperty(rng::Philox2x32, field::Symbol)
    threadId = threadIdx().x + (threadIdx().y - 1i32) * blockDim().x +
                               (threadIdx().z - 1i32) * blockDim().x * blockDim().y
    warpId = (threadId - 1i32) >> 0x5 + 1i32  # fld1

    if field === :seed
        @inbounds global_random_seed()[1]
    elseif field === :key
        @inbounds global_random_keys()[warpId]
    elseif field === :ctr1
        @inbounds global_random_counters()[warpId]
    elseif field === :ctr2
        blockId = blockIdx().x + (blockIdx().y - 1i32) * gridDim().x +
                                 (blockIdx().z - 1i32) * gridDim().x * gridDim().y
        globalId = threadId + (blockId - 1i32) * (blockDim().x * blockDim().y * blockDim().z)
        globalId%UInt32
    end::UInt32
end

@inline function Base.setproperty!(rng::Philox2x32, field::Symbol, x)
    threadId = threadIdx().x + (threadIdx().y - 1i32) * blockDim().x +
                               (threadIdx().z - 1i32) * blockDim().x * blockDim().y
    warpId = (threadId - 1i32) >> 0x5 + 1i32  # fld1

    if field === :key
        @inbounds global_random_keys()[warpId] = x
    elseif field === :ctr1
        @inbounds global_random_counters()[warpId] = x
    end
end

@device_override @inline Random.default_rng() = Philox2x32()

"""
    Random.seed!(rng::Philox2x32, seed::Integer, [counter::Integer=0])

Seed the on-device Philox2x32 generator with an UInt32 number.
Should be called by at least one thread per warp.
"""
function Random.seed!(rng::Philox2x32, seed::Integer, counter::Integer=0)
    rng.key = seed % UInt32
    rng.ctr1 = counter
    return
end

if VERSION >= v"1.7-"
@device_override Random.seed!(::Random._GLOBAL_RNG, seed) =
    Random.seed!(Random.default_rng(), seed)
end

"""
    Random.rand(rng::Philox2x32, UInt32)

Generate a byte of random data using the on-device Tausworthe generator.
"""
function Random.rand(rng::Philox2x32{R},::Type{UInt64}) where {R}
    ctr1, ctr2, key = rng.ctr1, rng.ctr2, rng.key

    if R > 0                               ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 1  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 2  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 3  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 4  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 5  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 6  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 7  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 8  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 9  key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 10 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 11 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 12 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 13 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 14 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end
    if R > 15 key = philox2x_bumpkey(key); ctr1, ctr2 = philox2x_round(ctr1, ctr2, key); end

    # update the warp counter
    # NOTE: this performs the same update on every thread in the warp, but each warp writes
    #       to a unique location so the duplicate writes are innocuous
    # XXX: what if this overflows? we can't increment ctr2. bump the key?
    rng.ctr1 += 1i32

    # NOTE: it's too expensive to keep both numbers around in case the user only wanted one,
    #       so just make our 2x32 generator return 64-bit numbers by default.
    return (ctr1 % UInt64) << 32 | (ctr2 % UInt64)
end

# normally distributed random numbers using Ziggurat algorithm
#
# copied from Base because we don't support its global tables

# a hacky method of exposing constant tables as constant GPU memory
function emit_constant_array(name::Symbol, data::AbstractArray{T}) where {T}
    Context() do ctx
        T_val = convert(LLVMType, T; ctx)
        T_ptr = convert(LLVMType, LLVMPtr{T,AS.Constant}; ctx)

        # define function and get LLVM module
        llvm_f, _ = create_function(T_ptr)
        mod = LLVM.parent(llvm_f)

        # create a global memory global variable
        # TODO: global_var alignment?
        T_global = LLVM.ArrayType(T_val, length(data))
        # XXX: why can't we use a single name like emit_shmem
        gv = GlobalVariable(mod, T_global, "gpu_$(name)_data", AS.Constant)
        linkage!(gv, LLVM.API.LLVMInternalLinkage)
        initializer!(gv, ConstantArray(data; ctx))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            ptr = gep!(builder, gv, [ConstantInt(0; ctx), ConstantInt(0; ctx)])

            untyped_ptr = bitcast!(builder, ptr, T_ptr)

            ret!(builder, untyped_ptr)
        end

        call_function(llvm_f, LLVMPtr{T,AS.Constant})
    end
end

for var in [:ki, :wi, :fi, :ke, :we, :fe]
    val = getfield(Random, var)
    gpu_var = Symbol("gpu_$var")
    @eval @inline $gpu_var() = CuDeviceArray($(size(val)), $(emit_constant_array(var, val)))
end

## randn

@device_override Random.randn(rng::AbstractRNG) =
    _randn(rng, Random.rand(rng, Random.UInt52Raw()))

@inline function _randn(rng::AbstractRNG, r::UInt64)
    @inbounds begin
        r &= 0x000fffffffffffff
        rabs = Int64(r>>1) # One bit for the sign
        idx = rabs & 0xFF
        x = ifelse(r % Bool, -rabs, rabs)*gpu_wi()[idx+1]
        rabs < gpu_ki()[idx+1] && return x # 99.3% of the time we return here 1st try
        return randn_unlikely(rng, idx, rabs, x)
    end
end

# this unlikely branch is put in a separate function for better efficiency
@noinline function randn_unlikely(rng, idx, rabs, x)
    @inbounds if idx == 0
        while true
            xx = -Random.ziggurat_nor_inv_r*log(Random.rand(rng))
            yy = -log(Random.rand(rng))
            yy+yy > xx*xx &&
                return (rabs >> 8) % Bool ? -Random.ziggurat_nor_r-xx : Random.ziggurat_nor_r+xx
        end
    elseif (gpu_fi()[idx] - gpu_fi()[idx+1])*Random.rand(rng) + gpu_fi()[idx+1] < exp(-0.5*x*x)
        return x # return from the triangular area
    else
        return Random.randn(rng)
    end
end

## randexp

@device_override Random.randexp(rng::AbstractRNG) =
    _randexp(rng, Random.rand(rng, Random.UInt52Raw()))

function _randexp(rng::AbstractRNG, ri::UInt64)
    @inbounds begin
        ri &= 0x000fffffffffffff
        idx = ri & 0xFF
        x = ri*gpu_we()[idx+1]
        ri < gpu_ke()[idx+1] && return x # 98.9% of the time we return here 1st try
        return randexp_unlikely(rng, idx, x)
    end
end

@noinline function randexp_unlikely(rng, idx, x)
    @inbounds if idx == 0
        return Random.ziggurat_exp_r - log(Random.rand(rng))
    elseif (gpu_fe()[idx] - gpu_fe()[idx+1])*Random.rand(rng) + gpu_fe()[idx+1] < exp(-x)
        return x # return from the triangular area
    else
        return Random.randexp(rng)
    end
end


text:
export DCAT

"""
`DCAT(A::AbstractOperator...)`

Block-diagonally concatenate `AbstractOperator`s.

```julia
julia> D = DCAT(HCAT(Eye(2),Eye(2)),DFT(3))
[[I,I],0;0,ℱ]  ℝ^2  ℝ^2  ℝ^4 -> ℝ^2  ℂ^3

julia> DCAT(Eye(10),Eye(10),FiniteDiff((4,4)))
DCAT  ℝ^10  ℝ^10  ℝ^(4, 4) -> ℝ^10  ℝ^10  ℝ^(3, 4)
```

To evaluate `DCAT` operators multiply them with a `Tuple` of `AbstractArray` of the correct domain size and type. The output will consist as well of a `Tuple` with the codomain type and size of the `DCAT`.

```julia
julia> D*ArrayPartition(ones(2),ones(2),ones(3))
([2.0, 2.0], Complex{Float64}[3.0+0.0im, 0.0+0.0im, 0.0+0.0im])

```

"""
struct DCAT{N,
	    L  <: NTuple{N,AbstractOperator},
	    P1 <: NTuple{N,Union{Int,Tuple}},
	    P2 <: NTuple{N,Union{Int,Tuple}}
	    } <: AbstractOperator
	A::L
	idxD::P1
	idxC::P2
	DCAT(A::L, idxD::P1, idxC::P2) where {L, P1, P2} = new{length(A), L, P1, P2}(A,idxD,idxC)
end

# Constructors
DCAT(A::AbstractOperator) = A

function DCAT(A::Vararg{AbstractOperator})

	# build H.idxs
	idxD,idxC = (),()
	for d in (1,2)
		idxs = []
		K = 0
		for i in eachindex(ndoms.(A,d))
			if ndoms(A[i],d) == 1 # flatten operator
				K += 1
				push!(idxs,K)
			else                   # stacked operator
				idxs = push!(idxs,(collect(K+1:K+ndoms(A[i],d))...,) )
				for ii = 1:ndoms(A[i],d)
					K += 1
				end
			end
		end
		d == 1 ? (idxC = (idxs...,)) : (idxD = (idxs...,))
	end

	return DCAT(A,idxD,idxC)

end

# Mappings
@generated function mul!(yy::ArrayPartition,
                         H::DCAT{N,L,P1,P2},
                         bb::ArrayPartition) where {N,L,P1,P2}

  # extract stuff
	ex = :(y = yy.x; b = bb.x )

	for i = 1:N

		if fieldtype(P2,i) <: Int
		# flatten operator
		# build mul!(y[H.idxC[i]], H.A[i], b)
			yy = :(y[H.idxC[$i]])
		else
		# stacked operator
		# build mul!(( y[H.idxC[i][1]], y[H.idxC[i][2]] ...  ), H.A[i], b)
        yy =  [ :(y[H.idxC[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P2,i)))]
        yy = :( ArrayPartition( $(yy...) ) )
		end

		if fieldtype(P1,i) <: Int
		# flatten operator
		# build mul!(H.buf, H.A[i], b[H.idxD[i]])
			bb = :(b[H.idxD[$i]])
		else
		# stacked operator
		# build mul!(H.buf, H.A[i],( b[H.idxD[i][1]], b[H.idxD[i][2]] ...  ))
        bb = [ :(b[H.idxD[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P1,i))) ]
        bb = :( ArrayPartition( $(bb...) ) )
		end

		ex = :($ex; mul!($yy,H.A[$i],$bb))

	end
	ex = :($ex; return y)
	return ex

end

@generated function mul!(yy::ArrayPartition,
                         A::AdjointOperator{DCAT{N,L,P1,P2}},
                         bb::ArrayPartition) where {N,L,P1,P2}

  # extract stuff
	ex = :(H = A.A; y = yy.x; b = bb.x )

	for i = 1:N

		if fieldtype(P1,i) <: Int
		# flatten operator
		# build mul!(y[H.idxD[i]], H.A[i]', b)
			yy = :(y[H.idxD[$i]])
		else
		# stacked operator
		# build mul!(( y[H.idxD[i][1]], y[H.idxD[i][2]] ...  ), H.A[i]', b)
        yy = [ :(y[H.idxD[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P1,i)))]
        yy = :( ArrayPartition( $(yy...) ))
		end

		if fieldtype(P2,i) <: Int
		# flatten operator
		# build mul!(H.buf, H.A[i]', b[H.idxC[i]])
			bb = :(b[H.idxC[$i]])
		else
		# stacked operator
		# build mul!(H.buf, H.A[i]',( b[H.idxC[i][1]], b[H.idxC[i][2]] ...  ))
        bb = [ :(b[H.idxC[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P2,i)))]
        bb = :( ArrayPartition( $(bb...) ) )
		end

		ex = :($ex; mul!($yy,H.A[$i]',$bb))

	end
	ex = :($ex; return y)
	return ex

end

# Properties
size(H::DCAT) = size(H,1),size(H,2)

function size(H::DCAT, i::Int)

	sz = []
	for s in size.(H.A,i)
		eltype(s) <: Int ? push!(sz,s) : push!(sz,s...)
	end
	p = vcat([[idx... ] for idx in (i == 1 ? H.idxC : H.idxD) ]...)
	invpermute!(sz,p)

	(sz...,)
end

fun_name(L::DCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*",0;0,"*fun_name(L.A[2])*"]" :
"DCAT"

function domainType(H::DCAT)
	domain = vcat([typeof(d)<:Tuple ? [d...] : d  for d in domainType.(H.A)]...)
	p = vcat([[idx... ] for idx in H.idxD]...)
	invpermute!(domain,p)
	return (domain...,)
end
function codomainType(H::DCAT)
	codomain = vcat([typeof(d)<:Tuple ? [d...] : d  for d in codomainType.(H.A)]...)
	p = vcat([[idx... ] for idx in H.idxC]...)
	invpermute!(codomain,p)
	return (codomain...,)
end

is_eye(L::DCAT) = all(is_eye.(L.A))
is_linear(L::DCAT) = all(is_linear.(L.A))
is_diagonal(L::DCAT) = all(is_diagonal.(L.A))
is_AcA_diagonal(L::DCAT) = all(is_AcA_diagonal.(L.A))
is_AAc_diagonal(L::DCAT) = all(is_AAc_diagonal.(L.A))
is_orthogonal(L::DCAT) = all(is_orthogonal.(L.A))
is_invertible(L::DCAT) = all(is_invertible.(L.A))
is_full_row_rank(L::DCAT) = all(is_full_row_rank.(L.A))
is_full_column_rank(L::DCAT) = all(is_full_column_rank.(L.A))

# utils
function permute(H::DCAT{N,L,P1,P2}, p::AbstractVector{Int}) where {N,L,P1,P2}

	unfolded = vcat([[idx... ] for idx in H.idxD]...)
	invpermute!(unfolded,p)

	new_part = ()
	cnt = 0
	for z in length.(H.idxD)
		new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...,))
		cnt += z
	end

	DCAT(H.A,new_part,H.idxC)
end

remove_displacement(D::DCAT) = DCAT(remove_displacement.(D.A), D.idxD, D.idxC)

# special cases
# Eye constructor
Eye(x::ArrayPartition) = DCAT(Eye.(x.x)...)
diag(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.
diag_AAc(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.
diag_AcA(L::DCAT{N,NTuple{N,E}}) where {N, E <: Eye} = 1.


text:
"""
Sets the tolerance for the iLQR solver to the intermediate tolerances, until the last
outer loop iteration.
"""
function set_tolerances!(solver::ALSolver{T}, i::Int, cost_tol=solver.opts.cost_tolerance,
                         grad_tol=solver.opts.gradient_tolerance) where T
    opts = options(solver)
    if i != solver.opts.iterations_outer
        opts.cost_tolerance = opts.cost_tolerance_intermediate
        opts.gradient_tolerance = opts.gradient_tolerance_intermediate
    else
        opts.cost_tolerance = cost_tol
        opts.gradient_tolerance = grad_tol
    end

    return nothing
end

function solve!(solver::ALSolver)
    reset!(solver)

    conset = get_constraints(solver)
    if !is_constrained(solver)
        ilqr = get_ilqr(solver)
        solve!(ilqr)
        terminate!(solver)
        return solver
    end

    # Terminal tolerances
    cost_tol = solver.opts.cost_tolerance
    grad_tol = solver.opts.gradient_tolerance

    Z̄ = solver.ilqr.Z̄
    for al_iter = 1:solver.opts.iterations_outer
        # Set potentially looser tolerances for inner iLQR solve
        set_tolerances!(solver, al_iter, cost_tol, grad_tol)

        # Solve problem with iLQR
        solve!(solver.ilqr)

        # Check solver status
        status(solver) > SOLVE_SUCCEEDED && break

        # Evaluate the constraints and the cost for the current trajectory
        J = TO.cost(solver, Z̄)  # NOTE: this evaluates the constraints
        c_max = max_violation(conset)
        μ_max = max_penalty(conset)
        record_iteration!(solver, J, c_max, μ_max)

        # Check if it's converged
        isconverged = evaluate_convergence(solver)
        if isconverged
            break
        end

        # Outer loop updates
        dualupdate!(conset)
        penaltyupdate!(conset)

        # Reset iLQR solver
        # TODO: is this necessary? it gets reset at the beginning of the iLQR solve method
        reset!(solver.ilqr)
    end
    # Reset tolerances to their original values
    solver.opts.cost_tolerance = cost_tol
    solver.opts.gradient_tolerance = grad_tol
    terminate!(solver)
    return solver
end

function record_iteration!(solver::ALSolver, J, c_max, μ_max)
    stats = solver.stats
    record_iteration!(stats, c_max=c_max, penalty_max=μ_max, is_outer=true)
    lg = getlogger(solver)
    @log lg "iter" stats.iterations
    @log lg "AL iter" stats.iterations_outer
    @log lg "cost" J
    @log lg "||v||" c_max
    @log lg μ_max
end

function evaluate_convergence(solver::ALSolver)
    lg = getlogger(solver)
    iter = solver.stats.iterations
    isconverged = false
    if solver.stats.c_max[iter] < solver.opts.constraint_tolerance
        @log lg "info" "Constraint tolerance met."
        solver.stats.status = SOLVE_SUCCEEDED
        isconverged = true
    end
    if solver.opts.kickout_max_penalty && solver.stats.penalty_max[i] >= solver.opts.penalty_max
        @log lg "info" "Hit max penalty."
        isconverged = true
    end
    if iter >= solver.opts.iterations
        @log lg "info" "Hit max iterations."
        solver.stats.status = MAX_ITERATIONS
        isconverged = true
    end
    if solver.stats.iterations_outer >= solver.opts.iterations_outer
        @log lg "info" "Hit max AL iterations."
        solver.stats.status = MAX_ITERATIONS_OUTER
        isconverged = true
    end
    return isconverged
end


text:
immutable QueryableWhere{T,Provider} <: Queryable{T,Provider}
	source
	filter::Expr
end

function where{T,Provider}(source::Queryable{T,Provider}, filter::Function, filter_expr::Expr)
    return QueryableWhere{T,Provider}(source, filter_expr)
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# Handles Taylor series of arbitrary but finite order

__precompile__(true)

"""
TaylorSeries

Main module for `TaylorSeries.jl` -- a package for Taylor expansions in one
or more independent variables.
"""
module TaylorSeries

using Compat
import Compat.String

import Base: ==, +, -, *, /, ^

import Base: zero, one, zeros, ones, isinf, isnan,
    convert, promote_rule, promote, eltype, length, show,
    real, imag, conj, ctranspose,
    rem, mod, mod2pi, abs, gradient,
    sqrt, exp, log, sin, cos, tan,
    asin, acos, atan, sinh, cosh, tanh,
    A_mul_B!

export Taylor1, TaylorN, HomogeneousPolynomial

export get_coeff, derivative, integrate,
    evaluate, evaluate!,
    show_params_TaylorN,
    get_order, get_numvars,
    set_variables, get_variables,
    ∇, jacobian, hessian

# one variable:
include("Taylor1.jl")

# several variables:
include("parameters.jl")
include("hash_tables.jl")
include("TaylorN.jl")

include("printing.jl")

end # module


text:
using Dice
using Dice: num_flips, num_nodes, to_dice_ir

code = @dice begin
    function uniform(b::Int)
        a = b/2
        d = true
        bits = Vector(undef, b)
        for i=1:b
            rep = flip(0.5)
            bits[i] = rep
            if i > a
                d &= rep
            else
                d &= true
            end
        end
        return DistInt(bits), d
    end

    a, b = uniform(4)
    CondInt(a, b)
end

# BDD analysis
bdd = compile(code)
infer(code, :bdd)

# IR analysis
# to_dice_ir(code)
# has_dice_binary() && rundice(code)
# has_dice_binary() && infer(code, :ocaml)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    SMPSScenario

Conveniece type that adheres to the [`AbstractScenario`](@ref) abstraction. Obtained when reading scenarios specified in SMPS format.

See also: [`SMPSSampler`](@ref)
"""
struct SMPSScenario{T <: AbstractFloat, M <: AbstractMatrix} <: AbstractScenario
    probability::Probability
    Δq::Vector{T}
    ΔT::M
    ΔW::M
    Δh::Vector{T}
    ΔC::M
    Δd₁::Vector{T}
    Δd₂::Vector{T}
end

function Base.zero(::Type{SMPSScenario{T}}) where T <: AbstractFloat
    return SMPSScenario(Probability(1.0),
                        AdditiveZeroArray{T,1}(),
                        AdditiveZeroArray{T,2}(),
                        AdditiveZeroArray{T,2}(),
                        AdditiveZeroArray{T,1}(),
                        AdditiveZeroArray{T,2}(),
                        AdditiveZeroArray{T,1}(),
                        AdditiveZeroArray{T,1}())
end

function StochasticPrograms.expected(scenarios::Vector{<:SMPSScenario{T}}) where T <: AbstractFloat
    isempty(scenarios) && return zero(SMPSScenario{T})
    expected = reduce(scenarios) do ξ₁, ξ₂
        SMPSScenario(Probability(1.0),
                     probability(ξ₁) * ξ₁.Δq + probability(ξ₂) * ξ₂.Δq,
                     probability(ξ₁) * ξ₁.ΔT + probability(ξ₂) * ξ₂.ΔT,
                     probability(ξ₁) * ξ₁.ΔW + probability(ξ₂) * ξ₂.ΔW,
                     probability(ξ₁) * ξ₁.Δh + probability(ξ₂) * ξ₂.Δh,
                     probability(ξ₁) * ξ₁.ΔC + probability(ξ₂) * ξ₂.ΔC,
                     probability(ξ₁) * ξ₁.Δd₁ + probability(ξ₂) * ξ₂.Δd₁,
                     probability(ξ₁) * ξ₁.Δd₂ + probability(ξ₂) * ξ₂.Δd₂)
    end
    return StochasticPrograms.ExpectedScenario(expected)
end

function StochasticPrograms.scenariotext(io::IO, scenario::SMPSScenario)
    print(io, " and underlying data:\n\n")
    println(io, "Δq = $(scenario.Δq)")
    println(io, "ΔT = $(scenario.ΔT)")
    println(io, "ΔW = $(scenario.ΔW)")
    println(io, "Δh = $(scenario.Δh)")
    println(io, "ΔC = $(scenario.ΔC)")
    println(io, "Δd₁ = $(scenario.Δd₁)")
    print(io, "Δd₂ = $(scenario.Δd₂)")
    return io
end


text:
using Plots
@testset "plotting" begin
    @testset "plotsdf" begin
        @test_throws ArgumentError plotsdf(1.0)
        @test_throws ArgumentError plotsdf(OU,1:2)
        @test_throws ArgumentError plotsdf(1.0,1:2)
        @test_throws ArgumentError plotsdf(OU(1.0,1.0),1)
        @test_throws ArgumentError plotsdf(OU(1.0,1.0),1:2,1)
    end
    @testset "plotasdf" begin
        @test_throws ArgumentError plotasdf(1.0)
        @test_throws ArgumentError plotasdf(OU,1:2,1)
        @test_throws ArgumentError plotasdf(OU(1.0,1.0),1:2)
        @test_throws ArgumentError plotasdf(OU(1.0,1.0),1:2,-1)
        @test_throws ArgumentError plotasdf(OU(1.0,1.0),1:2,-1,1)
    end
    @testset "plotacv" begin
        @test_throws ArgumentError plotacv(1.0)
        @test_throws ArgumentError plotacv(OU,1:2)
        @test_throws ArgumentError plotacv(OUUnknown{1}(1.0,1.0),1:2)
        @test_throws ArgumentError plotacv(OU,10,1)
        @test_throws ArgumentError plotacv(OU(1.0,1.0),1:2,1)
        @test_throws ArgumentError plotacv(OU(1.0,1.0),10,-1)
        @test_throws ArgumentError plotacv(OU(1.0,1.0),-10,1)
        @test_throws ArgumentError plotacv(OU(1.0,1.0),10,1im)
        @test_throws ArgumentError plotacv(OU(1.0,1.0),10,1,1)
    end
    @testset "plotei" begin
        @test_throws ArgumentError plotei(1.0)
        @test_throws ArgumentError plotei(OU,10,1)
        @test_throws ArgumentError plotei(OU(1.0,1.0),10)
        @test_throws ArgumentError plotei(OU(1.0,1.0),1:2,1)
        @test_throws ArgumentError plotei(OU(1.0,1.0),10,-1)
        @test_throws ArgumentError plotei(OU(1.0,1.0),-10,1)
        @test_throws ArgumentError plotei(OU(1.0,1.0),10,1im)
        @test_throws ArgumentError plotei(OU(1.0,1.0),10,1,1)
    end
end


text:
export DelaunayTriangulation

import StaticArrays: SVector, MVector
import .Simplices

include("addnoise.jl")

struct DelaunayTriangulation <: AbstractDelaunayTriangulation
    indices::Vector{Vector{Int32}}

    function DelaunayTriangulation(vertices::Vector{T}; joggle = 1e-8) where T <: Union{AbstractVector, SVector, MVector}
        # Slightly joggle points to avoid problems with QHull
        if joggle > 0
            addnoise!(vertices; joggle_factor = joggle)
        end
        simplex_inds = Simplices.Delaunay.delaunay(vertices)
        new(simplex_inds)
    end

    DelaunayTriangulation(vertices::Dataset; joggle = 1e-8) = DelaunayTriangulation(vertices.data)
end


text:
using DEBTrait, Test

el = "C4H9NO2"
chemFormBiom = [1, 1.8, 0.2, 0.5, 0, 0, 0]

chemical_composition = DEBTrait.ThermoStoichWizard.extract_composition(el) # CHNOSP
@test chemical_composition == [4,9,1,2,0,0]

stoich_electron_donor = DEBTrait.ThermoStoichWizard.get_stoich_electron_donor(el)
@test stoich_electron_donor == [-1, -10, 4, 1, 0, 0, 21, 18, 0, 0]

stoich_electron_acceptor = DEBTrait.ThermoStoichWizard.get_stoich_electron_acceptor()
@test stoich_electron_acceptor == [0., 2., 0., 0., 0., 0., -4., -4., -1., 0.]

stoich_cat_rxns = DEBTrait.ThermoStoichWizard.get_stoich_catabolic_reaction(el)
@test stoich_cat_rxns == [-1., -1., 4., 1., 0., 0., 3., 0., -4.5, 0.]

stoich_anabolic_O2, stoich_anabolic_HCO3 = DEBTrait.ThermoStoichWizard.get_stoich_anabolic_reaction(el, chemFormBiom)
@test stoich_anabolic_O2  ≈ [-0.25, 0.15, 0.0,  0.05,  0., 0., -0.05, 0.,  -0.075, 1.0] atol =1e-4
@test stoich_anabolic_HCO3  ≈ [-0.23333333, 0.16666667, -0.06666667,  0.03333333,  0., 0., -0.1, 0.,  0., 1.] atol = 1e-4

out = DEBTrait.ThermoStoichWizard.get_lambda(el, chemFormBiom)


text:
using Test
using PerronFrobenius

pts = rand(3, 10000)
μF = [4, 4, 4]

μ_avg = PerronFrobenius.average_measure(pts, μF , [μF, μF])
μ_induced = PerronFrobenius.induced_measure(pts, μF, μF)

@testset "Average measure" begin
	@testset "Average measure sums to 1" begin
		@test sum(μ_avg) ≈ 1
	end

	@testset "Measure induced from multiple same resolutions matches original" begin
		μ_induced = induced_measure(pts, μF, μF)
		@test all(μ_avg .≈ μ_induced)
	end
end


text:
a = rand(10^7) # Vector 1D uniforme en [0,1)

sum(a)

#Pkg.add("BenchmarkTools")

using BenchmarkTools

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # Haz un directorio temporario

# compila a una biblioteca compartida pipeando C_code a gcc
# (funciona sólo con gcc instalado):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

# define una funcion de Julia que llama a la función de C:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

c_sum(a)

c_sum(a) ≈ sum(a) # teclea \approx y luego <TAB> para obtener el símbolo ≈

≈  # alias para la función `isapprox`

?isapprox

c_bench = @benchmark c_sum($a)

println("C: Tiempo más rápido fue $(minimum(c_bench.times) / 1e6) msec")

d = Dict()  # un diccionario a.k.a un arreglo asociativo
d["C"] = minimum(c_bench.times) / 1e6  # en milisegundos
d

using Plots
gr()

t = c_bench.times / 1e6 # tiempos en milisegundos
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")

#Pkg.add("PyCall")

using PyCall

# Llama una función de bajo nivel de PyCall para obtener una lista de Python
# porque por default PyCall convertirá a un arreglo de NumPy en vez (benchmarkeamos NumPy más abajo)

apy_list = PyCall.array2py(a, 1, 1)

# obtener el "sum" interno de Python
pysum = pybuiltin("sum")

pysum(a)

pysum(a) ≈ sum(a)

py_list_bench = @benchmark $pysum($apy_list)

d["Python interno"] = minimum(py_list_bench.times) / 1e6
d

using Conda
#Conda.add("numpy")

numpy_sum = pyimport("numpy")["sum"]
apy_numpy = PyObject(a) # convierte a un arreglo de NumPy por default

py_numpy_bench = @benchmark $numpy_sum($apy_numpy)

numpy_sum(apy_list) # python thing

numpy_sum(apy_list) ≈ sum(a)

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

py"""
def py_sum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""

sum_py = py"py_sum"

py_hand = @benchmark $sum_py($apy_list)

sum_py(apy_list)

sum_py(apy_list) ≈ sum(a)

d["Python hecho a mano"] = minimum(py_hand.times) / 1e6
d

@which sum(a)

j_bench = @benchmark sum($a)

d["Julia interno"] = minimum(j_bench.times) / 1e6
d

function mysum(A)
    s = 0.0  # s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end

j_bench_hand = @benchmark mysum($a)

d["Julia hecho a mano"] = minimum(j_bench_hand.times) / 1e6
d

for (key, value) in sort(collect(d))
    println(rpad(key, 20, "."), lpad(round(value, 1), 8, "."))
end

for (key, value) in sort(collect(d), by=x->x[2])
    println(rpad(key, 20, "."), lpad(round(value, 2), 10, "."))
end

max_stars_repo_path: es-es/intro-to-julia-ES/plutos/12. Julia es rapido.jl

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
export asarray

"""
    asarray(x::Number[, parent::AbstractArray]) -> AbstactArray

Return a 0-dimensional array with item `x`, otherwise, do nothing.
If a `parent` is supplied, it will try to match the parent array type.
"""
asarray(x::Number) = fill(x, ())
asarray(x::Number, arr::Array) = fill(x, ())
asarray(x::AbstractArray, args...) = x

tsetdiff(t::Tuple, b) = setdiff!(collect(t), b)
tunique(t::Tuple) = unique!(collect(t))

"""
    nopermute(ix,iy)

check that all values in `iy` that are also in `ix` have the same relative order,

# example

```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: nopermute

julia> nopermute((1,2,3),(1,2))
true

julia> nopermute((1,2,3),(2,1))
false
```
e.g. `nopermute((1,2,3),(1,2))` is true while `nopermute((1,2,3),(2,1))` is false
"""
function nopermute(ix::NTuple, iy::NTuple)
    i, j, jold = 1, 1, 0
    # find each element of iy in ix and check that the order is the same
    for i in 1:length(iy)
        j = findfirst(==(iy[i]), ix)
        (j === nothing || j <= jold) && return false
        jold = j
    end
    return true
end

"""
    allunique(ix::Tuple)

return true if all elements of `ix` appear only once in `ix`.

# example

```jldoctest; setup = :(using OMEinsum)
julia> using OMEinsum: allunique

julia> allunique((1,2,3,4))
true

julia> allunique((1,2,3,1))
false
```
"""
allunique(ix::NTuple) = all(i -> count(==(i), ix) == 1, ix)

function conditioned_permutedims(A::AbstractArray{T,N}, perm, ind=()) where {T,N}
    if any(i-> (@inbounds perm[i]!=i), 1:N)
        @debug "conditioned_permutedims" size(A) Tuple(perm) Tuple(ind)
        return tensorpermute(A, perm)
    else
        return A
    end
end

function align_eltypes(xs::AbstractArray...)
    T = promote_type(eltype.(xs)...)
    return map(x->eltype(x)==T ? x : T.(x), xs)
end

function align_eltypes(xs::AbstractArray{T}...) where T
    xs
end

"""
    tensorpermute(A, perm)

Like `permutedims(A, perm)`, but calls the faster `TensorOperations.tensorcopy` when possible.
"""
function tensorpermute(A::StridedArray{T,N}, perm) where {T,N}
    TensorOperations.tensorcopy(A, ntuple(identity,N), perm)
end
tensorpermute(A::AbstractArray, perm) = permutedims(A, perm)


text:
# const LKH_VERSION = "LKH-2.0.9"
# const LKH_SRC_URL = "http://webhotel4.ruc.dk/~keld/research/LKH/$(LKH_VERSION).tgz"
# const LKH_WIN_EXE_URL = "http://webhotel4.ruc.dk/~keld/research/LKH/LKH-2.exe"

const LKH_VERSION = "LKH-3.0.6"
const LKH_SRC_URL = "http://webhotel4.ruc.dk/~keld/research/LKH-3/$(LKH_VERSION).tgz"
const LKH_WIN_EXE_URL = "http://webhotel4.ruc.dk/~keld/research/LKH-3/LKH-3.exe"

function download_win_exe()
    lkh_exe = joinpath(@__DIR__, "LKH.exe")
    download(LKH_WIN_EXE_URL, lkh_exe)
    return lkh_exe
end

function build_LKH()
    src_tarball = joinpath(@__DIR__, "lkh_src.tgz")
    download(LKH_SRC_URL, src_tarball)
    src_dir = joinpath(@__DIR__, "src")
    run(`tar zxvf $(src_tarball)`)
    cd(LKH_VERSION)
    run(`make`)
    mv("LKH", joinpath(@__DIR__, "LKH"), force=true)

    lkh_exe = joinpath(@__DIR__, "LKH")
    return lkh_exe
end

function download_LKH()
    if Sys.islinux() && Sys.ARCH == :x86_64
        return build_LKH()
    elseif Sys.isapple()
        return build_LKH()
    elseif Sys.iswindows()
        return download_win_exe()
    end
    error("Unsupported operating system. Only 64-bit linux, macOS, and Windows are supported.")
end

function install_LKH()
    executable = get(ENV, "LKH_EXECUTABLE", nothing)
    if !haskey(ENV, "LKH_EXECUTABLE")
        executable = download_LKH()
        ENV["LKH_EXECUTABLE"] = executable
    end

    if executable === nothing
        error("Environment variable `LKH_EXECUTABLE` not found.")
    else
        # pr2392 = joinpath(@__DIR__, LKH_VERSION, "pr2392.par")
        # run(`$(executable) $(pr2392)`)
    end

    open(joinpath(@__DIR__, "deps.jl"), "w") do io
        write(io, "const LKH_EXECUTABLE = \"$(escape_string(executable))\"\n")
    end
end

install_LKH()


text:
"""
# Right
In access control systems, a right is an action or operation (read, modify, hunt, fish, ...) that may be performed on a target object.

## Fields
- `label::string`: A free textual label or name to help identify the object.
"""
mutable struct Right <: AbstractRight
    label::String;
end;
function Right(;n::Int = 1)::Vector{Right}
    useless(x) = Right("a right");
    broadcast(useless, repeat([0], n));
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
"""
    zs = smoothconv(z,nas)

Smoothen field `z(ns,ns)` with a circular kernel of diameter `ns/nas` using convolution.
Takes into account missing values.
"""
	function smoothconv(zi,nas)

        @compat iinan=findall(isnan.(zi))
        @compat iinotnan=findall(.~isnan.(zi))
        zi[iinan].=0.

        nss=size(zi);
        ns=nss[1];
        sdim=div(ns,nas)/2; # the smoothing sigma is half large scale pixel wide

        mask=zeros(ns,ns);
        for i=1:ns
           for j=1:ns
               kx=i-1;
               ky=j-1;
               if(i>ns/2+1)
                  kx=i-ns-1 ;
               end
               if(j>ns/2+1)
                  ky=j-ns-1 ;
               end
               r2=kx*kx+ky*ky;
               mask[i,j]=exp(-(r2/(sdim*sdim))/2);
           end
        end
        fm=fft(mask)
        zf=real(ifft(fm.*fft(zi)))/sum(mask)
        if length(iinan)>0
           zi1=deepcopy(zi)
           zi1[iinotnan].=1.0
           zf=zf./(real(ifft(fm.*fft(zi1)))/sum(mask))
        end
	zf[iinan].=NaN
        return zf
   end


text:
using Functors, Test

struct Foo
  x
  y
end
@functor Foo

struct Bar
  x
end
@functor Bar

struct Baz
  x
  y
  z
end
@functor Baz (y,)

struct NoChildren
  x
  y
end

@static if VERSION >= v"1.6"
  @testset "ComposedFunction" begin
    f1 = Foo(1.1, 2.2)
    f2 = Bar(3.3)
    @test Functors.functor(f1 ∘ f2)[1] == (outer = f1, inner = f2)
    @test Functors.functor(f1 ∘ f2)[2]((outer = f1, inner = f2)) == f1 ∘ f2
    @test fmap(x -> x + 10, f1 ∘ f2) == Foo(11.1, 12.2) ∘ Bar(13.3)
  end
end

@testset "Nested" begin
  model = Bar(Foo(1, [1, 2, 3]))

  model′ = fmap(float, model)

  @test model.x.y == model′.x.y
  @test model′.x.y isa Vector{Float64}
end

@testset "Exclude" begin
  f(x::AbstractArray) = x
  f(x::Char) = 'z'

  x = ['a', 'b', 'c']
  @test fmap(f, x)  == ['z', 'z', 'z']
  @test fmap(f, x; exclude = x -> x isa AbstractArray) == x

  x = (['a', 'b', 'c'], ['d', 'e', 'f'])
  @test fmap(f, x)  == (['z', 'z', 'z'], ['z', 'z', 'z'])
  @test fmap(f, x; exclude = x -> x isa AbstractArray) == x
end

@testset "Walk" begin
  model = Foo((0, Bar([1, 2, 3])), [4, 5])

  model′ = fmapstructure(identity, model)
  @test model′ == (; x=(0, (; x=[1, 2, 3])), y=[4, 5])
end

@testset "Property list" begin
  model = Baz(1, 2, 3)
  model′ = fmap(x -> 2x, model)

  @test (model′.x, model′.y, model′.z) == (1, 4, 3)
end

@testset "fcollect" begin
  m1 = [1, 2, 3]
  m2 = 1
  m3 = Foo(m1, m2)
  m4 = Bar(m3)
  @test all(fcollect(m4) .=== [m4, m3, m1, m2])
  @test all(fcollect(m4, exclude = x -> x isa Array) .=== [m4, m3, m2])
  @test all(fcollect(m4, exclude = x -> x isa Foo) .=== [m4])

  m1 = [1, 2, 3]
  m2 = Bar(m1)
  m0 = NoChildren(:a, :b)
  m3 = Foo(m2, m0)
  m4 = Bar(m3)
  @test all(fcollect(m4) .=== [m4, m3, m2, m1, m0])
end

struct FFoo
  x
  y
  p
end
@flexiblefunctor FFoo p

struct FBar
  x
  p
end
@flexiblefunctor FBar p

struct FBaz
  x
  y
  z
  p
end
@flexiblefunctor FBaz p

@testset "Flexible Nested" begin
  model = FBar(FFoo(1, [1, 2, 3], (:y, )), (:x,))

  model′ = fmap(float, model)

  @test model.x.y == model′.x.y
  @test model′.x.y isa Vector{Float64}
end

@testset "Flexible Walk" begin
  model = FFoo((0, FBar([1, 2, 3], (:x,))), [4, 5], (:x, :y))

  model′ = fmapstructure(identity, model)
  @test model′ == (; x=(0, (; x=[1, 2, 3])), y=[4, 5])

  model2 = FFoo((0, FBar([1, 2, 3], (:x,))), [4, 5], (:x,))

  model2′ = fmapstructure(identity, model2)
  @test model2′ == (; x=(0, (; x=[1, 2, 3])))
end

@testset "Flexible Property list" begin
  model = FBaz(1, 2, 3, (:x, :z))
  model′ = fmap(x -> 2x, model)

  @test (model′.x, model′.y, model′.z) == (2, 2, 6)
end

@testset "Flexible fcollect" begin
  m1 = 1
  m2 = [1, 2, 3]
  m3 = FFoo(m1, m2, (:y, ))
  m4 = FBar(m3, (:x,))
  @test all(fcollect(m4) .=== [m4, m3, m2])
  @test all(fcollect(m4, exclude = x -> x isa Array) .=== [m4, m3])
  @test all(fcollect(m4, exclude = x -> x isa FFoo) .=== [m4])

  m0 = NoChildren(:a, :b)
  m1 = [1, 2, 3]
  m2 = FBar(m1, ())
  m3 = FFoo(m2, m0, (:x, :y,))
  m4 = FBar(m3, (:x,))
  @test all(fcollect(m4) .=== [m4, m3, m2, m0])
end


text:
module Fields
    abstract Field

    immutable Q <: Field end
    immutable V <: Field end
    immutable P <: Field end
    immutable F <: Field end
    immutable U <: Field end
    immutable R <: Field end
    immutable Λ <: Field end
    immutable Φ <: Field end
end


text:
module Types

    using RandomExtensions
    import Dates

    export Entity, EntityData, GAmodel, targetParams, GAParams

    abstract type Entity end
    # -------

    mutable struct EntityData
        entity
        generation::Int
        scores :: Array

        EntityData(entity, generation::Int) = new(entity, generation)
        EntityData(entity, model) = new(entity, model.params.currentGeneration)
    end

    mutable struct GAParams
        populationSize ::  Int
        generations ::  Int
        genomeSize ::  Int
        maxGenomeSize :: Int
        crossoverRate :: Float64
        mutationRate :: Float64
        elitism  :: Bool # Keep previous generation's fittest individual in place of worst in current
        historyPath :: IOStream  # Path to save log of fitness history at each generation. Can be used to plot on Excel chart etc.

        totalFitness :: Float64
        targetFitness :: Float64
        targetFitnessCount ::  Int
        currentGeneration ::  Int

        #thisGeneration :: Array
        #nextGeneration :: Array
        #fitnessTable :: Array{Float64,1}
    end

    # -------

    mutable struct GAmodel
        initial_pop_size::Int
        gen_num::Int

        population::Array
        #pop_data::Array{EntityData}
        #freezer::Array{EntityData}

        rng::AbstractRNG

        ga
        specific_fitness

        scores :: Array
        params :: GAParams

        instructionsSet #:: AbstractDict

        GAmodel() = new(0, 1, Any[], MersenneTwister(time_ns()), nothing, nothing, [],
                GAParams(136, 1000 , 150, 150, 0.7, 0.01, true, open("output.txt", "w") ,  0.0 , 0, 0.0 , 0 ),
                Dict() )
        GAmodel(params :: GAParams) = new(0, 1, Any[], MersenneTwister(time_ns()),
                nothing, nothing, [], params, Dict())
    end

    struct Function
        ## <summary>
        ## Controls how functions read input (,) from parent memory: either at the current memory data pointer or from the start of memory.
        ## If true, input will be read from position 0 from the parent. Meaning, the first input value that the parent read will be the first input value the function gets, regardless of the parent's current memory data position. This may make it easier for the GA to run the function, since it does not require an exact memory position before calling the function.
        ## If false (default), input will be read from the current memory data position of the parent. Meaning, if the parent has shifted the memory pointer up 3 slots, the function will begin reading from memory at position 3.
        ## </summary>
        ReadInputAtMemoryStart :: Bool
        ## <summary>
        ## Custom max iteration counts for functions.
        ## </summary>
        MaxIterationCount :: Int
    end

    struct FunctionInst #:: Function
        ## <summary>
        ## Starting instruction index for this function, within the program code.
        ## </summary>
         InstructionPointer
    end

    struct GAStatus
        ## <summary>
        ## Best fitness so far.
        ## </summary>
        Fitness :: Float64
        ## <summary>
        ## Best true fitness so far, used to determine when a solution is found.
        ## </summary>
        TrueFitness :: Float64
        ## <summary>
        ## Best program so far.
        ## </summary>
        Program :: String
        ## <summary>
        ## Best program output so far.
        ## </summary>
        Output :: String
        ## <summary>
        ## Current iteration (generation) count.
        ## </summary>
        Iteration :: Int
        ## <summary>
        ## Count of status prompts.
        ## </summary>
        StatusCount :: Int
        ## <summary>
        ## Number of instructions executed by the best program.
        ## </summary>
        Ticks :: Int
        ## <summary>
        ## Number of instructions executed by the best program, including function calls.
        ## </summary>
        TotalTicks :: Int
        ## <summary>
        ## Time of last improved evolution.
        ## </summary>
        LastChangeDate :: Dates.DateTime

        GAStatus() = new(0.0, 0.0, "", "", 0, 0, 0, 0, Dates.now())
    end

    #=struct targetParams
        targetString :: String
        targetFitness :: String
    end =#

end  # module Types


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 = "FFMPEG"
version = v"4.4"

# Collection of sources required to build FFMPEG
sources = [
    ArchiveSource("https://ffmpeg.org/releases/ffmpeg-$(version.major).$(version.minor).tar.xz",
                  "06b10a183ce5371f915c6bb15b7b1fffbe046e8275099c96affc29e17645d909"),
]

# Bash recipe for building across all platforms
# TODO: Theora once it's available
function script(; ffplay=false)
    "FFPLAY=$(ffplay)\n" * raw"""
cd $WORKSPACE/srcdir
cd ffmpeg-*/
sed -i 's/-lflite"/-lflite -lasound"/' configure
apk add coreutils yasm

if [[ "${target}" == *-linux-* ]]; then
    export ccOS="linux"
elif [[ "${target}" == *-apple-* ]]; then
    export ccOS="darwin"
elif [[ "${target}" == *-w32-* ]]; then
    export ccOS="mingw32"
elif [[ "${target}" == *-w64-* ]]; then
    export ccOS="mingw64"
elif [[ "${target}" == *-unknown-freebsd* ]]; then
    export ccOS="freebsd"
else
    export ccOS="linux"
fi

if [[ "${target}" == x86_64-* ]]; then
    export ccARCH="x86_64"
elif [[ "${target}" == i686-* ]]; then
    export ccARCH="i686"
elif [[ "${target}" == arm-* ]]; then
    export ccARCH="arm"
elif [[ "${target}" == aarch64-* ]]; then
    export ccARCH="aarch64"
elif [[ "${target}" == powerpc64le-* ]]; then
    export ccARCH="powerpc64le"
else
    export ccARCH="x86_64"
fi

if [[ "${target}" == arm-* ]]; then
    export CUDA_ARGS=""
elif [[ "${target}" == *-apple-* ]]; then
    export CUDA_ARGS=""
elif [[ "${target}" == *-unknown-freebsd* ]]; then
    export CUDA_ARGS=""
else
    export CUDA_ARGS="--enable-nvenc --enable-cuda-llvm"
fi

EXTRA_FLAGS=()
if [[ "${target}" == *-darwin* ]]; then
    EXTRA_FLAGS+=(--objcc="${CC} -x objective-c")
fi
if [[ "${FFPLAY}" == "true" ]]; then
    EXTRA_FLAGS+=("--enable-ffplay")
fi

# Remove `-march` flags
sed -i 's/cpuflags="-march=$cpu"/cpuflags=""/g' configure

./configure            \
  --enable-cross-compile \
  --cross-prefix=/opt/${target}/bin/${target}- \
  --arch=${ccARCH}     \
  --target-os=${ccOS}  \
  --cc="${CC}"         \
  --cxx="${CXX}"       \
  --dep-cc="${CC}"     \
  --ar=ar              \
  --nm=nm              \
  --sysinclude=${prefix}/include \
  --pkg-config=$(which pkg-config) \
  --pkg-config-flags=--static \
  --prefix=$prefix     \
  --sysroot=/opt/${target}/${target}/sys-root \
  --extra-libs=-lpthread \
  --enable-gpl         \
  --enable-version3    \
  --enable-nonfree     \
  --disable-static     \
  --enable-shared      \
  --enable-pic         \
  --disable-debug      \
  --disable-doc        \
  --enable-avresample  \
  --enable-libass      \
  --enable-libfdk-aac  \
  --enable-libfreetype \
  --enable-libmp3lame  \
  --enable-libopus     \
  --enable-libvorbis   \
  --enable-libx264     \
  --enable-libx265     \
  --enable-libvpx      \
  --enable-encoders    \
  --enable-decoders    \
  --enable-muxers      \
  --enable-demuxers    \
  --enable-parsers     \
  --enable-openssl     \
  --disable-schannel   \
  --extra-cflags="-I${prefix}/include" \
  --extra-ldflags="-L${libdir}" ${CUDA_ARGS} \
  "${EXTRA_FLAGS[@]}"
make -j${nproc}
if [[ "${FFPLAY}" == "true" ]]; then
    # Manually install only the FFplay binary
    mv "ffplay${exeext}" "${bindir}/ffplay${exeext}"
else
    # Install all FFMPEG stuff: libraries, executables, header files, etc...
    make install
fi
install_license LICENSE.md COPYING.*
"""
end

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

preferred_gcc_version = v"8"

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Symbolics
import Symbolics: getsource, getdefaultval, wrap, unwrap, getname
import SymbolicUtils: Term, symtype, FnType
using Test

@variables t
Symbolics.@register fff(t)
@test isequal(fff(t), Symbolics.Num(Symbolics.Term{Real}(fff, [Symbolics.value(t)])))

## @variables

many_vars = @variables t=0 a=1 x[1:4]=2 y[1:4](t)=3 w[1:4] = 1:4 z[1:4](t) = 2:5 p[1:4](..)

@test all(t->getsource(t)[1] === :variables, many_vars)
@test getdefaultval(t) == 0
@test getdefaultval(a) == 1
@test_throws ErrorException getdefaultval(x)
@test getdefaultval(x[1]) == 2
@test getdefaultval(y[2]) == 3
@test getdefaultval(w[2]) == 2
@test getdefaultval(w[4]) == 4
@test getdefaultval(z[3]) == 4

nxt = Namespace(unwrap(x), unwrap(t))
nxt_1 = Namespace(nxt, unwrap(x))
@test getname(nxt) == Symbol(:x, :(.), :t)
@test getname(nxt_1) == Symbol(:x, :(.), :t, :(.), :x)

@test p[1] isa Term
@test symtype(p[1]) <: FnType{Tuple, Real}
@test p[1](t) isa Symbolics.Num

## Wrapper types

abstract type AbstractFoo{T} end

struct Foo{T}<:AbstractFoo{T}
end
@symbolic_wrap struct FooWrap{T} <: AbstractFoo{T}
    val::Foo{T}
end

Symbolics.unwrap(r::FooWrap) = r.val

@test Symbolics.wrapper_type(AbstractFoo) == FooWrap
@test Symbolics.wrapper_type(AbstractFoo{Int}) == FooWrap # removes the type param.
                                                          # If required,
                                                          # a new method must be defined
                                                          # by the user of @symbolic_wrap

x = Foo{Int}()
@test wrap(x) isa FooWrap{Int}
@test unwrap(wrap(x)) isa Foo{Int}

@wrapped function foo(f::AbstractFoo, x::Real)
    x > 1 ? (x > 5 ? "hi" : Foo{Int}()) : 0
end

@test !applicable(foo, x, 2)
@test applicable(foo, x, wrap(2))
@test applicable(foo, wrap(x), 2)
@test applicable(foo, wrap(x), wrap(2))

@test foo(x, wrap(2)) isa FooWrap
@test foo(x, wrap(1)) isa Num
@test foo(x, wrap(6)) isa String


text:
module LambdaCalculus
    import Base: show

    export
        BOOL,
        TRUE,
        FALSE,
        And,
        Or,
        Not,
        Eq,
        Xor,

        ZERO,
        SUC,
        Add,
        Sub

    include("bool.jl")
    include("operators.jl")
    include("numbers.jl")

end # module


text:
# Define the nearest neighbor term `S⋅S` for the Heisenberg model
using ITensors
using ITensorMPOCompression
using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%1.3f", f)

function ITensors.op(::OpName"expS⋅S", ::SiteType"S=1/2",
    s1::Index, s2::Index; τ::Number)
    O = 0.5 * op("S+", s1) * op("S-", s2) +
        0.5 * op("S-", s1) * op("S+", s2) +
              op("Sz", s1) * op("Sz", s2)
    return exp(τ * O)
end

function fix_link_tags!(W::ITensor)
    is=filterinds(inds(W),tags="Site")[1]
    ts2=String(tags(is)[3])
    @assert ts2[1:2]=="n="
    nsite::Int64=tryparse(Int64,ts2[3:end])
    ils=filterinds(inds(W),"Link")
    n1::Int64=tryparse(Int64,String(tags(ils[1])[2])[3:end])
    if length(ils)>1
        n2::Int64=tryparse(Int64,String(tags(ils[2])[2])[3:end])
    else
        n2=0
    end
    #@show nsite n1 n2
    for i in 1:length(ils)
        tn=String(tags(ils[i])[2])
        if tn[1]=='n'
            if i==1
                if n2==nsite
                    tl="l=$(nsite-1)"
                end
                if n2==nsite-1
                    tl="l=$nsite"
                end
                if n2==0
                    if n1==nsite
                        tl="l=$nsite"
                    else
                        tl="l=$(nsite-1)"
                    end
                end
            end
            if i==2
                if n1==nsite
                    tl="l=$(nsite-1)"
                end
                if n1==nsite-1
                    tl="l=$nsite"
                end

            end
            replacetags!(W,tn,tl,tags="Link")
        end
    end
end
function fix_link_tags!(H::MPO)
    for W in H
        fix_link_tags!(W)
    end
end

N=10
τ = -0.01im
even = [("expS⋅S", (n, n+1), (τ = τ,)) for n in 1:2:N-1]
odd  = [("expS⋅S", (n, n+1), (τ = τ,)) for n in 2:2:N-1]
s = siteinds("S=1/2", N)
even_expτH = ops(even, s)
odd__expτH = ops(odd , s)

H=make_Heisenberg_AutoMPO(s,N)
orthogonalize!(H,orth=right)
pprint(H)
He = apply(even_expτH, H)
fix_link_tags!(He)
orthogonalize!(He,orth=right)
get_Dw(He)
pprint(He)
# W=He[1]
# pprint(W)
# d,n,r,c,lnr,lnc=parse_links(W)
# for j in 1:dim(c)
#     @show slice(W,c=>j)
# end
# # Ho = apply(odd__expτH, He)
# # pprint(Ho)
# Nothing


text:
const ST = ScientificTypes

ST.scitype(::PersistenceDiagram, ::ST.DefaultConvention; kwargs...) = PersistenceDiagram


text:
"""
for expanding macros
"""
module MLBootstrap
    using MLStyle
    using MLStyle.Render
    using MLStyle.Modules.AST
end
include("CodeDumper.jl")
using MLStyle
export bootstrap_gen_pipeline

struct Discard end
const discard = Discard()

function conditional_macro_expand(enable_macros)
    inner = @λ begin
        (Expr(:macrocall, name && if name in enable_macros end, args...) && e) -> inner(macroexpand(MLBootstrap, e))
        e :: Expr           -> Expr(e.head, map(inner, e.args)...)
        a                   -> a
    end
end

function bootstrap(expr, macros_to_expand::Vector{Any} = convert(Vector{Any}, [:match, :matchast]))
    (modname, toplevels) = @match expr begin
        :(module $modname; $(toplevels...) end) => (modname, toplevels)
    end

    macros_to_expand = [Symbol("@", each) for each in macros_to_expand]
    io = IOBuffer()
    println(io, "# This file is automatically generated by MLStyle Bootstrap Tools.")
    println(io, "module $(modname)")
    println(io, "using MLStyle")
    expr_to_str = conditional_macro_expand(macros_to_expand)
    for toplevel in toplevels
        print(io, "@eval \$(")
        pprint(io, expr_to_str(toplevel))
        println(io, ")")
    end
    println(io, "end")
    String(take!(io))
end

gen_file(fname, node) = open(fname, "w") do f
        write(f, bootstrap(node))
end


text:
struct FastAIDataset
    name
    subfolder
    extension
    description
    checksum
    datadepname
    size
end

const ROOT_URL = "https://s3.amazonaws.com/fast-ai-"

function FastAIDataset(
        name, subfolder, checksum = ""; extension = "tgz", description = "", datadepname = name, size = "???")
    return FastAIDataset(name, subfolder, extension, description, checksum, datadepname, size)
end

const DATASETCONFIGS = [
    # imageclas
    FastAIDataset("CUB_200_2011", "imageclas"),
    FastAIDataset("bedroom", "imageclas"),
    FastAIDataset("caltech_101", "imageclas"),
    FastAIDataset("cifar10", "imageclas", "637c5814e11aefcb6ee76d5f59c67ddc8de7f5b5077502a195b0833d1e3e4441"),
    FastAIDataset("cifar100", "imageclas", "085ac613ceb0b3659c8072143ae553d5dd146b3c4206c3672a56ed02d0e77d28"),
    FastAIDataset("food-101", "imageclas"),
    FastAIDataset("imagenette-160", "imageclas"),
    FastAIDataset("imagenette-320", "imageclas"),
    FastAIDataset("imagenette", "imageclas"),
    FastAIDataset("imagenette2-160", "imageclas", "64d0c4859f35a461889e0147755a999a48b49bf38a7e0f9bd27003f10db02fe5"),
    FastAIDataset("imagenette2-320", "imageclas", "569b4497c98db6dd29f335d1f109cf315fe127053cedf69010d047f0188e158c"),
    FastAIDataset("imagenette2", "imageclas"),
    FastAIDataset("imagewang-160", "imageclas"),
    FastAIDataset("imagewang-320", "imageclas"),
    FastAIDataset("imagewang", "imageclas"),
    FastAIDataset("imagewoof-160", "imageclas"),
    FastAIDataset("imagewoof-320", "imageclas"),
    FastAIDataset("imagewoof", "imageclas"),
    FastAIDataset("imagewoof2-160", "imageclas", "663c22f69c2802d85e2a67103c017e047096702ffddf9149a14011b7002539bf"),
    FastAIDataset("imagewoof2-320", "imageclas"),
    FastAIDataset("imagewoof2", "imageclas"),
    FastAIDataset("mnist_png", "imageclas", "9e18edaa3a08b065d8f80a019ca04329e6d9b3e391363414a9bd1ada30563672"),
    FastAIDataset("mnist_var_size_tiny", "imageclas", "8a0f6ca04c2d31810dc08e739c7fa9b612e236383f70dd9fc6e5a62e672e2283"),
    FastAIDataset("oxford-102-flowers", "imageclas"),
    FastAIDataset("oxford-iiit-pet", "imageclas"),
    FastAIDataset("stanford-cars", "imageclas"),

    # nlp
    FastAIDataset("ag_news_csv", "nlp"),
    FastAIDataset("amazon_review_full_csv", "nlp"),
    FastAIDataset("amazon_review_polarity_csv", "nlp"),
    FastAIDataset("dbpedia_csv", "nlp"),
    FastAIDataset("giga-fren", "nlp"),
    FastAIDataset("imdb", "nlp"),
    FastAIDataset("sogou_news_csv", "nlp"),
    FastAIDataset("wikitext-103", "nlp"),
    FastAIDataset("wikitext-2", "nlp"),
    FastAIDataset("yahoo_answers_csv", "nlp"),
    FastAIDataset("yelp_review_full_csv", "nlp"),
    FastAIDataset("yelp_review_polarity_csv", "nlp"),

    # imagelocal
    FastAIDataset("biwi_head_pose", "imagelocal"),
    FastAIDataset("camvid", "imagelocal"),
    FastAIDataset("pascal-voc", "imagelocal"),
    FastAIDataset("pascal_2007", "imagelocal"),
    FastAIDataset("pascal_2012", "imagelocal"),
    FastAIDataset("siim_small", "imagelocal"),
    FastAIDataset("skin-lesion", "imagelocal"),
    FastAIDataset("tcga-small", "imagelocal"),

    # sample
    FastAIDataset("adult_sample", "sample"),
    FastAIDataset("biwi_sample", "sample"),
    FastAIDataset("camvid_tiny", "sample", "cd42a9bdd8ad3e0ce87179749beae05b4beb1ae6ab665841180b1d8022fc230b"),
    FastAIDataset("dogscats", "sample"),
    FastAIDataset("human_numbers", "sample"),
    FastAIDataset("imdb_sample", "sample"),
    FastAIDataset("mnist_sample", "sample"),
    FastAIDataset("mnist_tiny", "sample"),
    FastAIDataset("movie_lens_sample", "sample"),
    FastAIDataset("planet_sample", "sample"),
    FastAIDataset("planet_tiny", "sample"),

    # coco
    FastAIDataset("coco_sample", "coco", "56960c0ac09ff35cd8588823d37e1ed0954cb88b8bfbd214a7763e72f982911c", size = "3GB"),
    FastAIDataset("train2017", "coco", datadepname="coco-train2017", extension="zip"),
    FastAIDataset("val2017", "coco", datadepname="coco-val2017", extension="zip"),
    FastAIDataset("test2017", "coco", datadepname="coco-test2017", extension="zip"),
    FastAIDataset("unlabeled2017", "coco", datadepname="coco-unlabeled2017", extension="zip"),
    FastAIDataset("image_info_test2017", "coco", datadepname="coco-image_info_test2017", extension="zip"),
    FastAIDataset("image_info_unlabeled2017", "coco", datadepname="coco-image_info_unlabeled2017", extension="zip"),
    FastAIDataset("annotations_trainval2017", "coco", datadepname="coco-annotations_trainval2017", extension="zip"),
    FastAIDataset("stuff_annotations_trainval2017", "coco", datadepname="coco-stuff_annotations_trainval2017", extension="zip"),
    FastAIDataset("panoptic_annotations_trainval2017", "coco", datadepname="coco-panoptic_annotations_trainval2017", extension="zip"),
]

const DATASETS = [d.datadepname for d in DATASETCONFIGS]
const DATASETS_IMAGECLASSIFICATION = vcat(
    [d.datadepname for d in DATASETCONFIGS if d.subfolder == "imageclas"],
    ["mnist_sample", "mnist_tiny", "dogscats"],

)

function DataDeps.DataDep(d::FastAIDataset)
    return DataDep(
        "fastai-$(d.datadepname)",
        """
        "$(d.name)" from the fastai dataset repository (https://course.fast.ai/datasets)

        $(d.description)

        Download size: $(d.size)
        """,
        "$(ROOT_URL)$(d.subfolder)/$(d.name).$(d.extension)",
        d.checksum,
        post_fetch_method = DataDeps.unpack,
    )
end

function initdatadeps()
    for d in DATASETCONFIGS
        DataDeps.register(DataDep(d))
    end
end


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

using Test, SimilaritySearch, InvertedFiles, LinearAlgebra, UnicodePlots

function runtest(; dim, n, m,
    numcenters=5ceil(Int, sqrt(n)), k=10, centersrecall=0.95, kbuild=1, ksearch=1,
    parallel_block=1000, ordering=DistanceOrdering(), minrecall=1.0, initial=:dnet, maxiters=0)
    A = randn(Float32, dim, n)
    B = randn(Float32, dim, m)
    X = MatrixDatabase(A)
    Q = MatrixDatabase(B)
    dist = SqL2Distance()
    seq = ExhaustiveSearch(dist, X)
    @info "creating gold standard"
    gsearchtime = @elapsed Igold, Dgold = searchbatch(seq, Q, k)
    @info "creating the KnrIndex"
    indextime = @elapsed index = KnrIndex(dist, X; kbuild, ksearch, parallel_block, centersrecall, initial, maxiters, ordering)
    @test length(index) == length(X)
    @info "searching in the index"
    @show dim, n, m , numcenters, k, centersrecall, ordering
    tsearchtime = @elapsed Ires, Dres = searchbatch(index, Q, k)
    recall = macrorecall(Igold, Ires)
    @info "before optimization: $(index)" (recall=recall, qps=1/tsearchtime, gold_qps=1/gsearchtime)
    @info "searchtime: gold: $(gsearchtime * m), index: $(tsearchtime * m), index-construction: $indextime"

    @info "**** optimizing ParetoRadius() ****"
    opttime = @elapsed optimize!(index, ParetoRadius(); verbose=false)
    tsearchtime = @elapsed Ires, Dres = searchbatch(index, Q, k)
    recall = macrorecall(Igold, Ires)
    @info "AFTER optimization: $(index)" (recall=recall, qps=1/tsearchtime, gold_qps=1/gsearchtime)
    @info "searchtime: gold: $(gsearchtime * m), index: $(tsearchtime * m), optimization-time: $opttime"

    @info "**** optimizing ParetoRecall() ****"
    opttime = @elapsed optimize!(index, ParetoRecall(); verbose=false)
    tsearchtime = @elapsed Ires, Dres = searchbatch(index, Q, k)
    recall = macrorecall(Igold, Ires)
    @info "AFTER optimization: $(index)" (recall=recall, qps=1/tsearchtime, gold_qps=1/gsearchtime)
    @info "searchtime: gold: $(gsearchtime * m), index: $(tsearchtime * m), optimization-time: $opttime"
    @test recall >= min(0.5, minrecall)

    @info "**** optimizing MinRecall(0.95) ****"
    opttime = @elapsed optimize!(index, MinRecall(0.95); verbose=false)
    tsearchtime = @elapsed Ires, Dres = searchbatch(index, Q, k)
    recall = macrorecall(Igold, Ires)
    @info "AFTER optimization: $(index)" (recall=recall, qps=1/tsearchtime, gold_qps=1/gsearchtime)
    @info "searchtime: gold: $(gsearchtime * m), index: $(tsearchtime * m), optimization-time: $opttime"
    @test recall >= min(0.7, minrecall)

    @info "********************* generic searches *******************"
    res = KnnResult(10)
    @time search(index, Q[1], res)
    res = reuse!(res)
    @time search(index, Q[2], res)
    F = sort!(length.(index.invfile.lists), rev=true)
    println(histogram(F, nbins=20))
    println(lineplot(F))

    #=
    @info "********************* SearchGraph ***************"
    neighborhood = SimilaritySearch.Neighborhood(reduce=IdentityNeighborhood(), logbase=2)
    G = SearchGraph(; dist, db=X, neighborhood)
    G.neighborhood.reduce = IdentityNeighborhood()
    index!(G, parallel_block=1000)
    opttime = @elapsed optimize!(G, MinRecall(0.95); verbose=true)
    @time Ires, Dres, tsearchtime = timedsearchbatch(G, Q, k)
    recall = macrorecall(Igold, Ires)
    @info "AFTER optimization: $(G)" (recall=recall, qps=1/tsearchtime, gold_qps=1/gsearchtime)
    @info "searchtime: gold: $(gsearchtime * m), G: $(tsearchtime * m), optimization-time: $opttime"
    # @test recall >= min(0.7, minrecall)
    =#
end

@testset "KnrIndex" begin
    @info "********************* JIT warming *********************"  # for fast benchmarking
    centersrecall = 0.95
    runtest(; dim=2, n=10^3, m=10, k=10, centersrecall, kbuild=1, ksearch=1, ordering=InternalDistanceOrdering(), minrecall=0)
    runtest(; dim=2, n=10^3, m=10, k=10, centersrecall, kbuild=1, ksearch=1, ordering=DistanceOnTopKOrdering(10), minrecall=0)
    runtest(; dim=2, n=10^3, m=10, k=10, centersrecall, kbuild=1, ksearch=1, minrecall=0)
    # DistanceOnTopKOrdering and InternalDistanceOrdering need high `kbuild` and `ksearch` to generate
    # some discrimination power by the internal distance
    # useful for costly distance functions or whenever the dataset is pretty large
    runtest(; dim=4, n=10^5, m=1000, numcenters=10, k=100, centersrecall, kbuild=5, ksearch=5, ordering=InternalDistanceOrdering(), minrecall=0)
    @info "********************* Real search (top-k) *********************"
    # useful for costly distance functions
    runtest(; dim=4, n=10^5, m=1000, k=100, centersrecall, kbuild=5, ksearch=5, ordering=DistanceOnTopKOrdering(1000), minrecall=0.6)
    @info "********************* Real search *********************"
    # most usages
    runtest(; dim=4, n=10^6, m=1000, k=100, centersrecall, kbuild=1, ksearch=1, ordering=DistanceOrdering())
end