text:
"""
# `PREM`

1-D model containing density, radial anisotropy and attenuation.

## Reference

Dziewoński, A., Anderson, D., 1981. Preliminary reference Earth model.
Phys. Earth Planet. Inter. 25, 297–356. https://doi.org/10.1016/0031-9201(81)90046-7
"""
const PREM = PREMPolyModel(
    r =
    [1221.5, 3480.0, 3630.0, 5600.0, 5701.0, 5771.0, 5971.0, 6151.0,
     6291.0, 6346.0, 6356.0, 6368.0, 6371.0],
    vp =
    [11.2622 11.0487 15.3891 24.952 29.2766 19.0957 39.7027 20.3926 4.1875 4.1875 6.8 5.8 1.45;
    0.0 -4.0362 -5.3181 -40.4673 -23.6027 -9.8672 -32.6166 -12.2569 3.9382 3.9382 0.0 0.0 0.0;
    -6.364 4.8023 5.5242 51.4832 5.5242 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
    0.0 -13.5732 -2.5514 -26.6419 -2.5514 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    vs =
    [3.6678 0.0 6.9254 11.1671 22.3459 9.9839 22.3512 8.9496 2.1518 2.1519 3.9 3.2 0.0;
     0.0 0.0 1.4672 -13.7818 -17.2473 -4.9324 -18.5956 -4.4597 2.3481 2.3481 0.0 0.0 0.0;
     -4.4475 0.0 -2.0834 17.4575 -2.0834 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.9783 -9.2777 0.9783 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    density =
    [13.0885 12.5815 7.9565 7.9565 7.9565 5.3197 11.2494 7.1089 2.691 2.691 2.9 2.6 1.02;
     0.0 -1.2638 -6.4761 -6.4761 -6.4761 -1.4836 -8.0298 -3.8045 0.6924 0.6924 0.0 0.0 0.0;
     -8.8381 -3.6426 5.5283 5.5283 5.5283 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.0 -5.5281 -3.0807 -3.0807 -3.0807 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    vph =
    [11.2622 11.0487 15.3891 24.9520 29.2766 19.0957 39.7027 20.3926 3.5908 3.5908 6.8000 5.8000 1.4500;
    0.0000 -4.0362 -5.3181 -40.4673 -23.6027 -9.8672 -32.6166 -12.2569 4.6172 4.6172 0.0000 0.0000 0.0000;
    -6.364 4.8023 5.5242 51.4832 5.5242 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
    0.0 -13.5732 -2.5514 -26.6419 -2.5514 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    vpv =
    [11.2622 11.0487 15.3891 24.9520 29.2766 19.0957 39.7027 20.3926 0.8317 0.8317 6.8000 5.8000 1.4500;
     0.0000 -4.0362 -5.3181 -40.4673 -23.6027 -9.8672 -32.6166 -12.2569 7.2180 7.2180 0.0000 0.0000 0.0000;
     -6.364 4.8023 5.5242 51.4832 5.5242 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.0 -13.5732 -2.5514 -26.6419 -2.5514 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    vsh =
    [3.6678 0.0000 6.9254 11.1671 22.3459 9.9839 22.3512 8.9496 -1.0839 -1.0839 3.9000 3.2000 0.0000;
     0.0000 0.0000 1.4672 -13.7818 -17.2473 -4.9324 -18.5956 -4.4597 5.7176 5.7176 0.0000 0.0000 0.0000;
     -4.4475 0.0 -2.0834 17.4575 -2.0834 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.9783 -9.2777 0.9783 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    vsv =
    [3.6678 0.0000 6.9254 11.1671 22.3459 9.9839 22.3512 8.9496 5.8582 5.8582 3.9000 3.2000 0.0000;
     0.0000 0.0000 1.4672 -13.7818 -17.2473 -4.9324 -18.5956 -4.4597 -1.4678 -1.4678 0.0000 0.0000 0.0000;
     -4.4475 0.0 -2.0834 17.4575 -2.0834 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.9783 -9.2777 0.9783 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0],
    eta =
    [1. 1. 1. 1. 1. 1. 1. 1.  3.3687  3.3687 1. 1. 1.;
     0. 0. 0. 0. 0. 0. 0. 0. -2.4778 -2.4778 0. 0. 0.],
    Qμ =
    [84.6 -1. 312. 312. 312. 143. 143. 143. 80. 600. 600. 600. -1.],
    Qκ =
    [1327.7 57823. 57823. 57823. 57823. 57823. 57823. 57823. 57823. 57823. 57823. 57823. 57823.],
    fref =
    1.0
    )

"""
# `PREM_NOOCEAN`

Version of [`PREM`](@ref) where the liquid ocean layer is replaced by the
crust beneath.  This model may be useful when the surface must be a solid.
"""
const PREM_NOOCEAN = let inds = 1:(PREM.n - 1)
    PREMPolyModel(
        r = [PREM.r[1:(PREM.n - 2)]; PREM.a],
        vp = PREM.vp[:,inds],
        vs = PREM.vs[:,inds],
        density = PREM.density[:,inds],
        vph = PREM.vph[:,inds],
        vpv = PREM.vpv[:,inds],
        vsh = PREM.vsh[:,inds],
        vsv = PREM.vsv[:,inds],
        eta = PREM.eta[:,inds],
        Qμ = PREM.Qμ[:,inds],
        Qκ = PREM.Qκ[:,inds],
        fref = PREM.fref
    )
end


text:
"""
Adapted from sandreza/Learning/sandbox/oceananigans_converter.jl
https://github.com/sandreza/Learning/blob/master/sandbox/oceananigans_converter.jl
"""

using JLD2, NetCDF

struct LESbraryData{𝒮, 𝒯, 𝒰, 𝒱}
    # initial conditions, 4
    T⁰::𝒮
    # S⁰::𝒮
    U⁰::𝒮
    V⁰::𝒮

    # fields at each moment in time, 4
    T::𝒯
    # S::𝒯
    U::𝒯
    V::𝒯

    # some second order statistics at each moment in time, 5
    wT::𝒯
    # wS::𝒯
    uu::𝒯
    vv::𝒯
    ww::𝒯
    # tke::𝒯

    # simulation constants, 8
    # ρ::𝒰
    α::𝒰
    β::𝒰
    # cᵖ::𝒰
    f⁰::𝒰
    g::𝒰
    L::𝒰
    # κₑ::𝒰

    # time and grid, 2
    t::𝒮
    z::𝒮

    # #top boundary condition data, see string for type of boundary condition, 4
    # top_T::𝒰
    # top_S::𝒰
    # top_U::𝒰
    # top_V::𝒰
    #
    # #bottom boundary condition data, see string for type of boundary condtion,4
    # bottom_T::𝒰
    # bottom_S::𝒰
    # bottom_U::𝒰
    # bottom_V::𝒰

    #info about the simulation, 1
    info::𝒱
end

"""
LESbraryData(filename)
# Description
- Constructor for Oceananigans data type. Loads data from LES
# Fields for the output are
    # initial conditions
    T⁰::𝒮
    S⁰::𝒮
    U⁰::𝒮
    V⁰::𝒮
    # fields at each moment in time
    T::𝒯
    S::𝒯
    U::𝒯
    V::𝒯
    # some second order statistics at each moment in time
    wT::𝒯
    wS::𝒯
    uu::𝒯
    vv::𝒯
    ww::𝒯
    tke::𝒯
    # simulation constants
    ρ::𝒰
    α::𝒰
    β::𝒰
    cᵖ::𝒰
    f⁰::𝒰
    g::𝒰
    # time and grid
    t::𝒮
    z::𝒮
    # #top boundary condition data, see string for type of boundary condition
    # top_T::𝒰
    # top_S::𝒰
    # top_U::𝒰
    # top_V::𝒰
    # #bottom boundary condition data, see string for type of boundary condtion
    # bottom_T::𝒰
    # bottom_S::𝒰
    # bottom_U::𝒰
    # bottom_V::𝒰
    #info about the simulation
    info::𝒱
"""
function ReadJLD2_LESbraryData(filename)
    les_data = jldopen(filename, "r")
    les_keys = keys(les_data)
    timeseries_keys = keys(les_data["timeseries"]["t"])

    # hold the entries for easy constructor creation
    container = []

    # println(collect(les_data["grid"]["zC"]))

    # size of arrays
    Nz = length(collect(les_data["grid"]["zC"]))-6

    Nt = length(timeseries_keys)

    ## construct arrays
    #Initial Conditions
    T⁰ = zeros(Nz)
    # S⁰ = zeros(Nz)
    U⁰ = zeros(Nz)
    V⁰ = zeros(Nz)

    #Timeseries
    T = zeros(Nz, Nt)
    # S = zeros(Nz, Nt)
    U = zeros(Nz, Nt)
    V = zeros(Nz, Nt)
    t = zeros(Nt)

    #Second Order Statistics
    wT  = zeros(Nz+1, Nt)
    # wS  = zeros(Nz, Nt)
    uu  = zeros(Nz, Nt)
    vv  = zeros(Nz, Nt)
    ww  = zeros(Nz+1, Nt)

    # grab arrays
    for j in 1:Nt
        # Fields
        key = timeseries_keys[j]
        @. T[:,j] = les_data["timeseries"]["T"][key][1,1,:]
        # @. S[:,j] = les_data["timeseries"]["S"][key][1,1,:]
        @. U[:,j] = les_data["timeseries"]["u"][key][1,1,:]
        @. V[:,j] = les_data["timeseries"]["v"][key][1,1,:]
        # Second Order Statistics
        @. wT[:,j] = les_data["timeseries"]["wT"][key][1,1,:]
        # @. wS[:,j] = les_data["timeseries"]["wS"][key][1,1,:]
        @. uu[:,j] = les_data["timeseries"]["uu"][key][1,1,:]
        @. vv[:,j] = les_data["timeseries"]["vv"][key][1,1,:]
        @. ww[:,j] = les_data["timeseries"]["ww"][key][1,1,:]

        t[j] = les_data["timeseries"]["t"][key]
    end

    # tke = 0.5*(uu .+ vv .+ ww)

    # Set initial Conditions
    @. T⁰ = T[:,1]
    # @. S⁰ = S[:,1]
    @. U⁰ = U[:,1]
    @. V⁰ = V[:,1]

    # Push initial conditions current stuff into container
    # push!(container, T⁰, S⁰, V⁰, U⁰)
    push!(container, T⁰, V⁰, U⁰)
    # Push fields into container
    # push!(container, T, S, U, V)
    push!(container, T, U, V)
    # Push second order statistics into container
    # push!(container, wT, wS, uu, vv, ww, tke)
    push!(container, wT, uu, vv, ww)

    # Now grab parameter
    # ρ = les_data["closure"]["density"]
    α = les_data["buoyancy"]["equation_of_state"]["α"]
    β = les_data["buoyancy"]["equation_of_state"]["β"]
    # cᵖ = les_data["parameters"]["specific_heat_capacity"]
    f⁰ = les_data["coriolis"]["f"]
    g = les_data["buoyancy"]["gravitational_acceleration"]
    L = les_data["grid"]["Lz"]
    # κₑ = les_data["parameters"]["diffusivity_T"]

    # b = T .* (g * α)

    # Push parameters to container
    # push!(container, ρ, α, β, cᵖ, f⁰, g, L, κₑ)
    push!(container, α, β, f⁰, g, L)

    # grab domain data
    z = collect(les_data["grid"]["zC"])[4:end-3] # padding of 3 on each side

    # push
    push!(container, t, z)

    print(les_keys)
    # now grab boundary condition data
    # top_T = les_data["boundary_conditions"]["top"]["FT"]
    # top_S = 0.0
    # top_U = les_data["boundary_conditions"]["top"]["Fu"]
    # top_V = 0.0
    #bottom boundary condition data, see string for type of boundary condtion
    # bottom_T = les_data["boundary_conditions"]["bottom"]["dTdz"]
    # bottom_S = 0.0
    # bottom_U = 0.0
    # bottom_V = 0.0

    # push to container
    # push!(container, top_T, top_S, top_U, top_V, bottom_T, bottom_S, bottom_U, bottom_V)

    # Now construct types
    𝒮 = typeof(T⁰)
    𝒯 = typeof(T)
    𝒰 = typeof(α)
    𝒱 = typeof("string")

    # now create data string
    info_string = "The top boundary conditions are flux boundary conditions \n"
    info_string *= "The  bottom boundary condition for temperature is a gradient boundary condition \n"
    info_string *= "The grid data is assumed to be evenly spaced and a power of two \n"

    # push to container
    push!(container, info_string)
    #return container
    close(les_data)
    return LESbraryData{𝒮, 𝒯, 𝒰, 𝒱}(container...)
end

# avg = "src/les/data/2daySuite/three_layer_constant_fluxes_Qu0.0e+00_Qb1.0e-07_f1.0e-04_Nh128_Nz128/three_layer_constant_fluxes_Qu0.0e+00_Qb1.0e-07_f1.0e-04_Nh128_Nz128_averaged_statistics.jld2"
# ReadJLD2_LESbraryData(avg)

# avg = "src/les/data/general_strat_sandreza/general_strat_4_profiles/general_strat_4_profiles.jld2"

# les_data = jldopen(avg, "r")
# les_keys = keys(les_data)
# timeseries_keys = keys(les_data["timeseries"]["t"])

# stats = "src/les/data/2daySuite/three_layer_constant_fluxes_Qu1.0e-03_Qb0.0e+00_Nh128_Nz128/three_layer_constant_fluxes_Qu1.0e-03_Qb0.0e+00_Nh128_Nz128_statistics.jld2"
# ReadJLD2_OceananigansData(stats)

ext: jl

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
# Naive Bayes Module

module NB
export NaiveBayes, train, predict

"""
Data structure for Naive Bayes\\
`pY` is probability for each target\\
`mean` is the corresponding mean of Xi for each taregt\\
`var` is the corresponding variance of Xi for each target
"""
mutable struct NaiveBayes
    n_features::Integer
    pY::Dict{Number,AbstractFloat}
    mean::Dict{Number,Array{AbstractFloat}}
    var::Dict{Number,Array{AbstractFloat}}
end

"""
Training function for Naive Bayes\\
`X_data` has shape (M,N)\\
`Y_data` has shape (M,)\\
Returns the computed data structure `NaiveBayes`
"""
function train(X_data::Array{T} where T<:Number, Y_data::Array{K} where K<:Number)::NaiveBayes
    @assert ndims(X_data) == ndims(Y_data) + 1
    @assert size(X_data)[1] == size(Y_data)[1]
    n_features = size(X_data)[2]
    # compute probability of Y
    pY = Dict{Number,AbstractFloat}()
    for Y_val in Y_data
        if haskey(pY, Y_val)
            pY[Y_val] += 1.0
        else
            pY[Y_val] = 1.0
        end
    end
    for Y_val in keys(pY)
        pY[Y_val] /= size(Y_data)[1]
    end
    # compute Xi mean and var for each Y
    mean = Dict{Number,Array{AbstractFloat}}()
    var = Dict{Number,Array{AbstractFloat}}()
    for Y_val in keys(pY)
        X_part = X_data[Y_data .== Y_val, :]
        mean[Y_val] = vec(sum(X_part, dims=1) ./ size(X_part)[1])
        var[Y_val] = vec(sum((X_part .- reshape(mean[Y_val], 1, :)).^2, dims=1) ./ size(X_part)[1])
    end
    return NaiveBayes(n_features, pY, mean, var)
end

"""
Vectorized gaussian function\\
Returns the product of P(xi|y)
"""
function gaussian(X_vec::Array, mean::Array, var::Array)::AbstractFloat
    @assert ndims(X_vec) == ndims(mean) == ndims(var) == 1
    @assert length(X_vec) == length(mean) == length(var)
    left = 1.0 ./ sqrt.((2.0 * pi) .* var)
    right = exp.(-(X_vec .- mean).^2 ./ (2.0 .* var))
    p = left .* right
    return prod(p)
end

"""
Predict function for Naive Bayes\\
`X_data` has ndims of 1 or 2
"""
function predict(X_data::Array{T} where T<:Number, data::NaiveBayes)::Array
    if ndims(X_data) == 1
        X_data = reshape(X_data, 1, :)
    end
    @assert ndims(X_data) == 2
    @assert size(X_data)[2] == data.n_features
    prediction = Array{Number}(undef, size(X_data)[1])
    for i in 1:size(X_data)[1]
        X_vec = X_data[i, :]
        Y_dict = Dict{Number,AbstractFloat}()
        for Y_val in keys(data.pY)
            Y_dict[Y_val] = data.pY[Y_val] * gaussian(X_vec, data.mean[Y_val], data.var[Y_val])
        end
        prediction[i] = sort(collect(Y_dict), by=m->m[2], rev=true)[1][1]
    end
    return prediction
end

end

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
module PlasmaBenchmarks

using Plasma
using BSON: @save

include("two-stream.jl")
include("weibel.jl")
# include("1d1v.jl")
# include("3d3v.jl")

end


text:
randuniform(dims::Base.Dims) = randuniform(Float64, dims)
randuniform(::Type{T}, dims::Base.Dims) where {T<:Number} = rand(T, dims)

randnormal(dims::Base.Dims) = randnormal(Float64, dims)
randnormal(::Type{T}, dims::Base.Dims) where {T<:Number} = randn(T, dims)

randisometry(dims::Base.Dims{2}) = randisometry(Float64, dims)
randisometry(::Type{T}, dims::Base.Dims{2}) where {T<:Number} =
    dims[1] >= dims[2] ? _leftorth!(randnormal(T, dims), QRpos(), 0)[1] :
        throw(DimensionMismatch("cannot create isometric matrix with dimensions $dims; isometry needs to be tall or square"))

const randhaar = randisometry


text:
using Test
include("module-2-exercise-2.jl")

@testset "sin tests" begin
    x = pi
    @test isapprox(sin(x), calculate_series_sin(x); atol = 0.05)
    x = 2*pi
    @test isapprox(sin(x), calculate_series_sin(x); atol = 0.05)
    x = 0
    @test isapprox(sin(x), calculate_series_sin(x); atol = 0.05)
    x = pi/7
    @test isapprox(sin(x), calculate_series_sin(x); atol = 0.05)
end


text:
isdefined(Base, :__precompile__) && __precompile__()
module PowerSimulationsDynamics

####################################### Structs Exports ####################################

# Base Exports
export ResidualModel
export MassMatrixModel

export Simulation
export Simulation!
export SimulationResults

export execute!
export read_results

# Frequency reference options
export ReferenceBus
export ConstantFrequency

# export perturbations
export NetworkSwitch
export ControlReferenceChange
export BranchTrip
export BranchImpedanceChange
export SourceBusVoltageChange
export GeneratorTrip
export LoadTrip
export LoadChange
# export BusTrip

# Export for routines
export small_signal_analysis
export get_state_series
export get_voltage_magnitude_series
export get_voltage_angle_series
export show_states_initial_value
export read_initial_conditions
export get_real_current_series
export get_imaginary_current_series
export get_activepower_series
export get_reactivepower_series
export get_setpoints

####################################### Package Imports ####################################
import Logging
import InfrastructureSystems
import SciMLBase
import DataStructures: OrderedDict
import Random
import DiffEqBase
import ForwardDiff
import SparseArrays
import LinearAlgebra
import Base.to_index
import NLsolve
import PrettyTables
import Base.ImmutableDict
import PowerSystems
import TimerOutputs
const PSY = PowerSystems
const IS = InfrastructureSystems
const PSID = PowerSimulationsDynamics

using DocStringExtensions

#@template (FUNCTIONS, METHODS) = """
#                                 $(TYPEDSIGNATURES)
#                                 $(DOCSTRING)
#                                 """

#Structs for General Devices and System
include("base/definitions.jl")
include("base/ports.jl")
include("base/bus_categories.jl")
include("base/load_categories.jl")
include("base/device_wrapper.jl")
include("base/branch_wrapper.jl")
include("base/frequency_reference.jl")
include("base/file_system.jl")
include("base/simulation_model.jl")
include("base/simulation_inputs.jl")
include("base/perturbations.jl")
include("base/caches.jl")
include("base/system_model.jl")
include("base/jacobian.jl")
include("base/mass_matrix.jl")
include("base/simulation_results.jl")
include("base/simulation.jl")
include("base/global_variables.jl")
include("base/supplemental_accesors.jl")
include("base/nlsolve_wrapper.jl")
include("base/simulation_initialization.jl")
include("base/small_signal.jl")

#Common Models
include("models/branch.jl")
include("models/device.jl")
include("models/network_model.jl")
include("models/dynline_model.jl")
include("models/ref_transformations.jl")

#Generator Component Models
include("models/generator_models/machine_models.jl")
include("models/generator_models/pss_models.jl")
include("models/generator_models/avr_models.jl")
include("models/generator_models/tg_models.jl")
include("models/generator_models/shaft_models.jl")

#Inverter Component Models
include("models/inverter_models/DCside_models.jl")
include("models/inverter_models/filter_models.jl")
include("models/inverter_models/frequency_estimator_models.jl")
include("models/inverter_models/outer_control_models.jl")
include("models/inverter_models/inner_control_models.jl")
include("models/inverter_models/converter_models.jl")

#Injection Models
include("models/load_models.jl")
include("models/source_models.jl")

#Saturation Models
include("models/saturation_models.jl")

#Initialization Parameters
include("initialization/init_device.jl")

#Initialization Generator
include("initialization/generator_components/init_machine.jl")
include("initialization/generator_components/init_shaft.jl")
include("initialization/generator_components/init_avr.jl")
include("initialization/generator_components/init_tg.jl")
include("initialization/generator_components/init_pss.jl")

#Initialization Inverter
include("initialization/inverter_components/init_filter.jl")
include("initialization/inverter_components/init_DCside.jl")
include("initialization/inverter_components/init_converter.jl")
include("initialization/inverter_components/init_frequency_estimator.jl")
include("initialization/inverter_components/init_inner.jl")
include("initialization/inverter_components/init_outer.jl")

#System Model
include("models/system.jl")

#PostProcessing
include("post_processing/post_proc_common.jl")
include("post_processing/post_proc_generator.jl")
include("post_processing/post_proc_inverter.jl")
include("post_processing/post_proc_results.jl")

#Utils
include("utils/psy_utils.jl")
include("utils/immutable_dicts.jl")
include("utils/print.jl")
include("utils/kwargs_check.jl")
include("utils/logging.jl")

end # module


text:
#
# Proof of concept using artificial data.
# Illustrate usefulness of the PCCF on artificially generated signal.
#

include("./lib/metadata.jl")
include("./2bayes_detector_replica.jl")
include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")
include("./lib/pccf.jl")

using SignalGenerator
using BayesChpOrig
using Metadata
using Roc
using Pccf

function run_sig_gen()
    ret = genSig(5, 100.0, 1.0, 2.0, 0.5, 1)
    sig = ret["sig"]
    chps= ret["changes"]
    writedlm(outdir*"testsig.txt", sig)
    writedlm(outdir*"testsigchps.txt", chps)
    print("run1")
end

function run_proof3(outputfolder)
    println("... run_proof3 > output folder = ", outputfolder)
    # genSig(n, mu, sigma, changeStep, sigStd)
    mu    = 100.0
    sigma = 10.0
    srand(1234)
    ret   = genSig(10, mu, sigma, 1.3, 0.75)
    sig   = ret["sig"]
    n     = length(sig)
    chps  = ret["changes"]
    probs = scale01(pccf(n, mu, sigma))
    # pccf_matrix_max_vals, pccf_matrix = pccf_bayes_detector(sig, 200 * probs, 200, 200 + 30, 0.15) # one outlier less
    # start INTERESTING EDGE CASE
    #pccf_matrix_max_vals, pccf_matrix = pccf_bayes_detector(sig, 200 * probs, 200, 200 + 150, 0.3) #150-max,min-110.
    # -- end INTERESTING EDGE CASE
    # start One outlier is removed case
    pccf_matrix_max_vals, pccf_matrix = pccf_bayes_detector(sig, probs, 200+100, 200, 0.3)
    # -- end One outlier is removed case
    matrix_max_vals, matrix           = bayes_detector(sig, lambda = 200)
    detections                        = extract_changes(matrix_max_vals)
    pccf_detections                   = extract_changes(pccf_matrix_max_vals)
    #
    writedlm(outputfolder * "pccf.txt", probs)
    writedlm(outputfolder * "changes.txt", chps)
    writedlm(outputfolder * "signal.txt", sig)
    writedlm(outputfolder * "detections.txt", detections)
    writedlm(outputfolder * "pccf_detections.txt", pccf_detections)
    if 1 == 0
        max_delay = 30
        tpr = tpRate(detections, chps, max_delay)
        fpr = faRateScaled(detections, chps, max_delay, length(sig), tn_price_scale=0.1)
        println("... test_bd > TPrate = ", tpr, " / FArate = ", fpr)
        ##writedlm(outdir*"concept_proof/matrix.txt", matrix, ",")
        ##writedlm(outdir*"concept_proof/matrix_max_vals.txt", matrix_max_vals, ",")
        # println("... test_bd > changes:", chps)
        # println("... test_bd > detections:", detections)
    end
end

run_proof3(outdir * "concept_proof3/")

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
# test show
@testset "show" begin
(path, io) = mktemp()
# show
redirect_stdout(io) do
@test !isa(try show(RawData()) catch ex ex end, Exception)
@test !isa(try show(FFTData()) catch ex ex end, Exception)
@test !isa(try show(CorrData()) catch ex ex end, Exception)
# summary
@test !isa(try summary(RawData()) catch ex ex end, Exception)
@test !isa(try summary(FFTData()) catch ex ex end, Exception)
@test !isa(try summary(CorrData()) catch ex ex end, Exception)
end
end


text:
module SSA

    # -------
    using RandomExtensions
    using Distributed

    include("Types.jl")
    using .Types

    import Dates

    using Random

    import Base, Base.show, Base.isless

    export  runssa,
            isless,
            freeze,
            defrost,
            generation_num,
            population

    # -------

    #isless(lhs::Entity, rhs::Entity) = lhs.fitness < rhs.fitness

    fitness!(ent, fitness_score) = ent.fitness = fitness_score

    distance(ind1 , ind2 ) = abs(ind1.fitness - ind2.fitness)

    global _g_model

    # -------

    function freeze(model::GAmodel, entity::EntityData)
        push!(model.freezer, entity)
        println("Freezing: ", entity)
    end

    function freeze(model::GAmodel, entity)
        entitydata = EntityData(entity, model.params.currentGeneration)
        freeze(model, entitydata)
    end

    freeze(entity) = freeze(_g_model, entity)

    function defrost(model::GAmodel, generation::Int)
        filter(model.freezer) do entitydata
            entitydata.generation == generation
        end
    end

    defrost(generation::Int) = defrost(_g_model, generation)

    generation_num(model::GAmodel = _g_model) = model.params.currentGeneration

    population(model::GAmodel = _g_model) = model.population

    function show_simulation(model :: GAmodel, ent)
        printed = model.specific_fitness.simulate_entity(ent, model.instructionsSet)
        @show ent.fitness
        #@show model.specific_fitness.fitness(ent)
        #println(printed)
        return "fitness : $(ent.fitness) \n $printed"
    end

    function findBestNotIn(ens, i)
        if i in ens
            return findBestNotIn(ens, i+1)
        else
            return i
        end
    end

    function subSetsGeneration(n, refSet)
        subSets = []
        currentTypeSet = []
        oldTypeSet = []
        lref = length(refSet)
        #Type2
        for i in 1:(n-1)
            map( a -> (push!(subSets, [a,i]); push!(currentTypeSet, [i, a]) ), i+1:n )
        end
        # Type3
        oldTypeSet = currentTypeSet
        oldTypeSet = map( a -> (i = findBestNotIn(a, 1); push!(subSets, [a;i]); [a;i]), oldTypeSet )
        #Type 4
        oldTypeSet = map( a -> (i = findBestNotIn(a, 1); push!(subSets, [a;i]); [a;i]), oldTypeSet )
        # Type 5
        map(a -> push!(subSets, 1:a), 5:n)
        map( s -> map(e -> refSet[(lref - e +1)] , s) , subSets)
    end

    function searchInCurrentRefSet_Basic(model)
        newSolutions = true
        j = 0
        while newSolutions
            #=index = model.ga.addNew(model.refSet, model.population)
            newE = model.population[index]
            deleteat!(model.population, index)
            push!(model.population, popfirst!(model.refSet))
            push!(model.refSet, newE)=#
            # 3. Generate NewSubsets with the subset generation method. Make NewSolutions = FALSE.
            lrefSet = length(model.refSet)
            model.subSets = subSetsGeneration(lrefSet, model.refSet)
            newSolutions = false
            while !isempty(model.subSets)
                #evaluate_refSet(model)
                #popfirst!(model.population, )
                # 4. Select the next subset s in NewSubsets
                s = pop!(model.subSets)
                # 5. Apply the solution combination method to s to obtain one or more new trial solutions x.
                # Apply the improvement method to the trial solutions.
                pool = []

                #=ls = length(s)
                for i in 1:ls-1
                    map( a -> ((child1, child2) = model.ga.crossover(Any[model.refSet[a],model.refSet[i]]);
                            push!(pool, child1); push!(pool, child2) ) , i+1:ls )
                end=#
                #group = map( j -> model.refSet[j], s)
                (child1, child2) = model.ga.crossoverGroup(s)
                push!(pool, child1)
                push!(pool, child2)

                pmap(
                    ent -> (model.ga.entityToBfInstructions!(ent); model.specific_fitness.fitness(ent, model.instructionsSet))
                    , pool)
                #sort!(pool, rev=true)
                # 6. Apply the reference set update method.
                newEntAdded = []
                refSetHasChanged = false
                #@show length(model.refSet)
                pmap( ent -> model.ga.clearCode!(ent), pool)
                for e in pool
                    if !isempty(model.refSet) && (model.refSet[1]) < e
                        model.refSet[1] = e
                        evaluate_refSet(model)
                        sort!(model.refSet; rev = false)
                        #push!(newEntAdded, e)
                        #deleteat!(model.refSet, 1)
                        refSetHasChanged = true
                    end
                end
                #append!(model.refSet, newEntAdded)
                evaluate_refSet(model)
                sort!(model.refSet; rev = false)

                lastIdx = length(model.refSet)
                _fitness = model.refSet[lastIdx].fitness
                # displays
                j += 1
                displayStep!(model, j)
                if (model.params.targetFitness > 0 && _fitness >= model.params.targetFitness)
                    model.params.targetFitnessCount = model.params.targetFitnessCount +1
                    if (model.params.targetFitnessCount > 500) # default : set to 1000
                        # Stop the algo
                        newSolutions = false
                        break;
                    end
                else
                    # in case we expand and lose the best one
                    model.params.targetFitnessCount = 0;
                end

                if refSetHasChanged
                    # 7. Make NewSolutions = TRUE.
                    newSolutions = true
                end
            end  #while
            #evaluate_refSet(model)
        end #while
    end

    function displayStep!(model, j; msg = "Iteration")
        lastIdx = length(model.refSet)
        _fitness = model.refSet[lastIdx].fitness
        if j % 5 == 0
            #@show refSetHasChanged
            #_log = ""
            _log = "    $(Dates.now()) , "
            _log *= "$msg : $j, " #": $(model.params.currentGeneration) , "
            _log *= "BEST: $_fitness , bonus = $(model.refSet[lastIdx].bonus) \n"

            if j%40 == 0
                _log *= show_simulation(model, model.refSet[lastIdx]) * "\n"
            end
            print(_log)
            if model.params.historyPath != nothing
                write(model.params.historyPath, _log)
            end
        end

    end

    function searchInCurrentRefSet_RelinkingPath(model)
        NumImpl = 10 # To adapt after
        #1. already done in model.refSet
        # 2.
        lrefSet = length(model.refSet)
        #model.subSets = subSetsGeneration(lrefSet)
        j = 0
        newSolutions = true
        while newSolutions
            #3.
            model.subSets = []
            subSetsInt = []
            for i in 1:(lrefSet-1)
                map( a -> push!(subSetsInt, [a,i]), i+1:lrefSet )
            end
            model.subSets = map( s -> map( (e -> model.refSet[e]) , s) , subSetsInt)
            newSolutions = false
            # displays
            j += 1
            displayStep!(model, j)
            while !isempty(model.subSets)
                #evaluate_refSet(model)
                pool = []
                lastIdx = length(model.refSet)
                _fitness = model.refSet[lastIdx].fitness
                if (model.params.targetFitness > 0 && _fitness >= model.params.targetFitness)
                    model.params.targetFitnessCount = model.params.targetFitnessCount +1
                    if (model.params.targetFitnessCount > 50) # default : set to 1000
                        # Stop the algo
                        newSolutions = false
                        break;
                    end
                else
                    # in case we expand and lose the best one
                    model.params.targetFitnessCount = 0;
                end
                #4.
                xxp = pop!(model.subSets)
                #5.
                pool2 = []
                (child1, child2) = model.ga.crossover(xxp)
                push!(pool, child1)
                push!(pool, child2)
                relinking!(xxp[1], xxp[2], pool2, model)
                #relinking2!(xxp[1], xxp[2], pool2, model)
                append!(pool, pool2)
                r= length(pool2)
                for i in 1:NumImpl:r
                    # 6.
                    newS = model.ga.create_entity(pool2[i].dna)
                    model.ga.mutate(newS, 0.05)
                    push!(pool, newS)
                end # for
                #7.
                pool2 = []
                relinking!(xxp[2], xxp[1], pool2, model)
                append!(pool, pool2)
                s= length(pool2)
                for i in 1:NumImpl:s
                    #8..
                    newS = model.ga.create_entity(pool2[i].dna)
                    model.ga.mutate(newS, 0.05)
                    push!(pool, newS)
                end # for
                #= An experiment part
                pool2 = []
                relinking2!(xxp[1], xxp[2], pool2, model)
                append!(pool, pool2)
                s= length(pool2)
                for i in 1:NumImpl:s
                    #8..
                    newS = model.ga.create_entity(pool2[i].dna)
                    model.ga.mutate(newS, 0.05)
                    push!(pool, newS)
                end # for
                # An experiment part =#
                #append!(pool, pool2)
                pmap( ent -> model.ga.clearCode!(ent), pool)
                pmap(
                    ent -> (model.ga.entityToBfInstructions!(ent); model.specific_fitness.fitness(ent, model.instructionsSet))
                    , pool)
                #=println("pool")
                map(e -> println(e.program) , pool)
                println()=#
                #evaluate_refSet(model)
                sort!(model.refSet; rev = false)
                for x in pool
                    if model.refSet[1] < x && !model.ga.isPresent(x,model.refSet)
                        #print("hi")
                        # 9.
                        model.refSet[1] = x
                        sort!(model.refSet; rev = false)
                        #10.
                        newSolutions = true
                    end # if
                end # for
                #11. Already done
            end #while
        end #while
    end

    function runssa(mdl::Module, fit_mdl :: Module )
        # get the parameters of the specific fitness
        model = GAmodel(fit_mdl.getParams())
        model.specific_fitness = fit_mdl

        # save the ga (mutation and crossover units function)
        model.ga = mdl
        # save once the dictionary of instructions Set
        model.instructionsSet = mdl.getInstructionsSet()
        #run the ga
        runssa(model; resume = false)
    end

    function runssa(model::GAmodel; resume = false)
        stop = false

        #set _expandAmount & _expandRate, TODO I think they are not placed very well in case we want specifi paramters for each fitness
        _expandAmount = 5
        _expandRate = 500

        if (!resume)
            #  initialize params.
            model.params.totalFitness = 0;
            model.params.targetFitness = model.specific_fitness.getTargetFitness();
            model.params.targetFitnessCount = 0;
            model.params.currentGeneration = 0;
        end

        if model.params.historyPath != nothing
            #to save params in the history file
            write(model.params.historyPath, "params : " * string(model.params) * " \n")
        end

        # beginning of the scatter implementation
        #=
        1. Start with P = Ø. Use the diversification generation method to construct a solution
        and apply the improvement method. Let x be the resulting solution.
        If x ∉ P then add x to P (i.e., P = P ∪ x ), otherwise, discard x.
        Repeat this step until |P| = PSize.
        =#
        lastIdx = model.params.populationSize #length(model.population)
        reset_model(model)

        create_initial_population(model)
        #pmap( ent -> model.ga.clearCode!(ent), model.population)
        evaluate_population(model)
        for i in 1:5
            push!(model.refSet, pop!(model.population))
        end
        #=
        2. Use the reference set update method to build RefSet = { x 1 , ..., x b } with the “best” b solutions in P.
         Order the solutions in RefSet according to their objective function value such that x 1 is
         the best solution and x b the worst. Make NewSolutions = TRUE.
        =#

        found = false
        tour = 0
        j = 0
        while !found
            tour += 1
            _log = "** Tour n. $tour at $(Dates.now()) \n"
            displayStep!(model, tour; msg ="Tour")
            print(_log)
            if model.params.historyPath != nothing
                write(model.params.historyPath, _log)
            end
            if (_expandAmount > 0 && tour % _expandRate == 0 && model.params.genomeSize < model.params.maxGenomeSize )
                model.params.genomeSize +=  _expandAmount;
                #_bestStatus.Fitness = 0; # Update display of best program, since genome has changed and we have a better/worse new best fitness.
            end
            if model.population[1].m_length != model.params.genomeSize
                #println("expannndd")
                newGenomeSize = model.params.genomeSize
                for m in model.population
                    if m.m_length != newGenomeSize
                        model.ga.expand(m, newGenomeSize)
                    end
                end
                for m in model.refSet
                    if m.m_length != newGenomeSize
                        model.ga.expand(m, newGenomeSize)
                    end
                end

            end

            for i in 1:10
                index = model.ga.addNew(model.refSet, model.population[1:div(model.params.populationSize,2)])
                push!(model.refSet, model.population[index])
                deleteat!(model.population, index)
            end

            evaluate_refSet(model)
            #@show length(model.refSet)
            #=if rand()< 0.5
                searchInCurrentRefSet_Basic(model)
            else
                searchInCurrentRefSet_RelinkingPath(model)
            end=#
            searchInCurrentRefSet_RelinkingPath(model)
            best = model.refSet[length(model.refSet)]
            found = best.fitness >= model.params.targetFitness
            if !found
                #append!(model.population, model.refSet)
                empty!(model.population)
                create_initial_population(model)
                evaluate_population(model)
                model.refSet = model.refSet[11:15]
            end

        end

        print("\n **********DONE*******")
        if model.params.historyPath != nothing
            best = model.refSet[length(model.refSet)] #[length(model.population)]
            _log = show_simulation(model, best) * "\n"
            #_log *= "\n Generation : $(model.params.currentGeneration) \n"
            _log *= "fitness $(best.fitness) \n"
            _log *= "dna : $(best.dna) \n"
            print(_log)
            write(model.params.historyPath, _log)
            close(model.params.historyPath)
        end
        model

    end

    # -------

    function reset_model(model::GAmodel)
        global _g_model = model

        model.params.currentGeneration = 1
        empty!(model.population)
    end

    function create_initial_population(model::GAmodel)
        for i = 1:model.params.populationSize
            entity = model.ga.create_entity(i, model.params.genomeSize)
            push!(model.population, entity)
        end
    end

    function evaluate_population(model::GAmodel)
        scores = pmap(
            ent -> (model.ga.entityToBfInstructions!(ent); model.specific_fitness.fitness(ent, model.instructionsSet))
            , model.population)

        model.params.totalFitness = sum(scores)

        # the isless for the entities is defined in bfga
        sort!(model.population; rev = false)
        model.scores = sort(scores; rev = false)
    end

    function evaluate_refSet(model::GAmodel)
        scores = pmap(
            ent -> (model.ga.entityToBfInstructions!(ent); model.specific_fitness.fitness(ent, model.instructionsSet))
            , model.refSet)
        sort!(model.refSet; rev = false)
    end

    #=function crossover_population(model::GAmodel, groupings)
        thisGeneration = (population(model))
        _length = 1

        model.population = Any[]

        model.params.currentGeneration += 1
        g = model.ga.create_entity(model.params.currentGeneration, model.params.genomeSize)
        g2 = model.ga.create_entity(model.params.currentGeneration, model.params.genomeSize)
        #println("Generation n' $(model.params.currentGeneration) ")
        if model.params.elitism
            l = model.params.populationSize #length(thisGeneration)
            g = thisGeneration[l]
            l_1 = (model.params.populationSize) -1
            g2 = thisGeneration[l-1]

            _length += 2;
        end
        #println("_______________________________________________________________________________")
        #@show model.scores

        for i in _length:2:model.params.populationSize
            pidx1 = rouletteSelection(model)
            pidx2 = rouletteSelection(model)
            #println("$pidx1    $pidx2") # for debugging
            child1, child2 = nothing, nothing
            parent1 = thisGeneration[pidx1]
            parent2 = thisGeneration[pidx2]

            if (rand() < model.params.crossoverRate)
                child1, child2 = model.ga.crossover(Any[parent1, parent2])
            else
                child1, child2 = parent1, parent2
                parent1.bonus = 0
                parent2.bonus = 0
            end

            push!(model.population, child1)
            push!(model.population, child2)
        end

        if model.params.elitism
            push!(model.population, g)
            #push!(model.pop_data, EntityData(g, model.params.currentGeneration))
            push!(model.population, g2)
            #push!(model.pop_data, EntityData(g2, model.params.currentGeneration))
        end

        model.params.populationSize = length(model.population)
        #pmap(model.ga.clearCode!, model.population)

        # Expand genomes.
        if model.population[1].m_length != model.params.genomeSize
            #println("expannndd")
            newGenomeSize = model.params.genomeSize
            for m in model.population
                if m.m_length != newGenomeSize
                    model.ga.expand(m, newGenomeSize)
                end
            end
        end
    end=#

    #=function mutate_population(model::GAmodel) # deprecated, the mutation is performed in the crossover
        pmap(e -> model.ga.mutate( e, 0.01), model.population)
        #=for entity in model.population
            model.ga.mutate(entity)
        end=#
    end

    # We can add strangers in the population
    function add_stranger(model :: GAmodel, num)
        #println("-_-")
        for i in 1:num
            entity = model.ga.create_entity(model.params.currentGeneration, model.population[1].m_length)
            push!(model.population, entity)
            #@async push!(model.pop_data, EntityData(entity, model.params.currentGeneration))
        end
    end=#

    function rouletteSelection(model :: GAmodel)

        n = model.params.populationSize #length(model.scores)
        randomFitness = rand() * (model.scores[n] == 0 ? 1 : model.scores[n])

        idx = -1
        first = 1
        last = model.params.populationSize
        mid = trunc(Int , (last - first)/2)

        #  ArrayList's BinarySearch is for exact values only
        #  so do this by hand.
        while (idx == -1 && first <= last)
            if (randomFitness < model.scores[mid])
                last = mid;
            elseif (randomFitness > model.scores[mid])
                first = mid;
            end
            mid = trunc(Int, (first + last)/2)
            #  lies between i and i+1
            if ((last - first) == 1)
                idx = last;
            end
        end

        return idx
    end

    function relinking_aux!(x1, x2, pool, model, ens_val, d, nb)
        if !isempty(ens_val)
            #@show ens_val
            r = pop!(ens_val)
            newS = model.ga.create_entity(x1.dna)
            sameRange = false
            if r == d -1
                newS.dna[(r*nb+1):end] = x2.dna[(r*nb+1):end]
            else
                newS.dna[(r*nb+1):((r+1)*nb)] = x2.dna[(r*nb+1):((r+1)*nb)]
            end
            push!(pool, newS)
            relinking_aux!(newS, x2, pool, model, ens_val, d, nb)
        end
    end

    function relinking3!(x1, x2, pool, model)
        #println(length(pool))
        nb = 6
        d = div(x2.m_length, nb)
        X = collect(0:(d-1))
        ens_val = X[randperm(length(X))]
        relinking_aux!(x1, x2, pool, model, ens_val, d, nb)
    end

    function relinking!(x1, x2, pool, model)
           #println(length(pool))
           nb = 6
           d = div(x2.m_length, nb)
           r = rand(0:(d-1))
           same = x1 == x2
           if same
               return
           end
           j = 0
           while x1.dna[(r*nb+1):((r+1)*nb)] == x2.dna[(r*nb+1):((r+1)*nb)]
               #@show j
               r = rand(0:(d-1))
               j+=1
               #@show j
               if j >= nb +3
                   #@show "yes"
                   return
               end
           end

           newS = model.ga.create_entity(x1.dna)
           sameRange = false
           if sum((x1.dna-x2.dna).^2) > (x1.m_length*0.05)
               if r == d -1
                   newS.dna[(r*nb+1):end] = x2.dna[(r*nb+1):end]
                   #newS[(r*5+1):end] = x2[(d*5+1):end]
               else
                   newS.dna[(r*nb+1):((r+1)*nb)] = x2.dna[(r*nb+1):((r+1)*nb)]
                   #newS[(r*5+1):(r+1)*5] = x2[(r*5+1):(r+1)*5]
               end
               push!(pool, newS)
               relinking!(newS, x2, pool, model)
           end
       end

       function relinking2!(x1, x2, pool, model)
              newS = model.ga.create_entity(x1.dna)

              if sum((x1.dna-x2.dna).^2) > (x1.m_length*0.1)
                  for i in 1:x1.m_length
                      newS.dna[i] = x1.dna[i] + (x2.dna[i]-x1.dna[i])/10
                  end
                  push!(pool, newS)
                  relinking!(newS, x2, pool, model)
              end
          end

end


text:
module PrintfFunctions

using Printf

export printf, sprintf

# The following code is copied from
#
#     https://github.com/JuliaIO/Formatting.jl/pull/15/files
#
# and modified to work with Julia v1.x.
#
###########
#
# Implementation of printf() and sprintf() that uses generated functions
# and the existing @printf and @sprintf macros.
#
# Original work by Tim Holy. Minor edits by Daniel Carrera.
#
struct FormatString{S} end

FormatString(str::AbstractString) = FormatString{Symbol(str)}

@generated function Base.print(::Type{FormatString{format}}, args...) where format
    meta = Expr(:meta, :inline)
    fmt = string(format)
    allargs = [:(args[$d]) for d = 1:length(args)]
    quote
        @printf($fmt, $(allargs...))
    end
end
@generated function Base.sprint(::Type{FormatString{format}}, args...) where format
    meta = Expr(:meta, :inline)
    fmt = string(format)
    allargs = [:(args[$d]) for d = 1:length(args)]
    quote
        @sprintf($fmt, $(allargs...))
    end
end

#
# USAGE: printf("%6d  %7.2f", 220/7, 22/7)
#
function printf(s::AbstractString, args...)
    print(FormatString(s), args...)
end
function sprintf(s::AbstractString, args...)
    sprint(FormatString(s), args...)
end
###########

end


text:
## Deprecated

## old initialization algorithms
#
# deprecated at 0.2.7
#

function randseed_initialize!{T<:AbstractFloat}(X::Matrix{T}, centers::Matrix{T})
    Base.depwarn("randseed_initialize! is deprecated. Please use initseeds! instead.",
                 :randseed_initialize!)

    K = size(centers, 2)
    iseeds = initseeds(RandSeedAlg(), X, K)
    copyseeds!(centers, X, iseeds)
end

function kmeanspp_initialize!{T<:AbstractFloat}(X::Matrix{T}, centers::Matrix{T})
    Base.depwarn("kmeanspp_initialize! is deprecated. Please use initseeds! instead.",
                 :kmeanspp_initialize!)

    K = size(centers, 2)
    iseeds = initseeds(KmppAlg(), X, K)
    copyseeds!(centers, X, iseeds)
end

function initial_medoids{T<:AbstractFloat}(costs::Matrix{T}, k::Int)
    Base.depwarn("initial_medoids is deprecated. Please use initseeds_by_costs instead.",
                 :initial_medoids)

    initseeds_by_costs(KmCentralityAlg(), costs, k)
end

function affinity_propagation(x::Matrix; opts...)
    Base.depwarn("affinity_propagation is deprecated. Please use affinityprop instead.",
                 :affinity_propagation)
    affinityprop(x; opts...)
end


text:
using Random, Distributions, MosekTools, JuMP

mutable struct GameNode
    remain::Int
    pwin::Array{Float64}
    plose::Array{Float64}
    score::Array{Int}
    y::Array{Float64}
    children
end

function addNode(score, remain)
    N = length(score)
    node = GameNode(remain, zeros(Float64, N), zeros(Float64, N), score, zeros(Float64, N), nothing)
    addChildren(node)
    return node
end

function initializeGame(N, G)
    return addNode(zeros(Int, N), G)
end

function addChildren(parent)
    remain = parent.remain
    score = parent.score
    N = length(score)
    outcomes = collect(Iterators.product(repeat([[0,1]], N)...))
    pwin = zeros(N)
    plose = zeros(N)
    if remain == 0
        children = nothing
        losers = findall(==(minimum(score)), score)
        winners = findall(==(maximum(score)), score)
        pwin = [a in winners for a in 1:N ]./length(winners)
        plose = [ a in losers for a in 1:N ]./length(losers)
    else
        children = []
        for i in 1:length(outcomes)
            node = addNode(score .+ outcomes[i], remain-1)
            push!(children, node)
            pwin .+= 1/length(outcomes)*node.pwin
            plose .+= 1/length(outcomes)*node.plose
        end
    end
    parent.children = children
    parent.pwin = pwin
    parent.plose = plose
end

function getWeights!(node, wts, g)
    N = length(node.score)
    O = 2^N
    if node.remain== 0
        return 0
    end
    for child in node.children
        wts[g+1] += sum(abs.(node.y .- child.y))/(N*O^(g+1))
        getWeights!(child, wts, g+1)
    end
end
function get_objective(node)
    N = length(node.score)
    obj = 0
    if node.remain == 0
        losers = findall(==(minimum(node.score)), node.score)
        obj += sum((node.plose[l] - node.y[l])^2 - (node.y[l] - 0)^2 for l in losers)
        obj += sum((node.y[i] - 0)^2 for i in 1:N)
        return obj
    else
        for child in node.children
            obj += get_objective(child)/length(node.children)
        end
        return obj
    end
end

#recursive TIC calculation
function calculate_TIC(node::GameNode, pi, stopped = false, yfixed = nothing)
    if stopped
        node.y = copy(yfixed)
    else
        node.y = copy(node.plose)
    end
    # stop recursion
    if node.remain==0
        return 0
    end
    if stopped == false
        winIC, loseIC = calc_incentive(node)
        if minimum(winIC*pi .- loseIC) < 0
            stopped = true
            yfixed = copy(node.plose)
        end
    end
    for child in node.children
        calculate_TIC(child, pi, stopped, yfixed)
    end
end

#check conditional probabilities of winning market prize and getting aid
# when winning vs losing this game
function calc_incentive(node::GameNode)
    N = length(node.score)
    pwin_if_win = zeros(N); pwin_if_lose = zeros(N)
    plose_if_lose = zeros(N); plose_if_win = zeros(N)
    outcomes = collect(Iterators.product(repeat([[0,1]], N)...))
    for (out, child) in zip(outcomes, node.children)
        winners = collect(out .== 1)
        losers = collect(out .== 0)
        pwin_if_win[winners] .+=  child.pwin[winners]/(length(outcomes)/2)
        pwin_if_lose[losers] .+= child.pwin[losers]/(length(outcomes)/2)
        plose_if_lose[losers] .+= child.plose[losers]/(length(outcomes)/2)
        plose_if_win[winners] .+= child.plose[winners]/(length(outcomes)/2)
    end
    win_incentive = pwin_if_win - pwin_if_lose
    lose_incentive = plose_if_lose - plose_if_win
    return win_incentive, lose_incentive
end

function add_IC(node::GameNode, model, g, j, pi)
    N = length(node.score)
    O = 2^N
    outcomes = collect(Iterators.product(repeat([[0,1]], N)...))
    if node.remain == 0
        return 0
    else
        winIC, loseIC = calc_incentive(node)
        for i in 1:N
            record = collect(Iterators.flatten([out[i] for out in outcomes]))
            signs = ones(length(outcomes))
            signs[record.== 1] .= -1
            x = model[:x]
            @constraint(model, sum(1/O*x[g+1, k+O*(j-1), i]*signs[k] for k in 1:O) <= winIC[i]*pi)
        end

        for (child, jnew) in zip(node.children, (O*(j-1)+1):O*j)
            add_IC(child, model, g+1, jnew, pi)
        end
    end
end

function optimalRule(tree::GameNode, pi, incentives = true, weights=false)
    N = length(tree.score)
    G = tree.remain
    losers = findLosers(tree)
    model = Model(Mosek.Optimizer)
    O = 2^N
    #draft rule for all partial histories
    @variable(model, 0<=x[i=0:G, j=1:O^i, k=1:N]<=1)

    #x = model[:x]
    #FAIR constraint
    for k in 1:N
        @constraint(model, x[0, 1, k]==1/N)
    end

    #PROB constraint
    for i in 1:G
        for j in 1:O^i
            @constraint(model, sum(x[i, j, k] for k in 1:N)==1)
        end
    end

    #DC Constraint under equal ability
    for i in 0:(G-1)
        for j in 1:O^i
            for k in 1:N
                @constraint(model, x[i, j, k] == 1/O*sum(x[i+1, s, k] for s in (O*(j-1)+1):O*j))
            end
        end
    end

    if incentives
        add_IC(tree, model, 0, 1, pi)
    end

    @objective(model, Min, sum( sum((1/length(losers[j])-x[G, j, l])^2 - (0-x[G, j, l])^2 for l in losers[j]) +
                                sum( (0-x[G, j, k])^2 for k in 1:N)
                                for j in 1:O^G)/O^G)
    optimize!(model)

    opt_wts = zeros(G)
    for g in 0:(G-1)
        wt = 0
        for j in 1:O^g
            for k in (O*(j-1)+1):O*j
                for n in 1:N
                    wt += abs(value(x[g , j, n]) - value(x[g+1, k, n]))/(N*O^(g+1))
                end
            end
        end
        opt_wts[g+1] = wt
    end

    if weights
        return objective_value(model), opt_wts
    else
        return objective_value(model)
    end
end

function findLosers(node)
    loser_record = []
    loserPass(node, loser_record)
    return loser_record
end

function loserPass(node, loser_record)
    if node.remain == 0
        losers = findall(==(minimum(node.score)), node.score)
        push!(loser_record, losers)
    else
        for child in node.children
            loserPass(child, loser_record)
        end
    end
end

max_stars_repo_name: evanmunro/draft-policy

max_issues_repo_name: evanmunro/draft-policy

max_forks_repo_name: evanmunro/draft-policy


text:
using NMF
using Test

p = 5
n = 8
k = 3

for T in (Float64, Float32)
    Wg = max.(rand(T, p, k) .- T(0.5), zero(T))
    Hg = max.(rand(T, k, n) .- T(0.5), zero(T))
    X = Wg * Hg
    W = Wg .+ rand(T, p, k)*T(0.1)

    NMF.solve!(NMF.CoordinateDescent{T}(α=0.0, maxiter=1000, tol=1e-9), X, W, Hg)

    @test X ≈ W * Hg atol=1e-4
end


text:
import Base: Pair, Tuple, show, ==
import LightGraphs: AbstractEdge, src, dst, reverse

# abstract type AbstractSymbolEdge{Symbol} <: AbstractEdge{Symbol} end

struct SymbolEdge{Symbol} <: AbstractEdge{Symbol}
    src::Symbol
    dst::Symbol
end

SymbolEdge(t::Tuple) = SymbolEdge(t[1], t[2])
SymbolEdge(p::Pair) = SymbolEdge(p.first, p.second)

e = SymbolEdge(:x1,:x2)
e = SymbolEdge((:x1,:x2))

eltype(e::SymbolEdge) = Symbol

# Accessors
src(e::SymbolEdge) = e.src
dst(e::SymbolEdge) = e.dst

# I/O
show(io::IO, e::SymbolEdge) = print(io, "Edge $(e.src) => $(e.dst)")

# Conversions
Pair(e::SymbolEdge) = Pair(src(e), dst(e))
Tuple(e::SymbolEdge) = (src(e), dst(e))

SymbolEdge(e::SymbolEdge) = SimpleEdge(e.src, e.dst)

# Convenience functions
reverse(e::SymbolEdge) = SymbolEdge(dst(e), src(e))
==(e1::SymbolEdge, e2::SymbolEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2))


text:
struct Discontinuity{T<:AbstractFloat, PT<:PathIntegral, QN}
    path::PT
    quadrature::QuadratureRule{T,QN}
end

function Discontinuity(path::PT, quadrature::QuadratureRule{T,QN}) where {T,PT,QN}
    Discontinuity{T,PT,QN}(path, quadrature)
end


text:
struct SimulatorConstraints <: Constraints
    n::Int
    ineq
    w
    slack
end

function simulator_constraints(model, T;
    w = [zeros(model.d) for t = 1:T-1],
    slack = false)

    n = model.n * (T - 1)
    ineq = (1:0)

    return SimulatorConstraints(n, ineq, w, slack)
end

function constraints!(c, Z, con::SimulatorConstraints, model, idx, h, T)
    n = model.n
    nq = model.nq

    # note: x1 and xT constraints are handled as simple bound constraints
    #       e.g., x1 <= x <= x1, xT <= x <= xT

    for t = 1:T-1
        x = view(Z, idx.x[t])
        u = view(Z, idx.u[t])
        x⁺ = view(Z, idx.x[t + 1])

        q1 = view(x, 1:nq)
        q2 = view(x, nq .+ (1:nq))
        q2_next = view(x⁺, 1:nq)
        q3 = view(x⁺, nq .+ (1:nq))

        @show q_sim = (con.slack ? step_slack(q1, q2, u, h, tol = 1.0e-5)[1] : step(q1, q2, u, h, tol = 1.0e-5)[1])
        c[(t-1) * n .+ (1:nq)] = q3 - q_sim
        c[(t-1) * n + nq .+ (1:nq)] = q2_next - q2
    end

    return nothing
end

function constraints_jacobian!(∇c, Z, con::SimulatorConstraints, model, idx, h, T)
    n = model.n
    m = model.m
    nq = model.nq
    # note: x1 and xT constraints are handled as simple bound constraints
    #       e.g., x1 <= x <= x1, xT <= x <= xT

    shift = 0

    In = Diagonal(ones(model.nq))

    for t = 1:T-1
        x = view(Z, idx.x[t])
        u = view(Z, idx.u[t])
        x⁺ = view(Z, idx.x[t + 1])

        q1 = view(x, 1:nq)
        q2 = view(x, nq .+ (1:nq))
        q2_next = view(x⁺, 1:nq)
        q3 = view(x⁺, nq .+ (1:nq))
        _, _, _, Δq1, Δq2, Δu1, _ = (con.slack ? step_slack(q1, q2, u, h, tol = 1.0e-5) : step(q1, q2, u, h, tol = 1.0e-5))

        #c[(t-1) * n .+ (1:nq)] = q3 - step(q1, q2, u, h)[1]
        r_idx = (t-1) * n .+ (1:nq)

        c_idx = idx.x[t + 1][nq .+ (1:nq)]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = vec(In)
        shift += s

        c_idx = idx.x[t][1:nq]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = -1.0 * vec(Δq1)
        shift += s

        c_idx = idx.x[t][nq .+ (1:nq)]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = -1.0 * vec(Δq2)
        shift += s

        c_idx = idx.u[t]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = -1.0 * vec(Δu1)
        shift += s

        # c[(t-1) * n + nq .+ (1:nq)] = q2_next - q2
        r_idx = (t-1) * n + nq .+ (1:nq)

        c_idx = idx.x[t + 1][1:nq]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = vec(In)
        shift += s

        c_idx = idx.x[t][nq .+ (1:nq)]
        s = length(r_idx) * length(c_idx)
        ∇c[shift .+ (1:s)] = -1.0 * Diagonal(In)
        shift += s
    end

    return nothing
end

function constraints_sparsity(con::SimulatorConstraints, model, idx, T;
    shift_row = 0, shift_col = 0)
    n = model.n
    nq = model.nq
    row = []
    col = []

    for t = 1:T-1

        #c[(t-1) * n .+ (1:nq)] = q3 - step(q1, q2, u, h)[1]
        r_idx = (t-1) * n .+ (1:nq)

        c_idx = idx.x[t + 1][nq .+ (1:nq)]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)

        c_idx = idx.x[t][1:nq]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)

        c_idx = idx.x[t][nq .+ (1:nq)]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)

        c_idx = idx.u[t]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)

        # c[(t-1) * n + nq .+ (1:nq)] = q2_next - q2
        r_idx = (t-1) * n + nq .+ (1:nq)

        c_idx = idx.x[t + 1][1:nq]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)

        c_idx = idx.x[t][nq .+ (1:nq)]
        row_col!(row, col, shift_row .+ r_idx, shift_col .+ c_idx)
    end

    return collect(zip(row, col))
end


text:
@windows_only using WinRPM

function copy_libs(src, dst)
    files = readdir(src)

    for i = 1:length(files)
        file = files[i]
        if ismatch(r"\w*?-\w*?(-.)?.(so|dylib|dll)$", file)
            cp("$src/$file", "$dst/$file", follow_symlinks=true, remove_destination=true)
        end
    end
end

function symlink_files(dir, ext)
    cd(dir)
    files = readdir(dir)
    for i = 1:length(files)
        file = files[i]
        filename = file[1:search(file, '.') - 1]
        if ismatch(r"\w*?-\w*?.(so|dylib|dll)$", file)
            run(`ln -sf $filename.$ext $file`)
        end
    end
end

function mkdir_if_necessary(dir)
    if !isdir(dir)
        mkdir(dir)
    end
end

deps = dirname(@__FILE__)
cd(deps)

@osx_only begin
    sfml = "http://www.sfml-dev.org/files/SFML-2.2-osx-clang-universal.tar.gz"
    csfml = "http://www.sfml-dev.org/files/CSFML-2.2-osx-clang-universal.tar.gz"

    if !isfile("sfml.tar.gz")
        println("Downloading SFML...")
        download(sfml, "sfml.tar.gz")
    end
    if !isfile("csfml.tar.gz")
        println("Downloading CSFML...")
        download(csfml, "csfml.tar.gz")
    end

    mkdir_if_necessary("sfml")
    run(`tar -xzf sfml.tar.gz -C sfml --strip-components=1`)

    mkdir_if_necessary("csfml")
    run(`tar -xzf csfml.tar.gz -C csfml --strip-components=1`)

    symlink_files("$deps/csfml/lib", "2.2.0.dylib")

    copy_libs("$deps/sfml/lib", deps)
    copy_libs("$deps/csfml/lib", deps)

    cp("$deps/sfml/extlibs/freetype.framework", "$deps/freetype.framework", remove_destination=true)
    cp("$deps/sfml/extlibs/sndfile.framework", "$deps/sndfile.framework", remove_destination=true)

    cd(deps)
    modules = ["system", "network", "audio", "window", "graphics"]
    for i = 1:length(modules)
        run(`ln -sf libcsfml-$(modules[i]).dylib libcsfml-$(modules[i]).2.2.dylib`)
    end
end

@linux_only begin

    modules = ["system", "network", "audio", "window", "graphics"]

    # check if SFML is installed in the system
    useSystemSFML = false
    systemLibPath = "/usr/lib"
    if isfile("$(systemLibPath)/libcsfml-system.so")
        useSystemSFML = true
    end

    if !useSystemSFML  # get our own, but may not work on some platforms [deps]
        sfml = "http://www.sfml-dev.org/files/SFML-2.2-linux-gcc-$WORD_SIZE-bit.tar.gz"
        csfml = "http://www.sfml-dev.org/files/CSFML-2.2-linux-gcc-$WORD_SIZE-bit.tar.bz2"

        if !isfile("sfml.tar.gz")
            println("Downloading SFML...")
            download(sfml, "sfml.tar.gz")
        end
        if !isfile("csfml.tar.bz2")
            println("Downloading CSFML...")
            download(csfml, "csfml.tar.bz2")
        end

        mkdir_if_necessary("sfml")
        run(`tar -xzf sfml.tar.gz -C sfml --strip-components=1`)

        mkdir_if_necessary("csfml")
        run(`tar -xjf csfml.tar.bz2 -C csfml --strip-components=1`)

        symlink_files("$deps/csfml/lib", "so.2.2.0")

        copy_libs("$deps/sfml/lib", deps)
        copy_libs("$deps/csfml/lib", deps)

        cd(deps)
        for i = 1:length(modules)
            run(`ln -sf libcsfml-$(modules[i]).so libcsfml-$(modules[i]).so.2.2`)
        end
    else
        # use system SFML/CSFML
        for i = 1:length(modules)
            cd(deps)
            run(`ln -sf $(systemLibPath)/libcsfml-$(modules[i]).so libcsfml-$(modules[i]).so.2.2`)
            run(`ln -sf $(systemLibPath)/libsfml-$(modules[i]).so libsfml-$(modules[i]).so`)
        end
    end
end

@windows_only begin
    GCCPath = Pkg.dir("WinRPM","deps","usr","$(Sys.ARCH)-w64-mingw32","sys-root","mingw","bin","gcc.exe")
    if !isfile(GCCPath)
        println("Installing gcc...")
        WinRPM.install("gcc",yes=true)
    end

    cd(deps)
    deps_files = readdir(deps)

    for i = 1:length(deps_files)
        file = deps_files[i]
        if file != "build.jl"
            rm(file, recursive=true)
        end
    end

    println("Downloading SFML and CSFML")
    run(`git clone https://github.com/zyedidia/sfml-binaries.git`)

    mv("sfml-binaries/sfml/sfml-$WORD_SIZE", "sfml")
    mv("sfml-binaries/csfml/csfml-$WORD_SIZE", "csfml")

    copy_libs("$deps/sfml/bin", deps)
    copy_libs("$deps/csfml/bin", deps)
end

cd(joinpath(dirname(@__FILE__),"..","src","c"))
julia_exe = joinpath(JULIA_HOME, Base.julia_exename())
run(`$(julia_exe) createlib.jl`)

cd(deps)

if isfile("sfml")
    rm("sfml", recursive=true)
end
if isfile("csfml")
    rm("csfml", recursive=true)
end

if isfile("libjuliasfml.dylib") || isfile("libjuliasfml.so") || isfile("libjuliasfml.dll")
    println("Successfully built SFML.jl!")
else
    println("Building SFML.jl failed!")
end


text:
#!/usr/bin/env julia
using ProjectEuler
using Statistics
using BenchmarkTools

syms=filter(s->occursin(r"pe\d*",string(s)),names(ProjectEuler))

suite = BenchmarkGroup()

for f in syms
    id = parse(Int,string(f)[3:end])
    suite[id]=@benchmarkable $(eval(f))()
end

tune!(suite)

results = run(suite, verbose = true)

for (k,v) in results
    println(k,": ",median(v))
end

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_name: miRoox/ProjectEuler


text:
using Graphs
using Test

# Test the random graph generation
# We're just running the functions, no testing that the graphs
# actually have the characteristics that they are supposed to have.

n = 10
p = 0.2
let g = erdos_renyi_graph(n,p, is_directed=false)
    # A graph should have n vertices
    @test num_vertices(g) == n
    @test !is_directed(g)
end

p = 1
let g = erdos_renyi_graph(n,p, is_directed=false)
    # When p = 1, the graph should be complete.
    @test num_edges(g) == sum(1:n-1)
end

p = 1
let g = erdos_renyi_graph(n,p, is_directed=true)
    # A graph should have n vertices
    @test num_vertices(g) == n
    @test is_directed(g)
    # When p = 1, the graph should be complete.
    @test num_edges(g) == 2*sum(1:n-1)
end

# Watts-Strogatz small world graphs

# First check that the initial circle is created when beta=0
n = 100
k = 6
beta = 0
let g = watts_strogatz_graph(n, k, beta)
    @test num_edges(g) == n*(k/2)
    @test num_vertices(g) == n
    for i = 1:n
        on = out_neighbors(i,g)
        for j = 1:k/2
            @test (((i+j-1) % n) + 1 in on)
        end
    end
end

# Check that no edges are lost in rewiring:
beta = 0.1
let g = watts_strogatz_graph(n, k, beta)
    @test num_edges(g) == n*(k/2)
    @test num_vertices(g) == n
end


text:
@testset "$ArrayType" for ArrayType in (AxisArray, KeyedArray, NamedDimsArray)

    x = ArrayType([1 2 3; 4 5 6], foo=["a", "b"], bar=[:x, :y, :z])

    @test is_transformable(x)

    @testset "apply" begin

        @testset "OneToOne" begin
            T = FakeOneToOneTransform()
            @test FeatureTransforms.apply(x, T) == ones(2, 3)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 2)
        end

        @testset "OneToMany" begin
            T = FakeOneToManyTransform()
            @test FeatureTransforms.apply(x, T) == ones(2, 6)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2, 2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 4)
        end

        @testset "ManyToOne" begin
            T = FakeManyToOneTransform()
            @test FeatureTransforms.apply(x, T; dims=1) == ones(3)
            @test FeatureTransforms.apply(x, T; dims=1, inds=[1, 2]) == ones(3)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[1, 2]) == ones(2)
        end

        @testset "ManyToMany" begin
            T = FakeManyToManyTransform()
            @test FeatureTransforms.apply(x, T; dims=1) == ones(2, 6)
            @test FeatureTransforms.apply(x, T; inds=[1, 2]) == ones(2, 2)
            @test FeatureTransforms.apply(x, T; dims=2, inds=[2, 3]) == ones(2, 4)
        end

    end

    @testset "apply!" begin
        T = FakeOneToOneTransform()

        _x = copy(x)
        FeatureTransforms.apply!(_x, T)
        @test _x == ones(2, 3)

        # https://github.com/invenia/FeatureTransforms.jl/issues/68
        @test_broken FeatureTransforms.apply!(copy(x), T; inds=[1, 2])
    end

    @testset "apply_append" begin

        @testset "OneToOne" begin

            T = FakeOneToOneTransform()

            @test FeatureTransforms.apply_append(x, T; append_dim=1) == vcat(x, ones(2, 3))
            @test FeatureTransforms.apply_append(x, T; append_dim=2) == hcat(x, ones(2, 3))
            @test FeatureTransforms.apply_append(x, T; append_dim=3) == cat(x, ones(2, 3), dims=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=1, inds=[1], append_dim=1),
                vcat(x, ones(1, 3))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=2, inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 2))
            )
        end

        @testset "OneToMany" begin

            T = FakeOneToManyTransform()

            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; append_dim=2),
                hcat(x, ones(2, 6))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 2))
            )
        end

        @testset "ManyToOne" begin

            T = FakeManyToOneTransform()

            @test FeatureTransforms.apply_append(x, T; dims=1, append_dim=1) == vcat(x, ones(1, 3))
            @test FeatureTransforms.apply_append(x, T; dims=2, append_dim=2) == hcat(x, ones(2))
        end

        @testset "ManyToMany" begin

            T = FakeManyToManyTransform()

            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=1)
            @test_throws DimensionMismatch FeatureTransforms.apply_append(x, T; append_dim=3)

            @test isequal(
                FeatureTransforms.apply_append(x, T; append_dim=2),
                hcat(x, ones(2, 6))
            )

            @test isequal(
                FeatureTransforms.apply_append(x, T; inds=[1, 2], append_dim=2),
                hcat(x, ones(2, 2))
            )
        end

    end

    if ArrayType != AxisArray
        @testset "indexing with dims" begin

            T = FakeOneToOneTransform()

            # apply
            @test FeatureTransforms.apply(x, T; dims=:foo) == ones(2, 3)
            @test FeatureTransforms.apply(x, T; dims=:bar) == ones(2, 3)

            # apply to certain dims/inds
            @test FeatureTransforms.apply(x, T; dims=:bar, inds=Key([:x, :z])) == ones(2, 2)

            # apply_append
            @test isequal(
                FeatureTransforms.apply_append(x, T; dims=:foo, append_dim=:foo),
                cat(x, ones(2, 3); dims=:foo)
            )
        end
    end

end


text:
function DCPP_solve(data, graph)
    """Given data, solve a DRCPP problem with integer programming."""

    ### Constants linked to the problem ###

    # All arcs
    Ar = graph["req_arcs"]
    Anr = graph["noreq_arcs"]
    A = vcat(Ar, Anr)

    # Number of edges/arcs by category
    card_Ar = nrow(Ar)
    card_Anr = nrow(Anr)
    card_A = nrow(A)

    # Traversing cost
    trav_cost = A[:, :trav_cost]

    # Service cost
    serv_cost = A[:, :serv_cost]

    # Demand
    demand = A[:, :demand]

    ### Model definition ###
    model = Model(with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    ### Variables ###
    X_tot = vcat([1 for arc in 1:card_Ar], [0 for arc in 1:card_Anr])
    @variable(model, Y_tot[1:card_A], Int)

    ### Objective ###
    @objective(model, Min, sum(X_tot.*serv_cost+Y_tot.*trav_cost))

    ### Constraint ###
    # Positivity constraint
    @constraint(model, positivity[arc=1:card_A], Y_tot[arc] >= 0)

    # Connectivity constraint
    @constraint(model, connectivity[i=1:data[1,:nodes]],
        sum([X_tot[arc]+Y_tot[arc] for arc in 1:card_A if A[arc,:start_node]==i]) ==
        sum([X_tot[arc]+Y_tot[arc] for arc in 1:card_A if A[arc,:end_node]==i])
        )

    ### Solve problem ###

    start = now()
    optimize!(model)
    stop = now()

    solving_time = stop-start

    ### Format output ###

    tours = []

    cost = sum(X_tot.*serv_cost) + sum(value.(Y_tot).*trav_cost)

    X = Dict(
        "req_arcs" => [1 for arc in 1:card_Ar],
        "req_edges" => [],
        "req_edges_swap" => []
        )

    Y = Dict(
        "req_arcs" => [value.(Y_tot[arc]) for arc in 1:card_Ar],
        "noreq_arcs" => [value.(Y_tot[arc]) for arc in (1+card_Ar):card_A],
        "req_edges" => [],
        "noreq_edges" => [],
        "req_edges_swap" => [],
        "noreq_edges_swap" => []
        )

    tour = Dict(
        "depot" => data[1, :depots][1],
        "cost" => cost,
        "capacity" => sum(X_tot.*demand),
        "X" => X,
        "Y" => Y
        )

    push!(tours, tour)

    return(tours, solving_time)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
"""
This module deals with reading events from the terminal.

Mostly these events will be input from the user, but some of the functionality is about the _terminal_ sending messages as well.

The fundamental piece of functionality that this module offers is being able to a single event from an input stream.

```jldoctest
julia> read(IOBuffer("\e[7;2~"), Event)
Shift+Home
```

If Ansillary does not recognise an event it will return the bytes that it has read _and any remaining bytes in the stream_.

```jldoctest
julia> read(IOBuffer("\e[Q\e[7;2~"), Event)
Ansillary.Inputs.Unknown(UInt8[0x1b, 0x5b, 0x51, 0x1b, 0x5b, 0x37, 0x3b, 0x32, 0x7e])
```

This is done to minimise the risk of the application receiving a seemingly valid event that is actually part of an unknown event, which could lead to highly undesirable consequences. Note that **this is not completely foolproof**, for example Konsole sends `"\x18@sy"` when Super+y is pressed and Ansillary will happily accept that as four separate events.

```jldoctest
julia> let input = IOBuffer("\x18@sy"); [read(input, Event), read(input, Event), read(input, Event), read(input, Event)] end
4-element Array{Ansillary.Inputs.Key,1}:
 Ctrl+x
 @
 s
 y
```

[`Events`](@ref) is an iterator that wraps `read(::Any, ::Event)`, see the file `examples/keylogger.jl` to see how that works. [`EventLoop`](@ref) is an iterator that also sends a [`Tick`](@ref) event every `period`, see the file `examples/snake.jl` or `tests/interactive.jl` to see how that works (also make sure to pay attention to the warnings on that type's documentation).

This module also contains the [`paste`](@ref) function, which allows you to turn on bracketed paste mode. In bracketed paste mode when the user pastes something into the terminal, a [`PasteStart`](@ref) event will be sent, followed by the pasted text, then finally a [`PasteEnd`](@ref) event.

!!! warning

	The functionality of this module will only work properly if the terminal is in raw mode, e.g.

	```julia
	Screen.raw() do
		for key in Events(stdin)
			@show key
			if key == CTRL_C
				break
			end
		end
	end
	```
"""
module Inputs

import ..TERMINAL
import REPL.Terminals

using Dates: Millisecond, Period, Second, now

export Alt,
    Backspace,
    Character,
    Ctrl,
    CTRL_C,
    Delete,
    Down,
    End,
    Esc,
    Event,
    Events,
    EventLoop,
    F,
    Home,
    Insert,
    Left,
    Millisecond,
    Modified,
    Null,
    PageDown,
    PageUp,
    PasteStart,
    PasteEnd,
    Right,
    Second,
    Shift,
    Super,
    Tick,
    Up,
    paste

"""
Permanently set the terminal to bracketed paste mode.

Use [`nopaste!`](@ref) to return the terminal to normal.

!!! note

	Prefer using [`paste`](@ref) where possible as it's very easy to accidently leave the bracketed paste mode enabled using this method.
"""
paste!(terminal = TERMINAL[]) = print(terminal.out_stream, Terminals.CSI, "?2004h")

"""
Turn off bracketed paste mode.
"""
nopaste!(terminal = TERMINAL[]) = print(terminal.out_stream, Terminals.CSI, "?2004l")

"""
Temporarily set bracketed paste mode for the duration of the function.
"""
function paste(f, terminal = TERMINAL[])
    paste!(terminal)
    try
        f()
    finally
        nopaste!(terminal)
    end
end

abstract type Event end

"""
The location of the cursor.

See also: [`Cursor.location`](@ref).
"""
struct Location <: Event
    row::UInt16
    column::UInt16
end

"""
The event sent before text is pasted.

Only sent when bracketed paste mode ([`paste`](@ref)) is enabled.
"""
struct PasteStart <: Event end

"""
The event sent when text has finished being pasted.

Only sent when bracketed paste mode ([`paste`](@ref)) is enabled.
"""
struct PasteEnd <: Event end

"""
An event that is not currently recognised.
"""
struct Unknown <: Event
    data::Vector{UInt8}
end

# When trying to read an event, we want to soak up any extra bytes because it is better to miss an event that was sent quickly than it is to accidentally send an event that the user didn't wish to send.
Unknown(data, input) = Unknown([
    data
    bytesavailable(input) > 0 ? readavailable(input) : []
])

"""
Event that is sent by the event loop every period.
"""
struct Tick <: Event end

abstract type Input <: Event end

abstract type Key <: Input end

Base.show(io::IO, key::Key) = print(io, nameof(typeof(key)))

struct Backspace <: Key end

struct Backtab <: Key end

struct Delete <: Key end

struct Up <: Key end

struct Down <: Key end

struct Left <: Key end

struct Right <: Key end

struct Home <: Key end

struct End <: Key end

struct PageUp <: Key end

struct PageDown <: Key end

struct Insert <: Key end

struct Esc <: Key end

struct Null <: Key end

"""
A function key.
"""
struct F <: Key
    number::UInt8

    F(number) =
        if 1 <= number <= 12
            new(number)
        else
            throw(DomainError(number, "only function keys 12 or smaller are supported"))
        end
end

Base.show(io::IO, f::F) = print(io, 'F', f.number)

"""
A printable character.
"""
struct Character <: Key
    value::Char
end

Base.show(io::IO, character::Character) = print(io, character.value)

abstract type Modifier end

Base.show(io::IO, modifier::Modifier) = print(io, nameof(typeof(modifier)))

struct Super <: Modifier end

struct Ctrl <: Modifier end

struct Alt <: Modifier end

struct Shift <: Modifier end

"""
A key event that has been modified, e.g. Ctrl+c.

A `Modified` event can be created using the `+` operator:

```jldoctest
julia> Ctrl()+'c' == Modified(Character('c'), [Ctrl()])
true

julia> Ctrl()+Super()+Insert() == Modified(Insert(), [Super(), Ctrl()])
true
```

The order of modifiers given to the constructor is not important:

```jldoctest
julia> Modified(Delete(), [Ctrl(), Alt()]) == Modified(Delete(), [Alt(), Ctrl()])
true
```
"""
struct Modified <: Key
    key::Key
    modifiers::Vector{Modifier}

    Modified(key, modifiers) = new(key, unique(modifiers))
end

Base.show(io::IO, modified::Modified) =
    print(io, join(modified.modifiers, '+'), '+', modified.key)

# We actually want the vector in `Modified` to be a set, but actually using `Set` is just annoying. This vector can only have at most 4 elements, so performance not a worry (and probably better than using a `Set` TBH).
function Base.:(==)(left::Modified, right::Modified)
    left.key == right.key && all(modifier in right.modifiers for modifier in left.modifiers)
end

Base.:(+)(l::Modifier, r::Modifier) = [l; r]
Base.:(+)(m::Modifier, k::Key) = Modified(k, [m])
Base.:(+)(m::Modifier, c::Char) = Modified(Character(c), [m])
Base.:(+)(l::Vector{Modifier}, r::Modifier) = [l; r]
Base.:(+)(m::Vector{Modifier}, k::Key) = Modified(k, m)
Base.:(+)(m::Vector{Modifier}, c::Char) = Modified(Character(c), m)

const CTRL_C = Ctrl() + 'c'

# TODO: Rewrite this based on: http://www.inwap.com/pdp10/ansicode.txt

const CTRL_LOWER_RANGE = 0x01:0x1a
const CTRL_LOWER_OFFSET = UInt8('a') - first(CTRL_LOWER_RANGE)
const CTRL_UPPER_RANGE = 0x1c:0x1f
const CTRL_UPPER_OFFSET = UInt8('4') - first(CTRL_UPPER_RANGE)

const DIGITS = UInt8.('0':'9')

const ESC = UInt8(Terminals.CSI[1])
const CSI = UInt8(Terminals.CSI[2])

const FN_INDICATOR_UNKNOWN = UInt8('[')
const FN_RANGE_UNKNOWN = UInt8.('A':'E')
const FN_OFFSET_UNKNOWN = UInt8('A')
const FN_INDICATOR_XTERM = UInt8('O')
const FN_RANGE_XTERM = UInt8.('P':'S')
const FN_OFFSET_XTERM = UInt8('P')

const DIRECTIONS = Dict(
    UInt8('A') => Up,
    UInt8('B') => Down,
    UInt8('C') => Right,
    UInt8('D') => Left,
    UInt8('F') => End,
    UInt8('H') => Home,
)

function modifiers(parameter)
    modmap = parameter - 1
    modifiers = Modifier[]

    if modmap & 0b0001 != 0
        push!(modifiers, Shift())
    end

    if modmap & 0b0010 != 0
        push!(modifiers, Alt())
    end

    if modmap & 0b0100 != 0
        push!(modifiers, Ctrl())
    end

    if modmap & 0b1000 != 0
        push!(modifiers, Super())
    end

    modifiers
end

function Base.read(input::IO, ::Type{Event})
    first = read(input, UInt8)

    # Escape, must check if it's the key or a sequence.
    if first == ESC
        read_esc(input, first)
    elseif first == 0x00
        Null()
    elseif first == 0x7f
        Backspace()
        # Line feeds and carriage returns are indistinguishable from Ctrl+j and Ctrl+m so catch them before checking for Ctrl.
    elseif first == UInt8('\n') || first == UInt8('\r')
        # Enter key sends a carriage return not a line feed, which is surprising to most people.
        Character('\n')
        # Tabs are indistinguishable from Ctrl+i so catch them before checking for Ctrl.
    elseif first == UInt8('\t')
        Character('\t')
    elseif first in CTRL_LOWER_RANGE
        read_ctrl(input, first, CTRL_LOWER_OFFSET)
    elseif first in CTRL_UPPER_RANGE
        read_ctrl(input, first, CTRL_UPPER_OFFSET)
    else
        read_utf8(input, [], first)
    end
end

function read_alt(input, first, second)
    if second in CTRL_LOWER_RANGE
        ctrl = read_ctrl(input, second, CTRL_LOWER_OFFSET)
        Modified(ctrl.key, [ctrl.modifiers; Alt()])
    elseif second in CTRL_UPPER_RANGE
        ctrl = read_ctrl(input, second, CTRL_UPPER_OFFSET)
        Modified(ctrl.key, [ctrl.modifiers; Alt()])
    else
        char = read_utf8(input, [first], second)
        if char isa Character
            Modified(char, [Alt()])
        else
            char
        end
    end
end

function read_csi(input, first, second)
    third = read(input, UInt8)

    if third == FN_INDICATOR_UNKNOWN
        read_fn(input, [first, second, third], FN_RANGE_UNKNOWN, FN_OFFSET_UNKNOWN)
    elseif third in keys(DIRECTIONS)
        return DIRECTIONS[third]()
    elseif third == UInt8('Z')
        return Backtab()
    elseif third == UInt8('M')
        @warn("X10 mouse not implemented")
        Unknown([first, second, third], input)
    elseif third == UInt8('<')
        @warn("XTerm mouse not implemented")
        Unknown([first, second, third], input)
    elseif third in DIGITS
        read_numbered_escape(input, [first, second], third)
    else
        Unknown([first, second, third], input)
    end
end

read_ctrl(_, first, offset) = Modified(Character(Char(first + offset)), [Ctrl()])

function read_esc(input, first)
    if bytesavailable(input) == 0
        return Esc()
    end

    second = read(input, UInt8)

    if second == ESC
        # Double escape is equivalent to single.
        Esc()
        # XTerm-style function keys.
    elseif second == FN_INDICATOR_XTERM
        read_fn(input, [first, second], FN_RANGE_XTERM, FN_OFFSET_XTERM)
    elseif second == CSI
        read_csi(input, first, second)
    else
        read_alt(input, first, second)
    end
end

function read_fn(input, previous, range, offset)
    current = read(input, UInt8)

    if current in range
        F(1 + current - offset)
    else
        # TODO: This can have parameters.
        Unknown([previous; current], input)
    end
end

function read_location(input, previous, buffer, final)
    parameters = parse.(UInt16, split(String(copy(buffer)), ';'))
    if length(parameters) != 2
        return Unknown([previous; buffer; final], input)
    end

    Location(parameters...)
end

function read_modified_direction(input, previous, buffer, final)
    key = DIRECTIONS[final]()
    parameters = parse.(UInt16, split(String(copy(buffer)), ';'))
    Modified(key, modifiers(parameters[end]))
end

function read_numbered_escape(input, previous, final)
    buffer = UInt8[]
    while !(0x40 <= final <= 0x7e)
        push!(buffer, final)
        final = read(input, UInt8)
    end

    if final == UInt8('~')
        read_special_keycode(input, previous, buffer, final)
    elseif final in keys(DIRECTIONS)
        read_modified_direction(input, previous, buffer, final)
    elseif final == UInt8('M')
        @warn("rxvt mouse not implemented")
        Unknown([previous; buffer; final], input)
    elseif final == UInt8('R')
        read_location(input, previous, buffer, final)
    else
        Unknown([previous; buffer; final], input)
    end
end

function read_special_keycode(input, previous, buffer, final)
    parameters = parse.(UInt16, split(String(copy(buffer)), ';'))
    if length(parameters) < 1
        return Unknown([previous; buffer; final], input)
    end

    keycode = parameters[1]
    key = if keycode == 1 || keycode == 7
        Home()
    elseif keycode == 2
        Insert()
    elseif keycode == 3
        Delete()
    elseif keycode == 4 || keycode == 8
        End()
    elseif keycode == 5
        PageUp()
    elseif keycode == 6
        PageDown()
    elseif keycode in 11:15
        F(keycode - 10)
    elseif keycode in 17:21
        F(keycode - 11)
    elseif keycode in 23:24
        F(keycode - 12)
    elseif keycode == 200
        PasteStart()
    elseif keycode == 201
        PasteEnd()
    else
        return Unknown([previous; buffer; final], input)
    end

    if length(parameters) == 1 || parameters[2] == 1
        key
    elseif length(parameters) == 2 && parameters[2] <= 15
        Modified(key, modifiers(parameters[2]))
    else
        Unknown([previous; buffer; final], input)
    end
end

function read_utf8(input, previous, first)
    # 0XXXXXXX is the only byte of a one-byte UTF-8 character.
    # 10XXXXXX is a continuation byte that can't occur as the first byte.
    # 110XXXXX is the first byte of a two-byte UTF-8 character.
    # 1110XXXX is the first byte of a three-byte UTF-8 character.
    # 11110XXX is the first byte of a four-byte UTF-8 character.
    # 11111XXX is an invalid byte that can't occur anywhere.

    ones = leading_ones(first)
    bytes = if ones == 0
        1
    elseif ones == 1 || ones > 4
        return Unknown([previous; first], input)
    else
        ones
    end

    buffer = [first; read(input, bytes - 1)]
    string = String(copy(buffer))

    if isvalid(string) && length(string) == 1
        Character(Base.first(string))
    else
        Unknown([previous; buffer], input)
    end
end

"""
Iterate over input events.

# Examples

```
for key in Inputs.Events(stdin)
	@show key
	if key == Inputs.CTRL_C
		break
	end
end
```
"""
struct Events{I <: IO}
    input::I
end

Events(terminal::Terminals.TTYTerminal) = Events(terminal.in_stream)
Events() = Events(TERMINAL[].in_stream)

Base.iterate(e::Events, _ = nothing) = (read(e.input, Event), nothing)
Base.IteratorSize(::Type{Events{T}}) where {T} = Base.IsInfinite()
Base.IteratorEltype(::Type{Events{T}}) where {T} = Base.HasEltype()
Base.eltype(::Type{Events}) = Input

"""
Iterate over a simple event loop.

The event loop will send a [`Tick`](@ref) event every period (though this isn't guaranteed, for example computationally heavy loop bodies will cause the tick to be sent after the body finishes), and send any received input events between those ticks.

!!! warning

	There is currently a bug in this iterator which makes it "eat" the next byte sent to the stream. This shouldn't be an issue in the common use-case of an event loop, where the process exits when the event loop does.

!!! warning

	The current implementation is incompatible with the following methods because there is a background task constantly reading from the input stream:

	* [`Cursor.location`](@ref).
	* [`Cursor.move!`](@ref) methods that take a [`Cursor.Row`](@ref).
	* [`Screen.size`](@ref).

# Examples

```
for key in Inputs.EventLoop(Inputs.Second(1))
	@show key
	if key == CTRL_C
		break
	end
end
```
"""
struct EventLoop
    channel::Channel{Event}
end

EventLoop(period::Period) = EventLoop(TERMINAL[].in_stream, period)

function EventLoop(input, period::Period)
    channel = Channel{Event}(0)

    reader = @async for event in Events(input)
        put!(channel, event)
    end
    bind(channel, reader)

    ticker = @async begin
        last = now()
        while true
            remaining = last + period - now()
            if remaining > Millisecond(0)
                sleep(remaining)
            end

            put!(channel, Tick())
            last = now()
        end
    end
    bind(channel, ticker)

    EventLoop(channel)
end

Base.iterate(e::EventLoop) = Base.iterate(e, nothing)
Base.iterate(e::EventLoop, _) = (take!(e.channel), nothing)
Base.IteratorSize(::Type{EventLoop}) = Base.IsInfinite()
Base.IteratorEltype(::Type{EventLoop}) = Base.HasEltype()
Base.eltype(::Type{EventLoop}) = Event

end # module


text:
using Cbc
using PowerSimulations
using PowerSystems
using DataStructures
using InfrastructureSystems
const IS = InfrastructureSystems
const PSI = PowerSimulations
const PSY = PowerSystems
Cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer)

abstract type TestOpProblem <: PSI.AbstractOperationsProblem end

system = build_c_sys5_re(; add_reserves = true)
solver = optimizer_with_attributes(Cbc.Optimizer)

devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticBranch),
    :TT => DeviceModel(TapTransformer, StaticBranch),
);
services = Dict{Symbol, ServiceModel}();

template = PSI.OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services);

operations_problem = PSI.OperationsProblem(
    TestOpProblem,
    template,
    system;
    optimizer = solver,
    use_parameters = true,
);

set_services_template!(
    operations_problem,
    Dict(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :Down_Reserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
    ),
)

op_results = solve!(operations_problem)
re_results = PSI.run_economic_dispatch(system; optimizer = solver, use_parameters = true)


text:
function Base.push!(A::AxisVector, item)
    StaticRanges.can_set_last(axes(A, 1)) || throw(MethodError(push!, (A, item)))
    push!(parent(A), item)
    grow_last!(axes(A, 1), 1)
    return A
end

function Base.push!(A::AxisVector, item::Pair)
    axis = axes(A, 1)
    StaticRanges.can_set_last(axis) || throw(MethodError(push!, (A, item)))
    push!(parent(A), last(item))
    push_key!(axis, first(item))
    return A
end

function Base.pushfirst!(A::AxisVector, item)
    can_change_size(A) || throw(MethodError(pushfirst!, (A, item)))
    pushfirst_axis!(axes(A, 1))
    pushfirst!(parent(A), item)
    return A
end

function pushfirst_axis!(axis::Axis, key)
    pushfirst!(keys(axis), key)
    grow_last!(parent(axis), 1)
    return nothing
end

function pushfirst_axis!(axis::AbstractAxis)
    grow_last!(parent(axis), 1)
    return nothing
end

function pushfirst_axis!(axis::Axis)
    grow_first!(keys(axis), 1)
    grow_last!(parent(axis), 1)
    return nothing
end

function Base.pushfirst!(A::AxisVector, item::Pair)
    can_change_size(A) || throw(MethodError(pushfirst!, (A, item)))
    axis = axes(A, 1)
    pushfirst_axis!(axis, first(item))
    pushfirst!(parent(A), last(item))
    return A
end

function Base.mapslices(f, a::AxisArray; dims, kwargs...)
    return reconstruct_reduction(a, Base.mapslices(f, parent(a); dims=dims, kwargs...), dims)
end

macro def_equals(f, X,Y)
    if X === :AxisArray
        if Y === :AxisArray
            esc(quote
                function Base.$f(x::$X, y::$Y)
                    if is_dense_wrapper(x) && is_dense_wrapper(y)
                        return Base.$f(parent(x), parent(y))
                    else
                        for (x_i,y_i) in zip(x,y)
                            Base.$f(x_i, y_i) || return false
                        end
                        return true
                    end
                end
            end)
        else
            esc(quote
                function Base.$f(x::$X, y::$Y)
                    if is_dense_wrapper(x)
                        return Base.$f(parent(x), y)
                    else
                        for (x_i,y_i) in zip(x,y)
                            Base.$f(x_i, y_i) || return false
                        end
                        return true
                    end
                end
            end)
        end
    else
        esc(quote
            function Base.$f(x::$X, y::$Y)
                if is_dense_wrapper(y)
                    return Base.$f(x, parent(y))
                else
                    for (x_i,y_i) in zip(x,y)
                        Base.$f(x_i, y_i) || return false
                    end
                    return true
                end
            end
        end)
    end
end

@def_equals(==, AxisArray, AxisArray)
@def_equals(==, AbstractArray, AxisArray)
@def_equals(==, AxisArray, AbstractArray)
@def_equals(==, AxisArray, AbstractAxis)
@def_equals(==, AbstractAxis, AxisArray)
@def_equals(==, AxisArray, GapRange)
@def_equals(==, GapRange, AxisArray)

@def_equals(isequal, AxisArray, AxisArray)
@def_equals(isequal, AbstractArray, AxisArray)
@def_equals(isequal, AxisArray, AbstractArray)
@def_equals(isequal, AxisArray, AbstractAxis)
@def_equals(isequal, AbstractAxis, AxisArray)

Base.isapprox(a::AxisArray, b::AxisArray; kw...) = isapprox(parent(a), parent(b); kw...)
Base.isapprox(a::AxisArray, b::AbstractArray; kw...) = isapprox(parent(a), b; kw...)
Base.isapprox(a::AbstractArray, b::AxisArray; kw...) = isapprox(a, parent(b); kw...)

Base.copy(A::AxisArray) = AxisArray(copy(parent(A)), map(copy, axes(A)))

function Base.similar(A::AxisArray, ::Type{T}, dims::Tuple{Vararg{Int}}) where {T}
    p = similar(parent(A), T, dims)
    return unsafe_reconstruct(A, p; axes=SimpleAxis.(axes(p)))
end

@inline function Base.similar(A::AxisArray{T}, dims::Tuple{Vararg{Int}}) where {T}
    return similar(A, T, dims)
end
function Base.similar(A::AxisArray, ::Type{T}, dims::Tuple{Vararg{Union{Integer,OneTo}}}) where {T}
    p = similar(parent(A), T, dims)
    c = AxisArrayChecks{CheckedAxisLengths}()
    return AxisArray(p, map((key, axis) -> compose_axis(key, axis, c), dims, axes(p)); checks=c)
end

function Base.similar(a::AxisArray, ::Type{T}, dims::Tuple{Union{Integer, Base.OneTo},Vararg{Union{Integer, Base.OneTo}}}) where {T}
    p = similar(parent(A), T, dims)
    c = AxisArrayChecks{CheckedAxisLengths}()
    return AxisArray(p, map((key, axis) -> compose_axis(key, axis, c), dims, axes(p)); checks=c)
end

function Base.similar(A::AxisArray, ::Type{T}, ks::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}) where {T}
    p = similar(parent(A), T, map(length, ks))
    c = AxisArrayChecks{CheckedAxisLengths}()
    return AxisArray(p, map((key, axis) -> compose_axis(key, axis, c), ks, axes(p)); checks=c)
end

#=
similar(::Type{T}, ks::Tuple{Vararg{AbstractAxis,N} where N}) where T<:AbstractArray
similar(::Type{T}, dims::Tuple{Vararg{Int64,N}} where N) where T<:AbstractArray in Base at abstractarr

ay.jl:675)
26: similar
30: similar
60: similar
61: similar
=#
const DimOrAxes = Union{<:AbstractAxis,Base.DimOrInd}

function Base.similar(a::AbstractArray, ::Type{T}, dims::Tuple{Vararg{DimOrAxes}}) where {T}
    return _similar(a, T, dims)
end

function _similar(a::AxisArray, ::Type{T}, dims::Tuple) where {T}
    p = similar(parent(a), T, map(Base.to_shape, dims))
    axs = map((key, axis) -> compose_axis(key, axis, NoChecks), dims, axes(p))
    return AxisArray{eltype(p),ndims(p),typeof(p),typeof(axs)}(p, axs; checks=NoChecks)
end
function _similar(a, ::Type{T}, dims::Tuple) where {T}
    p = similar(a, T, map(Base.to_shape, dims))
    axs = map((key, axis) -> compose_axis(key, axis, NoChecks), dims, axes(p))
    return AxisArray{eltype(p),ndims(p),typeof(p),typeof(axs)}(p, axs; checks=NoChecks)
end

function Base.similar(A::AxisArray)
    p = similar(parent(A))
    return unsafe_reconstruct(A, p; axes=map(assign_indices, axes(A), axes(p)))
end

function Base.similar(::Type{T}, shape::Tuple{DimOrAxes,Vararg{DimOrAxes}}) where {T<:AbstractArray}
    p = similar(T, Base.to_shape(shape))
    axs = map((key, axis) -> compose_axis(key, axis, NoChecks), shape, axes(p))
    return AxisArray{eltype(p),ndims(p),typeof(p),typeof(axs)}(p, axs; checks=NoChecks)
end
#=
function Base.similar(::Type{T}, ks::Tuple{Vararg{<:AbstractAxis}}) where {T<:AbstractArray}
    p = similar(T, map(length, ks))
    c = AxisArrayChecks{CheckedAxisLengths}()
    return AxisArray(p, map((key, axis) -> compose_axis(key, axis, c), ks, axes(p)); checks=c)
end
=#

for (tf, T, sf, S) in (
    (parent, :AxisArray, parent, :AxisArray),
    (parent, :AxisArray, identity, :AbstractArray),
    (identity, :AbstractArray, parent, :AxisArray))

    @eval function Base.cat(A::$T, B::$S, Cs::AbstractArray...; dims)
        p = cat($tf(A), $sf(B); dims=dims)
        return cat(AxisArray(p, cat_axes(A, B, p, dims)), Cs..., dims=dims)
    end
end

for (tf, T, sf, S) in (
    (parent, :AxisVecOrMat, parent, :AxisVecOrMat),
    (parent, :AxisArray, identity, :VecOrMat),
    (identity, :VecOrMat, parent, :AxisArray))
    @eval function Base.vcat(A::$T, B::$S, Cs::VecOrMat...)
        p = vcat($tf(A), $sf(B))
        axs = vcat_axes(A, B, p)
        return vcat(AxisArray(p, axs; checks=NoChecks), Cs...)
    end

    @eval function Base.hcat(A::$T, B::$S, Cs::VecOrMat...)
        p = hcat($tf(A), $sf(B))
        axs = hcat_axes(A, B, p)
        return hcat(AxisArray(p, axs; checks=NoChecks), Cs...)
    end
end

function Base.hcat(A::AxisArray{T,N}) where {T,N}
    if N === 1
        return AxisArray(hcat(parent(A)), (axes(A, 1), axes(A, 2)))
    else
        return A
    end
end

Base.vcat(A::AxisArray{T,N}) where {T,N} = A

Base.cat(A::AxisArray{T,N}; dims) where {T,N} = A