text:
using Test
using Glob
using Random
using JLD2
using CalibrateEDMF
using CalibrateEDMF.DistributionUtils
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.ReferenceStats
using CalibrateEDMF.TurbulenceConvectionUtils
using CalibrateEDMF.Pipeline
using CalibrateEDMF.Pipeline: get_ref_model_kwargs
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
import TurbulenceConvection
cedmf = pkgdir(CalibrateEDMF)
test_dir = joinpath(cedmf, "test", "Pipeline")
using CalibrateEDMF.HelperFuncs
using Random
include(joinpath(test_dir, "config.jl"))

# Shared simulations
config = get_config()
# Generate reference data
ref_config = config["reference"]
ref_model = ReferenceModel(
    ref_config["y_names"][1],
    ref_config["y_dir"][1],
    ref_config["y_dir"][1],
    ref_config["case_name"][1],
    ref_config["t_start"][1],
    ref_config["t_end"][1],
)
namelist_args = config["scm"]["namelist_args"]
# Generate "true" data
run_reference_SCM(ref_model, run_single_timestep = false)

# Test a range of calibration initializations
@testset "init_calibration" begin

    # Different configurations
    prior_means = [
        Dict("entrainment_factor" => [0.15], "detrainment_factor" => [0.4]),
        nothing,
        nothing,
        Dict("entrainment_factor" => [0.1], "detrainment_factor" => [0.2]),
    ]
    batch_sizes = [nothing, 1, 1, nothing]
    augments = [true, true, true, false]
    l2_regs = [1.0, Dict("entrainment_factor" => [0.0], "detrainment_factor" => [0.1]), nothing, nothing]
    norms = [true, false, true, true]
    dim_scalings = [true, false, true, false]
    tikhonov_modes = ["relative", "absolute", "relative", "relative"]
    modes = ["pmap", "hpc", "hpc", "pmap"]
    validations = [config["reference"], config["reference"], config["reference"], nothing]
    algos = ["Sampler", "Unscented", "SparseInversion", "Inversion"]

    for (prior_mean, batch_size, augment, norm, dim_scaling, tikhonov_mode, mode, validation, algo, l2_reg) in
        zip(prior_means, batch_sizes, augments, norms, dim_scalings, tikhonov_modes, modes, validations, algos, l2_regs)
        config["prior"]["prior_mean"] = prior_mean
        config["reference"]["batch_size"] = batch_size
        config["process"]["augmented"] = augment
        config["regularization"]["l2_reg"] = l2_reg
        config["process"]["algorithm"] = algo
        config["regularization"]["normalize"] = norm
        config["regularization"]["dim_scaling"] = dim_scaling
        config["regularization"]["tikhonov_mode"] = tikhonov_mode
        config["validation"] = validation

        outdir_path = init_calibration(config; config_path = joinpath(test_dir, "config.jl"), mode = mode)

        @test isdir(outdir_path)
        @test isfile(joinpath(outdir_path, "prior.jld2"))
        @test isfile(joinpath(outdir_path, "ekobj_iter_1.jld2"))
        @test isfile(joinpath(outdir_path, "versions_1.txt"))
        !isnothing(batch_size) ? (@test isfile(joinpath(outdir_path, "ref_model_batch.jld2"))) : nothing

        versions = readlines(joinpath(outdir_path, "versions_1.txt"))
        for i in 1:config["process"]["N_ens"]
            @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
        end
        rm(outdir_path, recursive = true)
    end
end

@testset "Pipeline with restart" begin

    ###  Test HPC pipeline
    outdir_path = init_calibration(config; mode = "hpc", config_path = joinpath(test_dir, "config.jl"))

    # Check for output
    @test all(isdir.(config["reference"]["scm_parent_dir"]))
    @test isdir(outdir_path)
    @test isfile(joinpath(outdir_path, "prior.jld2"))
    @test isfile(joinpath(outdir_path, "ekobj_iter_1.jld2"))
    @test isfile(joinpath(outdir_path, "versions_1.txt"))

    ekobj = load(ekobj_path(outdir_path, 1))["ekp"]
    versions = readlines(joinpath(outdir_path, "versions_1.txt"))
    for i in 1:config["process"]["N_ens"]
        @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
    end

    # Run one simulation and perturb results to emulate ensemble
    scm_args = load(scm_init_path(outdir_path, versions[1]))
    batch_indices = scm_args["batch_indices"]
    priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))
    model_evaluator = scm_args["model_evaluator"]
    model_evaluator = precondition(model_evaluator, priors, namelist_args = namelist_args)
    sim_dirs, g_scm_orig, g_scm_pca_orig = run_SCM(model_evaluator, namelist_args = namelist_args)
    for version in versions
        g_scm = g_scm_orig .* (1.0 + rand())
        g_scm_pca = g_scm_pca_orig .* (1.0 + rand())
        jldsave(
            scm_output_path(outdir_path, version);
            sim_dirs,
            g_scm,
            g_scm_pca,
            model_evaluator,
            version,
            batch_indices,
        )
        @test isfile(joinpath(outdir_path, "scm_output_$version.jld2"))
    end

    ek_update(ekobj, priors, 1, config, versions, outdir_path)

    # Test ek_update output
    @test isfile(joinpath(outdir_path, "ekobj_iter_2.jld2"))
    @test isfile(joinpath(outdir_path, "versions_2.txt"))
    versions = readlines(joinpath(outdir_path, "versions_2.txt"))
    for i in 1:config["process"]["N_ens"]
        @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
    end

    restart_calibration(ekobj, priors, 2, config, outdir_path, mode = "hpc")
    @test isfile(joinpath(outdir_path, "ekobj_iter_3.jld2"))
    @test isfile(joinpath(outdir_path, "versions_3.txt"))
    versions = readlines(joinpath(outdir_path, "versions_3.txt"))
    for i in 1:config["process"]["N_ens"]
        @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
    end

    rm(outdir_path, recursive = true)
end

@testset "Pipeline_with_validation" begin

    config["reference"]["batch_size"] = 1
    config["reference"]["n_obs"] = [10]
    config["process"]["algorithm"] = "Unscented"
    config["validation"] = config["reference"]
    config["validation"]["batch_size"] = 1
    config["process"]["augmented"] = true
    config["regularization"]["l2_reg"] = 0.5
    init_calibration(config; mode = "hpc", config_path = joinpath(test_dir, "config.jl"))
    outdir_path_list = glob("results_*_SCM*", config["output"]["outdir_root"])
    outdir_path = outdir_path_list[1]
    versions = readlines(joinpath(outdir_path, "versions_1.txt"))

    # Run one simulation and perturb results to emulate ensemble
    scm_args = load(scm_init_path(outdir_path, versions[1]))
    batch_indices = scm_args["batch_indices"]
    model_evaluator = scm_args["model_evaluator"]
    sim_dirs, g_scm_orig, g_scm_pca_orig = run_SCM(model_evaluator, namelist_args = namelist_args)
    for version in versions
        g_scm = g_scm_orig .* (1.0 + rand())
        g_scm_pca = g_scm_pca_orig .* (1.0 + rand())
        jldsave(
            scm_output_path(outdir_path, version);
            sim_dirs,
            g_scm,
            g_scm_pca,
            model_evaluator,
            version,
            batch_indices,
        )
        jldsave(
            scm_val_output_path(outdir_path, version);
            sim_dirs,
            g_scm,
            g_scm_pca,
            model_evaluator,
            version,
            batch_indices,
        )
        @test isfile(joinpath(outdir_path, "scm_output_$version.jld2"))
        @test isfile(joinpath(outdir_path, "scm_val_output_$version.jld2"))
    end
    @test isfile(joinpath(outdir_path, "ref_model_batch.jld2"))
    ref_model_batch = load(joinpath(outdir_path, "ref_model_batch.jld2"))["ref_model_batch"]
    @test isa(ref_model_batch.ref_models, Vector{ReferenceModel})
    @test isa(ref_model_batch.eval_order, Vector{Int64})
    @test isfile(joinpath(outdir_path, "val_ref_model_batch.jld2"))
    val_ref_model_batch = load(joinpath(outdir_path, "ref_model_batch.jld2"))["ref_model_batch"]
    @test isa(val_ref_model_batch.ref_models, Vector{ReferenceModel})
    @test isa(val_ref_model_batch.eval_order, Vector{Int64})

    priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))
    ekobj = load(ekobj_path(outdir_path, 1))["ekp"]
    ek_update(ekobj, priors, 1, config, versions, outdir_path)

    # Test ek_update output
    @test isfile(joinpath(outdir_path, "ekobj_iter_2.jld2"))
    @test isfile(joinpath(outdir_path, "versions_2.txt"))
    versions = readlines(joinpath(outdir_path, "versions_2.txt"))
    for i in 1:config["process"]["N_ens"]
        @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
    end

    restart_calibration(ekobj, priors, 2, config, outdir_path, mode = "hpc")
    @test isfile(joinpath(outdir_path, "ekobj_iter_3.jld2"))
    @test isfile(joinpath(outdir_path, "versions_3.txt"))
    versions = readlines(joinpath(outdir_path, "versions_3.txt"))
    for i in 1:config["process"]["N_ens"]
        @test isfile(joinpath(outdir_path, "scm_initializer_$(versions[i]).jld2"))
    end
end


text:
using Pkg
Pkg.activate("$(@__DIR__)")
Pkg.update()

using Documenter
using PlateMotionRequests

DOCGEN = """

## API Reference

```@docs
PlateMotionRequests
platemotion
write_platemotion
read_platemotion
```
"""

cp("$(@__DIR__)/../README.md", "$(@__DIR__)/src/index.md", force = true)
@assert isfile("$(@__DIR__)/src/index.md")
open("$(@__DIR__)/src/index.md", write = true, append = true) do io
    write(io, DOCGEN)
end
makedocs(modules = [PlateMotionRequests], sitename = "PlateMotionRequests.jl")

ext: jl


text:
using Regression
using Base.Test

import Regression: Options

srand(654321)

## auxiliary functions

_zoverlap(w_e, w) = countnz((w_e .== 0) & (w .== 0)) / max(countnz(w .== 0), countnz(w_e .== 0))

## prepare for simulation

d = 100
num_zeros = 60
w = randn(d)
iz = sort(randperm(d)[1:num_zeros])
w[iz] = 0.0
@assert countnz(w) == d - num_zeros

n = 5000
X = randn(d, n)
y = X'w + 0.1 * randn(n)

## solve w using proximal gradient descent (i.e. ISTA)

println("    with solver ProxGD")

ret = Regression.solve(linearreg(X, y);
            reg=L1Reg(0.01 * n),
            solver=ProxGD(),
            options=Options(armijo=0.2, ftol=1.0e-8 * n))

w_e = ret.sol
relerr = sumabs2(w_e - w) / sumabs2(w)

@test relerr < 1.0e-3
@test _zoverlap(w_e, w) > 0.95


text:
mysource = Base.read(Base.source_path(), String)
println(Int(ARGS[1] == mysource))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text: Dict{String,Any}("albums" => Dict{String,Any}("songs" => Any[Dict{String,Any}("name" => Dict{String,Any}("value" => "Glory Days","type" => "string"))]))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
export spiketriggeredspectrum, pfcoherence, pfplv, pfppc0, pfppc1, pfppc2, pxcorr

#
# Perform tapered FFT of individual spikes embedded in a continuous signal
#
function spiketriggeredspectrum{T<:Integer,S<:Real}(points::AbstractVector{T}, field::AbstractVector{S},
                                                    window::UnitRange{Int};
                                                    nfft::Int=nextpow2(length(window)),
                                                    freqrange::UnitRange{Int}=1:(nfft >> 1 + 1),
                                                    tapers::Union{AbstractVector, AbstractMatrix}=hanning(length(window)))
    n = length(window)
    nfreq = length(freqrange)
    npoints = length(points)
    nfield = length(field)
    dtype = outputtype(S)

    freqoff = freqrange[1]-1
    winoff = window[1]-1

    fftin = zeros(dtype, nfft)
    sts = Array(Complex{dtype}, length(freqrange), size(tapers, 2), npoints)
    fftout = zeros(Complex{dtype}, nfft >> 1 + 1)

    p = FFTW.Plan(fftin, fftout, 1, FFTW.ESTIMATE, FFTW.NO_TIMELIMIT)

    # Calculate PSD for each point and taper
    @inbounds begin
        for j = 1:npoints
            point = points[j]
            if point-window[1] < 1
                error("Not enough samples before point $j")
            end
            if point+window[end] > nfield
                error("Not enough samples after point $j")
            end

            for i = 1:size(tapers, 2)
                for l = 1:n
                    fftin[l] = tapers[l, i]*field[point+winoff+l]
                end
                FFTW.execute(p.plan, fftin, fftout)
                for l = 1:nfreq
                    sts[l, i, j] = fftout[freqoff+l]
                end
            end
        end
    end

    sts
end

#
# Spike-field coherence
#
# See Fries, P., Roelfsema, P. R., Engel, A. K., König, P., & Singer,
# W. (1997). Synchronization of oscillatory responses in visual cortex
# correlates with perception in interocular rivalry. Proceedings of the
# National Academy of Sciences, 94(23), 12699–12704.
#
# For bias correction, see Grasse, D. W., & Moxon, K. A. (2010).
# Correcting the bias of spike field coherence estimators due to a
# finite number of spikes. Journal of Neurophysiology, 104(1), 548–558.
# doi:10.1152/jn.00610.2009
#
# This is not (quite) the same thing as the spike field coherence
# measure in Jarvis & Mitra, 2001 and implemented in Chronux!
function pfcoherence{T<:Real}(sts::Array{Complex{T},3}; debias::Bool=false)
    npoints = size(sts, 3)

    # Calculate PSD of spike triggered average for each taper
    meanpsd = sum(sts, 3)
    psdsum = Base.sumabs2!(sts, 3)

    # Calculate spike field coherence
    sfc = zeros(T, size(sts, 1))
    @inbounds begin
        for i = 1:size(psdsum, 2), l = 1:size(psdsum, 1)
            c = abs2(meanpsd[l, i])/(npoints*psdsum[l, i])
            if debias
                c = (npoints*c - 1)/(npoints - 1)
            end
            sfc[l] += c
        end
    end
    scale!(sfc, 1.0/size(sts, 2))
end

# Jackknife surrogates of point-field coherence
function pfcoherence_jn{T<:Real}(sts::Array{Complex{T},3}; debias::Bool=false)
    npoints = size(sts, 3)

    # Calculate PSD of spike triggered average for each taper
    meanpsd = sum(sts, 3)
    psdsum = Base.sumabs2!(sts, 3)

    # Calculate spike field coherence
    jn = zeros(T, size(sts, 1), size(sts, 3))
    for ispike = 1:1:size(sts, 3), itaper = 1:size(sts, 2), ifreq = 1:size(sts, 1)
        v = sts[ifreq, itaper, ispike]
        c = abs2(meanpsd[ifreq, itaper] - v)/((npoints-1)*(psdsum[ifreq, itaper] - abs2(v)))
        if debias
            c = ((npoints-1)*c - 1)/(npoints - 2)
        end
        jn[ifreq, ispike] += c
    end
    scale!(jn, 1.0/size(sts, 2))
end

#
# Point-field PLV/PPC
#
# See Vinck, M., Battaglia, F. P., Womelsdorf, T., & Pennartz, C.
# (2012). Improved measures of phase-coupling between spikes and the
# Local Field Potential. Journal of Computational Neuroscience, 33(1),
# 53–75. doi:10.1007/s10827-011-0374-4

# Computes the sum of complex phases over third dimension
function phasesum{T<:Real}(sts::Array{Complex{T},3})
    ps = zeros(Complex{T}, size(sts, 1), size(sts, 2))
    for i = 1:size(sts, 3), j = 1:size(sts, 2), k = 1:size(sts, 1)
        v = sts[k, j, i]
        ps[j, k] += v/abs(v)
    end
end

# Computes the sum of phases for each trial
function phasesumtrial{T<:Real, U<:Integer}(sts::AbstractArray{Complex{T}, 3}, trials::AbstractVector{U})
    if size(sts, 3) != length(trials)
        error("length of trials does not match third dimension of spike-triggered spectrum")
    end
    nutrial = length(unique(trials))
    pstrial = similar(sts, Complex{T}, size(sts, 1), size(sts, 2), nutrial)
    fill!(pstrial, zero(Complex{T}))
    nintrial = zeros(Int, nutrial)

    trial = 1
    lasttrial = trials[1]
    @inbounds begin
        for i = 1:size(sts, 3)
            trial += trials[i] != lasttrial
            lasttrial = trials[i]

            nintrial[trial] += 1
            for j = 1:size(sts, 2), k = 1:size(sts, 1)
                v = sts[k, j, i]
                v /= abs(v)
                pstrial[k, j, trial] += v
            end
        end
    end

    (pstrial, nintrial)
end

# Compute pfdiffsq with preallocated memory for out, pssum, and
# pssqsum. out is assumed to be zeroed; pssum and pssqsum are not.
function _pfdiffsq!{T}(out::Union{Vector{T}, Matrix{T}}, pssum::Array{Complex{T},2}, pssqsum::Array{T,2},
                       pstrial::AbstractArray{Complex{T},3}, trials=1:size(pstrial, 3), idx::Int=1)
    fill!(pssum, zero(Complex{T}))
    fill!(pssqsum, zero(T))
    @inbounds for itrial = trials, itaper = 1:size(pstrial, 2), ifreq = 1:size(pstrial, 1)
        x = pstrial[ifreq, itaper, itrial]
        pssum[ifreq, itaper] += x
        pssqsum[ifreq, itaper] += abs2(x)
    end
    @inbounds for itaper = 1:size(pstrial, 2), ifreq = 1:size(pstrial, 1)
        out[ifreq, idx] += abs2(pssum[ifreq, itaper]) - pssqsum[ifreq, itaper]
    end
    out
end

# Computes the difference between the squared sum of pstrial along the
# third dimension and the sum of squares of pstrial along the third
# dimension. If "true" is provided as the second argument, perform
# jackknifing. If a number is provided, perform bootstrapping.
function pfdiffsq{T<:Real}(pstrial::AbstractArray{Complex{T},3}, jackknife::Bool)
    pssum = Array(Complex{T}, size(pstrial, 1), size(pstrial, 2))
    pssqsum = Array(T, size(pstrial, 1), size(pstrial, 2))
    tval = _pfdiffsq!(zeros(T, size(pstrial, 1)), pssum, pssqsum, pstrial)

    if jackknife
        jnval = zeros(T, size(pstrial, 1), size(pstrial, 3))
        @inbounds begin
            for i = 1:size(pstrial, 3), j = 1:size(pstrial, 2), k = 1:size(pstrial, 1)
                v = pstrial[k, j, i]
                jnval[k, i] += abs2(pssum[k, j] - v) - (pssqsum[k, j] - abs2(v))
            end
        end
    else
        jnval = zeros(T, 0, 0)
    end
    (tval, jnval)
end

function pfdiffsq{T<:Real}(pstrial::AbstractArray{Complex{T},3}, nbootstraps::Int)
    pssum = Array(Complex{T}, size(pstrial, 1), size(pstrial, 2))
    pssqsum = Array(T, size(pstrial, 1), size(pstrial, 2))
    tval = _pfdiffsq!(zeros(T, size(pstrial, 1)), pssum, pssqsum, pstrial)

    bsval = similar(pstrial, T, size(pstrial, 1), nbootstraps)
    fill!(bsval, zero(T))
    @inbounds begin
        trials = Array(Int, size(pstrial, 3))
        for iboot = 1:nbootstraps
            rand!(1:size(pstrial, 3), trials)
            _pfdiffsq!(bsval, pssum, pssqsum, pstrial, trials, iboot)
        end
        retbsval = bsval
    end
    (tval, retbsval)
end

jnvar(jnval, dim) = scale!(var(jnval, dim), (size(jnval, dim) - 1)/size(jnval, dim))

# Computes the PLV
pfplv{T<:Real}(sts::AbstractArray{Complex{T},3}) =
    amean(scale!(phasesum(sts), 1.0/size(sts, 3)), 2)

# Computes the PPC0 statistic from Vinck et al. (2012), which assumes
# that the LFPs surrounding all spikes are independent
pfppc0{T<:Real}(sts::AbstractArray{Complex{T}, 3}) =
    (abs2(phasesum(sts)) - size(sts, 3))/(size(sts, 3)*(size(sts, 3) - 1))

# Computes the PPC1 statistic from Vinck et al. (2012), which accounts
# for statistical dependence of spike-LFP phases within a trial
function pfppc1{T<:Real, U<:Integer}(sts::AbstractArray{Complex{T}, 3}, trials::AbstractVector{U};
                                     estimatevar::Bool=false)
    pstrial, nintrial = phasesumtrial(sts, trials)
    tval, jnval = pfdiffsq(pstrial, estimatevar)

    nintrial2sum = sqsum(nintrial)
    scale!(tval, 1.0/(size(pmtrial, 2) * (abs2(size(sts, 3)) - nintrial2sum)))

    if estimatevar
        scale!(jnval, 1.0/(size(pmtrial, 2) * (abs2(size(sts, 3) - nintrial) -
                                               (nintrial2sum - abs2(nintrial)))))
        (tval, jnvar(jnval, 2))
    else
        tval
    end
end

# Computes the PPC2 statistic from Vinck et al. (2012), which accounts
# for statistical dependence of spike-LFP phases within a trial as well
# as dependence between spike count and phase
function pfppc2{T<:Real, U<:Integer}(sts::AbstractArray{Complex{T}, 3}, trials::AbstractVector{U};
                                     estimatevar::Bool=false, nbootstraps::Int=0)
    estimatevar == false || nbootstraps == 0 ||
        error("cannot both bootstrap and estimate variance by jackknifing")

    pmtrial, nintrial = phasesumtrial(sts, trials)
    @inbounds begin
        for i = 1:size(pmtrial, 3)
            @assert nintrial[i] != 0
            m = 1.0/nintrial[i]
            for j = 1:size(pmtrial, 2), k = 1:size(pmtrial, 1)
                pmtrial[k, j, i] *= m
            end
        end
    end

    tval, jnval = pfdiffsq(pmtrial, nbootstraps != 0 ? nbootstraps : estimatevar)
    normf = 1.0/(size(pmtrial, 2) * size(pmtrial, 3) * (size(pmtrial, 3) - 1))
    scale!(tval, normf)

    if nbootstraps != 0
        (tval, scale!(jnval, normf))
    elseif estimatevar
        scale!(jnval, 1.0/(size(pmtrial, 2) * (size(pmtrial, 3) - 1) * (size(pmtrial, 3) - 2)))
        (tval, jnvar(jnval, 2))
    else
        tval
    end
end

# Binned cross-correlation of event times
# x and y are monotonically increasing vectors of event times, bins is
# a monotonically increasing vector that specifies time bins for the
# cross-correlation function
function pxcorr(x::AbstractVector, y::AbstractVector, bins::AbstractVector)
    counts = zeros(Int, length(bins))
    nx = length(x)
    minbin = minimum(bins)
    maxbin = maximum(bins)
    nbins = length(bins)
    for yspk in y
        xlo = searchsortedfirst(x, yspk+minbin, Base.Sort.Forward)
        xhi = searchsortedlast(x, yspk+maxbin, xlo, nx, Base.Sort.Forward)
        xlo = xlo < 1 ? 1 : xlo
        xhi = xhi > nx ? nx : xhi
        for i = xlo:xhi
            yspkdiff = x[i] - yspk
            bin = searchsortedfirst(bins, yspkdiff)
            bin = bin > nbins ? nbins : bin
            counts[bin] += 1
        end
    end
    counts
end


text:
# providing the necessary arithmetic functions

"""
    +(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d}

Adds the values of `a` and `b` modulo `d` and returns a new ℤd{T, d}
"""
+(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d} = ℤd{T,d}(a.value + b.value)

"""
    -(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d}

Subtracts the value of `b` from `a` modulo `d` and returns a new ℤd{T, d} (always positive).
"""
-(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d} = ℤd{T,d}(a.value - b.value)

"""
    -(x::ℤd{T,d}) where {T <: Integer, d}

Negating x in ℤd yields again x.
"""
-(x::ℤd{T,d}) where {T <: Integer, d} = x

"""
    *(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d}

Multiplies the values of `a` and `b` modulo `d` and returns a new ℤd{T, d}
"""
*(a::ℤd{T, d}, b::ℤd{T, d}) where {T <: Integer, d} = ℤd{T,d}(a.value * b.value)

"""
    *(a::Integer, b::ℤd{T, d}) where {T <: Integer, d}

Enables ℤ × ℤd → ℤd multiplication.
"""
*(a::Integer, b::ℤd{T, d}) where {T <: Integer, d} = ℤd{T,d}(a * b.value)
*(a::ℤd{T, d}, b::Integer) where {T <: Integer, d} = b * a

"""
    convert(::Type{K}, x::ℤd{T, d}) where {K <: Number, T <: Integer, d}

Converts a given ℤd{T, d} to a general number without `mod d`-logic.

This is usefull for many more complex operations and definitions of functions that live outside ℤd.
Also, some external functions are not reimplemented, therefore they can only be used if ℤd{T, d} is treated like a normal number.
If the function-result should again map to ℤd{T, d}, remember to convert the result back to ℤd{T, d}.

!!! info
    If the external function returns floating point values, it's necessary to convert them back to Integers first using something like
    ```julia
    # rounding a float x to the next integer and converting it back to an integer.
    Int(round(x; digits=0))
    ```
"""
convert(::Type{K}, x::ℤd{T, d}) where {K <: Number, T <: Integer, d} = convert(K, x.value)

"""
Promote the `ℤd{T, d}` type over `T`.
"""
promote_rule(::Type{T}, ::Type{ℤd{T, d}}) where {T <: Integer, d} = ℤd{T,d}

"""
Promote the `ℤ₂` type over `Bool`.
"""
promote_rule(::Type{Bool}, ::Type{ℤ₂}) where {d} = ℤ₂

"""
    Int(x::ℤd{T,d}) where {T <: Integer, d}

Converts a given ℤd{T, d} to a general Integer. Used mostly internal by Julia. See `convert(::Type{K}, x::ℤd{T, d})` for more details on conversion.
"""
Int(x::ℤd{T,d}) where {T <: Integer, d} = x.value

"""
    AbstractFloat(x::ℤd{T,d}) where {T <: Integer, d}

Converts a given ℤd{T, d} to a general Float. Used mostly internal by Julia. See `convert(::Type{K}, x::ℤd{T, d})` for more details on conversion.
"""
AbstractFloat(x::ℤd{T,d}) where {T <: Integer, d} = AbstractFloat(x.value)

"""
    show(io::IO, x::ℤd{T, d}) where {T <: Integer, d}

Pretty-printing for ℤd as Integers.
"""
show(io::IO, x::ℤd{T, d}) where {T <: Integer, d} = print(io, x.value)

"""
    trailing_zeros(x::ℤd{T, d}) where {T <: Integer, d}

Internal Julia function needed in some cases, like when using Sets. Passes the call down to the value of x.
"""
trailing_zeros(x::ℤd{T, d}) where {T <: Integer, d} = trailing_zeros(Int(x.value))

# There may be special cases where we want to simplify the operations, these can be appended below
# The special case handled here is for d = 2

"""
    +(x::ℤ₂, y::ℤ₂)

Special case of addition in ℤ₂ using Booleans.
"""
+(x::ℤ₂, y::ℤ₂) = ℤ₂(x.value != y.value)

"""
    -(x::ℤ₂, y::ℤ₂)

Special case of subtraction in ℤ₂ using Booleans.
"""
-(x::ℤ₂, y::ℤ₂) = ℤ₂(x.value != y.value)

"""
    +(x::ℤ₂, y::ℤ₂)

Special case of multiplication in ℤ₂ using Booleans.
"""
*(x::ℤ₂, y::ℤ₂) = ℤ₂(x.value * y.value)

"""
    show(io::IO, x::ℤ₂)

Special case of pretty-printing in ℤ₂ using Booleans.
"""
show(io::IO, x::ℤ₂) = print(io, Int8(x.value))

"""
    inv(A::Matrix{ℤd{T, d}})

Invert a matrix `A` in ℤd by multiplying the inverse in ℝ with the determinant.
"""
inv(A::Matrix{ℤd{T, d}}) where {T <: Integer, d} = Matrix{ℤd{T, d}}(Int.(round.(inv(Matrix{Int}(A)) * Int(round(det(A); digits=0)); digits=0)))


text:
#!/bin/julia

using DelimitedFiles

#=
FAZER

- Criar funçao para gerar amostras com algumas palavras
- Criar um módulo de limpeza de texto chamado `Limpeza.jl` (copiar as Regex do projeto tfidf-simples)
- Criar ferramenta que descobre qual o idioma de uma frase com base na lista de bigramas coletados.

=#

home = pwd()
homecorpus = home * "/" * "corpus"

function get_linguas()::Array{String,1}
	lis = readdir(homecorpus)
	return lis
end

function get_textos(lingua::String)::Array{String,1}
	return readdir(homecorpus * "/" * lingua)
end

function ler(lingua::String, texto::String)::String
	arquivo = homecorpus * "/" * lingua * "/" * texto
	return read(arquivo, String)
end

function lerlinhas(lingua::String, texto::String)::Array{String,1}
	arquivo = homecorpus * "/" * lingua * "/" * texto
	return readlines(arquivo)
end

function pushdict!(dict::Dict, item::Any)
	if (item in keys(dict))
		dict[item] += 1
	else
		push!(dict, item => 1)
	end
end

function bigramas!(tokens::Dict, texto::String)::Dict{Any,Any}

	words = split(lowercase(texto))

	for w in words

		#   | margem |        meio        | margem |
		#   |    1   |        0..n        |   1    |
		#
		# a:     [[a_,1,1]]
		# de:    [[de,1,1]]
		# com:   [[co,1,0], [om,0,1]]
		# para:  [[pa,1,0], [ar,0,0], [ra,0,1]]
		# arara: [[ar,1,0], [ra,0,0], [ar,0,0], [ra,0,1]]
		chars = split(w, "")
		len = length(chars)
		margin = 1 # margen direita e esquerda

		if len == 1
			token = Any[]
			push!(token, chars[1])
			push!(token, nothing)
			push!(token, 1)
			push!(token, 1)
			pushdict!(tokens, token)
			continue
		end

		i = 1
		while i < len
			token = Any[]
			push!(token, chars[i])
			push!(token, chars[i+1])
			if i <= margin # margem direita
				push!(token, i)
			else
				push!(token, 0)
			end
			if i >= len - margin # margem esquerda
				push!(token, len - i)
			else
				push!(token, 0)
			end
			pushdict!(tokens, token)
			i += 1
		end
	end

	return tokens
end

function main()
	linguas = Dict{String, Dict}(map(x -> Pair(x, Dict()), get_linguas()))
	for lingua in keys(linguas)
		tokens = Dict()
		textos = get_textos(lingua)
		for texto in textos
			bigramas!(tokens, ler(lingua, texto))
		end
		linguas[lingua] = tokens
		println(lingua, ": ", length(tokens))
	end

	# CSV: [lingua, countsum]
	L = map(k -> [k, sum(values(linguas[k])) ], collect(keys(linguas)))
	writedlm( "linguas.csv",  L, ',')

	# CSV: [char1, char2, margin1, margin2, df, idf]
	d = Dict()
	foreach(x -> foreach(y -> pushdict!(d, y), keys(x[2])), collect(linguas))
	B = map(k -> vcat(k[1], k[2], log(length(linguas)/k[2])), collect(d))
	writedlm("bigramas.csv",  B, ',')

end

main()

max_stars_repo_name: fabiolimace/lincom

max_issues_repo_name: fabiolimace/lincom

max_forks_repo_name: fabiolimace/lincom


text:
struct DecodedPosterior
    pa::String
    po::String
    start::Int
    finish::Int
end

function Base.show(io::IO, x::DecodedPosterior)
    println(io,x.pa)
    println(io,x.po)
    print_with_color(:cyan,io,x.fds,"\n")
    print(io,x.fdp)
end

function check_assignment(P,verbose,N)

    count = 0
    L = length(P)
    if verbose == true
	println("Let us check the assignment...")
    end
    for i in L:-1:2
        maxP = -Inf
	idxn = []
	idxx = []
	for x in 0:1, n in 0:N+1
	   if P[i][x,n] > maxP
	      maxP = P[i][x,n]
	      idxx = x
	      idxn = n
	   end
	end
	if maxP == -Inf
	   println("Problem with node $i")
	end
	n =idxn
	x =idxx
        maxP = -Inf
        idxn = []
        idxx = []
        for xj in 0:1, nj in 0:N+1
	    if P[i-1][xj,nj] > maxP
	       maxP = P[i-1][xj,nj]
	       idxx = xj
	       idxn = nj
	    end
	 end
	 nj = idxn
	 xj = idxx
	 sat = x == 1 ? (n > nj) : (n == nj || n == N+1)
        if sat == false && verbose == true
	   println("- $i → ($x, $n) $i-1 → ($xj, $nj)")
        end
        count += 1 - convert(Int, sat)
    end

    if count == 0
	if verbose == true
	    println("The subsequence satisfies the constraints")
	end
        return true
    else
	if verbose == true
	    println("There are $count short-range constraints not satisfied")
	end
	return false
    end

end

#amino = ['A', 'C', 'D', 'E', 'F', 'G', 'H','I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y', '-']
#nbase = ['A', 'U', 'C', 'G', '-']

aminoalphabet = Dict{String,Int64}()
nbasealphabet = Dict{String,Int64}()
aminoalphabet = Dict('A'=> 1, 'B'=> 21, 'C' =>2, 'D'=> 3, 'E'=> 4, 'F'=> 5, 'G'=> 6, 'H'=> 7, 'I'=> 8, 'J'=> 21,'K'=> 9, 'L'=> 10, 'M'=> 11, 'N'=>12, 'O'=>21, 'P'=> 13, 'Q'=>14, 'R'=>15, 'S'=>16, 'T'=>17, 'U'=>21, 'V'=>18,'W'=>19, 'X'=>21, 'Y'=>20, 'Z'=> 21, '-'=> 21)
nbasealphabet = Dict('A'=> 1, 'U'=> 2, 'C'=> 3, 'G'=> 4, '-'=> 5 ,'T'=> 2, 'N'=> 5, 'R'=> 5, 'X' => 5, 'V'=> 5, 'H'=>5, 'D'=>5, 'B'=>5, 'M'=>5, 'W'=>5, 'S'=>5, 'Y'=>5, 'K'=>5)

function letter2num(c::Union{Char,UInt8}, ctype::Symbol)
    if ctype == :amino
	return aminoalphabet[c]
    end
    if ctype == :nbase
	return nbasealphabet[c]
    end
end

function compute_cost_function(J::Array{Float64,4}, h::Array{Float64,2}, seqins::String, L::Int, ctype::Symbol, λo::Vector{Float64}, λe::Vector{Float64}, μext::Float64, μint::Float64)

   seq = ""
   N = length(seqins)
   for i = 1:N
      if isuppercase(seqins[i]) || seqins[i] == '-'
	 seq *= seqins[i]
      end
   end
   #println(seqins)
   #println(seq)
   en = compute_potts_en(J,h,seq,L,ctype)
   extend = false
   idx = 0
   first = true
   nopen = 0
   next = 0
   for i = 1:N
      if isuppercase(seqins[i]) || seqins[i] == '-'
	 idx += 1
	 extend = false
      elseif islowercase(seqins[i]) & extend == false
	 extend = true
	 nopen += 1
	 en += λo[idx+1]
	 #println("open ", idx+1, " ", λo[idx+1])
      elseif islowercase(seqins[i]) & extend == true
	 en += λe[idx+1]
	 next += 1
	 #println("ext ", idx+1, " ", λe[idx+1])
      end
   end
   ngaps_int = 0
   ngaps_ext = 0
   first = true
   for i = 1:L
      if seq[i] == '-'
	 if first
	    ngaps_ext += 1
	 else
	    ngaps_int += 1
	 end
      end
      if seq[i] != '-'
	 first = false
      end
   end
   first = true
   for i = L:-1:1
      if seq[i] == '-'
	 if first
	    ngaps_ext += 1
	    ngaps_int -= 1
	 end
      end
      if seq[i] != '-'
	 first = false
      end
   end
   #println(ngaps_ext, " gaps esterni e ", ngaps_int, " interni")
   #println(nopen, " open insertions plus ", next, " extended")
   en += μext * ngaps_ext + μint * ngaps_int
   return en
end

function compute_potts_en(J::Array{Float64,4}, h::Array{Float64,2}, seq::String, L::Int,ctype::Symbol)

    en = 0
    for i in 1:L
        Aᵢ = letter2num(seq[i],ctype)
        en += -h[Aᵢ,i]
        for j in i+1:L
            Aⱼ = letter2num(seq[j],ctype)

            en += -J[Aᵢ,Aⱼ,i,j]
        end
    end
    return en
end

function count_gaps_ins(v::String)

    N = length(v)
    Ng = 0
	Ni = 0
	Nb = 0
    types = zeros(N)
    for i = 1:N
        if isuppercase(v[i])
            types[i] = 1
        end
        if islowercase(v[i])
            types[i] = -1
        end
        if v[i] == '-'
            types[i] = 0
        end
    end
    for i = 1:N-1
        if v[i] =='-'
            Ng += 1
        elseif islowercase(v[i])
            Ni += 1
        end
        if types[i] != types[i+1]
            Nb += 1
        end
    end
    if v[N] == '-'
        Ng += 1
    end
    if islowercase(v[N])
        Ni += 1
    end

    return Ng, Ni, Nb
end

function hammingdist(seqaux, seqtrue)
    length(seqaux) == length(seqtrue) || error("seqs of different lengths")
    ctr = 0
	moreg = 0
	lessg = 0
	diffm = 0
    for i in eachindex(seqaux)
        if seqaux[i] !== seqtrue[i]
            ctr += 1
        end
    	if seqaux[i] == '-' && seqtrue[i] != '-'
	    	moreg += 1
		end
		if seqaux[i] != '-' && seqtrue[i] == '-'
	    #display(["lg" seqaux[i] seqtrue[i]])
	    	lessg += 1
		end
		if seqaux[i] != '-' && seqtrue[i] != '-' && seqaux[i] != seqtrue[i]
	    	diffm +=1
		end
    end
    return ctr, moreg, lessg, diffm
end

function read_parameters(
    filename::String,
    q::Int,
    L::Int;
    gap::Int = 0,
    typel::Symbol = :bm,
)

    if typel == :plm
        @info "Assuming J a b i j and h a i format"
    else
        @info "Assuming J i j a b and h i a format"
    end
    @info "Output tersors: J[a b i j] and h[a i]"
    @info "Gap in input file $gap now in $q"
    J = zeros(q, q, L, L)
    h = zeros(q, L)

    if gap == q
        offset = 0
    else
        offset = 1
    end
    open(filename) do file
        for ln in eachline(file)
            line = split(ln, ' ')
            if occursin('J', ln)
                if typel == :bm
                    i = parse(Int64, line[2]) + offset
                    j = parse(Int64, line[3]) + offset
                    a = gap == q ? parse(Int64, line[4]) + offset :
                        parse(Int64, line[4])
                    b = gap == q ? parse(Int64, line[5]) + offset :
                        parse(Int64, line[5])
                else
                    i = parse(Int64, line[4]) + offset
                    j = parse(Int64, line[5]) + offset
                    a = gap == q ? parse(Int64, line[2]) + offset :
                        parse(Int64, line[2])
                    b = gap == q ? parse(Int64, line[3]) + offset :
                        parse(Int64, line[3])
                end
                if a == gap && gap == 0
                    a = q
                end
                if b == gap && gap == 0
                    b = q
                end
                J[a, b, i, j] = parse(Float64, line[6])
                J[b, a, j, i] = parse(Float64, line[6])
            end
            if occursin('h', ln)
                if typel == :bm
                    i = parse(Int64, line[2]) + offset
                    a = gap == q ? parse(Int64, line[3]) + offset :
                        parse(Int64, line[3])
                else
                    i = parse(Int64, line[3]) + offset
                    a = gap == q ? parse(Int64, line[2]) + offset :
                        parse(Int64, line[2])
                end
                if a == gap && gap == 0
                    a = q
                end
                h[a, i] = parse(Float64, line[4])
            end
        end
    end
    return J, h
end

function decodeposterior(P,strseq)

    N = length(strseq)
    L = length(P)
    pa = ""
    for i in 1:L
        maxP = -Inf
	idxx = []
	idxn = []
        for x in 0:1, n in 0:N+1
	    if P[i][x,n] > maxP
	       maxP = P[i][x,n]
	       idxx = x
	       idxn = n
	    end
        end
	if maxP == -Inf
	   println("Problem with $i")
	end
	n = idxn
	x = idxx
	if x == 0
	    pa *= '-'
        else
	    pa *= strseq[n]
        end
    end
    po = ""
    i = 1
    start = false
    nold = 0
	f = 0
	l = 0
    while i <= L
        maxP = -Inf
	idxx = []
	idxn = []
	for x in 0:1, n in 0:N+1
	   if P[i][x,n] > maxP
	      maxP = P[i][x,n]
	      idxx = x
	      idxn = n
	   end
	end
        n = idxn
		x = idxx
        if x == 1 && start == true
            delta = n - nold - 1
            if delta > 0
                for d in 1:delta
                    po *= lowercase(strseq[nold+d])
                end
            end
            po *= strseq[n]
            nold = n
            l = n
        end
        if x == 1 && start == false
            start = true
            po *= strseq[n]
            nold = n
            f = n
        end
        if x == 0
            po *= "-"
        end
        i += 1
    end
    DecodedPosterior(pa,po,f,l)
end


text:
#using SparseArrays

function mehortra_guess2( nlp, pars, timer, x::Vector, a::Vector, J, g::Vector )
    s = deepcopy(a);
    m = length(s)

    if true
        if pars.output_level >= 2
          println("estimating intial y")
      end
      start_advanced_timer(timer, "INIT/estimate_y_tilde")
      y = estimate_y_tilde( J, g, pars )
      if isbad(y)
        warning("y is bad!!!!!")
      end
      pause_advanced_timer(timer, "INIT/estimate_y_tilde")
    else
      y = ones(m)
    end

    return s, y
end

function mehortra_guess( nlp, pars, timer, x::Vector, a::Vector, J, g::Vector )
    s = deepcopy(a);
    m = length(s)

    if true
        if pars.output_level >= 2
      println("estimating intial y")
  end
      start_advanced_timer(timer, "INIT/estimate_y_tilde")
      y = estimate_y_tilde( J, g, pars )
      if isbad(y)
        warning("y is bad!!!!!")
      end
      pause_advanced_timer(timer, "INIT/estimate_y_tilde")
    else
      y = ones(m)
    end

    return mehortra_guarding( nlp, pars, timer, x, y, s, a, J, g )
end

function mehortra_guarding( nlp, pars, timer, x::Vector, y_tilde::Vector, s_tilde::Vector, a::Vector, J, g::Vector )
    start_advanced_timer(timer, "INIT/mehortra_guarding")

    bounds = true
    if bounds
      ais = cons_indicies(nlp)
      bis = bound_indicies(nlp)
      s_tilde[bis] = a[bis]
      #@show ais, bis
    else
      ais = 1:(length(cons_indicies(nlp)) + length(bound_indicies(nlp)))
      bis = bound_indicies(nlp)
    end

    if isbad(s_tilde)
      throw(Eval_NaN_error(getbad(s_tilde),x,"s"))
    end

    if isbad(y_tilde)
      throw(Eval_NaN_error(getbad(y_tilde),x,"y"))
    end
    d_s = max(-2.0 * minimum(s_tilde[ais]), 0.0) + LinearAlgebra.norm(g - J' * y_tilde,Inf) / (1.0 + LinearAlgebra.norm(y_tilde,Inf))
    #@show LinearAlgebra.norm(g - J' * y_tilde,Inf) / (1.0 + LinearAlgebra.norm(y_tilde,Inf))
    d_y = max(-2.0 * minimum(y_tilde), 0.0)

    s_tilde[ais] = s_tilde[ais] .+ d_s .+ 1e-8
    y_tilde = y_tilde .+ d_y
    if isbad(y_tilde)
      throw(Eval_NaN_error(getbad(y_tilde),x,"y"))
    end

    d_y_tilde = d_y .+ 0.5 * dot(s_tilde, y_tilde) / sum(s_tilde)
    y_tilde = y_tilde .+ d_y_tilde
    y_tilde = min.(pars.init.dual_max, max.(y_tilde, pars.init.dual_min))

    d_s_tilde = d_s .+ 0.5 * dot(s_tilde, y_tilde) / sum(y_tilde)
    s_tilde[ais] = s_tilde[ais] .+ d_s_tilde

    #@show d_s, d_y, LinearAlgebra.norm(g)
    #s_new, y = mehortra_guarding( deepcopy(s), deepcopy(y), threshold )
    if isbad(y_tilde)
      throw(Eval_NaN_error(getbad(y_tilde),x,"y"))
    end
    if isbad(s_tilde)
      throw(Eval_NaN_error(getbad(s_tilde),x,"s"))
    end

    #@show y, J

    ifree = _i_not_fixed(nlp.nlp)
    uvar = nlp.nlp.meta.uvar[ifree]
    lvar =  nlp.nlp.meta.lvar[ifree]

    @assert(length(bis) == sum(uvar .< Inf) + sum(lvar .> -Inf))
    #lb()

    for i in bound_indicies(nlp)
        if (s_tilde[i] <= 0.0)
            println("Error:")
            @show i, s_tilde[i], a[i] #lvar[i], uvar[i]
        end
        if bounds
          @assert(s_tilde[i] == a[i])
        end
    end
    pause_advanced_timer(timer, "INIT/mehortra_guarding")

    return s_tilde, y_tilde
end

function estimate_y_tilde( J::SparseMatrixCSC{Float64,Int64}, g::Array{Float64,1}, pars::Class_parameters )
    try
      if pars.output_level >= 2
          println("estimate y_tilde ...")
          @show densest_col(J)
          @show densest_row(J)
      end

      #@time H = Symmetric(J * J' + LinearAlgebra.norm(J,Inf) * 1e-4 * SparseArrays.speye( size(J,1) ))
      n = size(J,2); m = size(J,1);

      lambda = 1e-4
      if pars.output_level >= 2
          @show lambda
      end
      if false
        #@time H = [SparseArrays.speye(n) -J'; J lambda * SparseArrays.speye(m)]
        @time H = [SparseMatrixCSC{Float64}(LinearAlgebra.I, n, n) -J'; J lambda * SparseMatrixCSC{Float64}(LinearAlgebra.I, m, m)]
        @time M = lufact( H );
        rhs = [-g; zeros(m)];
        sol = M \ rhs
        y = sol[(n+1):end]
      else # cholesky factor
        scaling = 1.0 #LinearAlgebra.norm(J,Inf)^2
        #H = lambda * SparseArrays.speye(n) + J' * J / scaling;
	H = lambda * SparseMatrixCSC{Float64}(LinearAlgebra.I, n, n) + J' * J / scaling;
        #M = cholfact( H );
	M = cholesky( H );
        dx = scaling * (M \ -g)
        y = -J * dx
      end

      if pars.output_level >= 2
          println("linear system solved in init")
      end
      return y
    catch (e)
      println("error in estimate_y_tilde")
      @show e
      return ones(size(J,1))
    end
end

function mehortra_guarding( s_tilde::Array{Float64,1}, y_tilde::Array{Float64,1}, threshold::Float64 )
    d_s = max(-2.0 * minimum(s_tilde), threshold)
    d_y = max(-2.0 * minimum(y_tilde), threshold)

    d_s_tilde = d_s .+ 0.5 * dot(s_tilde .+ d_s, y_tilde .+ d_y) / sum(y_tilde .+ d_y)
    d_y_tilde = d_y .+ 0.5 * dot(s_tilde .+ d_s, y_tilde .+ d_y) / sum(s_tilde .+ d_s)

    s = s_tilde .+ d_s_tilde
    y = y_tilde .+ d_y_tilde

    return s, y
end
#=
kkt_solver = pick_KKT_solver(pars);
form_system!(kkt_solver, iter, timer)
fact_succeed, inertia_num_fac, new_delta = ipopt_strategy!(iter, kkt_solver, pars, timer)
kkt_associate_rhs!(kkt_solver, iter, Reduct_affine(), timer)
compute_direction!(kkt_solver, timer)
=#
function change_mu!(iter::Class_iterate,new_mu::Float64, pars::Class_parameters)
    iter.point.mu = new_mu
    #iter.primal_residual_intial = (get_cons(iter) - iter.point.s) / new_mu
    #iter.primal_residual_intial = (iter.point.s - get_cons(iter)) / new_mu

    y_c = new_mu ./ iter.point.s
    buffer = 2.0
    y = iter.point.y
    y = min.( y_c / (pars.ls.comp_feas * buffer), max.(y, pars.ls.comp_feas * y_c * buffer)) # project onto complementarity constraints

    iter.point.y = y
end

#=
GLOBAL_COUNT = 0
function KNITRO_guess(iter::Class_iterate, dir::Class_point, pars::Class_parameters)
    global GLOBAL_COUNT
    GLOBAL_COUNT += 1
    if GLOBAL_COUNT == 100
      y = iter.point.y
      s = iter.point.s
      y_temp = y + dir.y
      s_temp = s + dir.s
      y_est = y_temp - 2.0 * min(0.0, minimum(y))
      s_est = s_temp - 2.0 * min(0.0, minimum(s))

      mu = iter.point.mu
      mu_est = LinearAlgebra.norm(y,Inf) * LinearAlgebra.norm(get_primal_res(iter),Inf)
      #0.5 * (iter)+ Statistics.mean(-y_guess .* get_primal_res(iter))
      #mu_est = LinearAlgebra.norm(dir.x,1) #min(mu * 2.0, max(mu/2.0, mu_est)) #+ get_cons(iter) .* iter.point.y
      suggested_mu = min(100.0 * mu,max(mu_est,0.01 * mu)) #
      #@show suggested_mu
      change_mu!(iter,suggested_mu,pars)
    end
end=#


text:
## Uniform ##

struct TuringUniform{T} <: ContinuousUnivariateDistribution
    a::T
    b::T
end
TuringUniform() = TuringUniform(0.0, 1.0)
function TuringUniform(a::Int, b::Int)
    return TuringUniform{Float64}(Float64(a), Float64(b))
end
function TuringUniform(a::Real, b::Real)
    T = promote_type(typeof(a), typeof(b))
    return TuringUniform{T}(T(a), T(b))
end
Distributions.logpdf(d::TuringUniform, x::Real) = uniformlogpdf(d.a, d.b, x)
Distributions.logpdf(d::TuringUniform, x::AbstractArray) = uniformlogpdf.(d.a, d.b, x)
Base.minimum(d::TuringUniform) = d.a
Base.maximum(d::TuringUniform) = d.b

Distributions.Uniform(a::TrackedReal, b::Real) = TuringUniform{TrackedReal}(a, b)
Distributions.Uniform(a::Real, b::TrackedReal) = TuringUniform{TrackedReal}(a, b)
Distributions.Uniform(a::TrackedReal, b::TrackedReal) = TuringUniform{TrackedReal}(a, b)
Distributions.logpdf(d::Uniform, x::TrackedReal) = uniformlogpdf(d.a, d.b, x)
Distributions.logpdf(d::Uniform, x::TrackedArray) = uniformlogpdf.(d.a, d.b, x)
function uniformlogpdf(a, b, x)
    c = -log(b - a)
    if a <= x <= b
        return c
    else
        return oftype(c, -Inf)
    end
end
uniformlogpdf(a::Real, b::Real, x::TrackedReal) = track(uniformlogpdf, a, b, x)
uniformlogpdf(a::TrackedReal, b::TrackedReal, x::Real) = track(uniformlogpdf, a, b, x)
uniformlogpdf(a::TrackedReal, b::TrackedReal, x::TrackedReal) = track(uniformlogpdf, a, b, x)
@grad function uniformlogpdf(a, b, x)
    diff = data(b) - data(a)
    T = typeof(diff)
    if a <= data(x) <= b && a < b
        l = -log(diff)
        da = 1/diff^2
        return l, Δ -> (da * Δ, -da * Δ, zero(T) * Δ)
    else
        n = T(NaN)
        return n, Δ -> (n, n, n)
    end
end
ZygoteRules.@adjoint function uniformlogpdf(a, b, x)
    diff = b - a
    T = typeof(diff)
    if a <= x <= b && a < b
        l = -log(diff)
        da = 1/diff^2
        return l, Δ -> (da * Δ, -da * Δ, zero(T) * Δ)
    else
        n = T(NaN)
        return n, Δ -> (n, n, n)
    end
end
ZygoteRules.@adjoint function Distributions.Uniform(args...)
    return ZygoteRules.pullback(TuringUniform, args...)
end

## Beta ##

function _betalogpdfgrad(α, β, x)
    di = digamma(α + β)
    dα = log(x) - digamma(α) + di
    dβ = log(1 - x) - digamma(β) + di
    dx = (α - 1)/x + (1 - β)/(1 - x)
    return (dα, dβ, dx)
end
ZygoteRules.@adjoint function betalogpdf(α::Real, β::Real, x::Number)
    return betalogpdf(α, β, x), Δ -> (Δ .* _betalogpdfgrad(α, β, x))
end

## Gamma ##

function _gammalogpdfgrad(k, θ, x)
    dk = -digamma(k) - log(θ) + log(x)
    dθ = -k/θ + x/θ^2
    dx = (k - 1)/x - 1/θ
    return (dk, dθ, dx)
end
ZygoteRules.@adjoint function gammalogpdf(k::Real, θ::Real, x::Number)
    return gammalogpdf(k, θ, x), Δ -> (Δ .* _gammalogpdfgrad(k, θ, x))
end

## Chisq ##

function _chisqlogpdfgrad(k, x)
    hk = k/2
    d = digamma(hk)
    dk = (-log(oftype(hk, 2)) - d + log(x))/2
    dx = (hk - 1)/x - one(hk)/2
    return (dk, dx)
end
ZygoteRules.@adjoint function chisqlogpdf(k::Real, x::Number)
    return chisqlogpdf(k, x), Δ -> (Δ .* _chisqlogpdfgrad(k, x))
end

## FDist ##

function _fdistlogpdfgrad(v1, v2, x)
    temp1 = v1 * x + v2
    temp2 = log(temp1)
    vsum = v1 + v2
    temp3 = vsum / temp1
    temp4 = digamma(vsum / 2)
    dv1 = (log(v1 * x) + 1 - temp2 - x * temp3 - digamma(v1 / 2) + temp4) / 2
    dv2 = (log(v2) + 1 - temp2 - temp3 - digamma(v2 / 2) + temp4) / 2
    dx = v1 / 2 * (1 / x - temp3) - 1 / x
    return (dv1, dv2, dx)
end
ZygoteRules.@adjoint function fdistlogpdf(v1::Real, v2::Real, x::Number)
    return fdistlogpdf(v1, v2, x), Δ -> (Δ .* _fdistlogpdfgrad(v1, v2, x))
end

## TDist ##

function _tdistlogpdfgrad(v, x)
    dv = (digamma((v + 1) / 2) - 1 / v - digamma(v / 2) - log(1 + x^2 / v) + x^2 * (v + 1) / v^2 / (1 + x^2 / v)) / 2
    dx = -x * (v + 1) / (v + x^2)
    return (dv, dx)
end
ZygoteRules.@adjoint function tdistlogpdf(v::Real, x::Number)
    return tdistlogpdf(v, x), Δ -> (Δ .* _tdistlogpdfgrad(v, x))
end

## Semicircle ##

function semicircle_dldr(r, x)
    diffsq = r^2 - x^2
    return -2 / r + r / diffsq
end
function semicircle_dldx(r, x)
    diffsq = r^2 - x^2
    return -x / diffsq
end

logpdf(d::Semicircle{<:Real}, x::TrackedReal) = semicirclelogpdf(d.r, x)
logpdf(d::Semicircle{<:TrackedReal}, x::Real) = semicirclelogpdf(d.r, x)
logpdf(d::Semicircle{<:TrackedReal}, x::TrackedReal) = semicirclelogpdf(d.r, x)

semicirclelogpdf(r, x) = logpdf(Semicircle(r), x)
M, f, arity = DiffRules.@define_diffrule DistributionsAD.semicirclelogpdf(r, x) =
    :(semicircle_dldr($r, $x)), :(semicircle_dldx($r, $x))
da, db = DiffRules.diffrule(M, f, :a, :b)
f = :($M.$f)
@eval begin
    @grad $f(a::TrackedReal, b::TrackedReal) = $f(data(a), data(b)), Δ -> (Δ * $da, Δ * $db)
    @grad $f(a::TrackedReal, b::Real) = $f(data(a), b), Δ -> (Δ * $da, Tracker._zero(b))
    @grad $f(a::Real, b::TrackedReal) = $f(a, data(b)), Δ -> (Tracker._zero(a), Δ * $db)
    $f(a::TrackedReal, b::TrackedReal)  = track($f, a, b)
    $f(a::TrackedReal, b::Real) = track($f, a, b)
    $f(a::Real, b::TrackedReal) = track($f, a, b)
end
if VERSION < v"1.2"
    Base.inv(::Irrational{:π}) = 1/π
end

## Binomial ##

binomlogpdf(n::Int, p::TrackedReal, x::Int) = track(binomlogpdf, n, p, x)
@grad function binomlogpdf(n::Int, p::TrackedReal, x::Int)
    return binomlogpdf(n, data(p), x),
        Δ->(nothing, Δ * (x / p - (n - x) / (1 - p)), nothing)
end
ZygoteRules.@adjoint function binomlogpdf(n::Int, p::Real, x::Int)
    return binomlogpdf(n, p, x),
        Δ->(nothing, Δ * (x / p - (n - x) / (1 - p)), nothing)
end

function binomlogpdf(n::Int, p::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(p)
    Δ = ForwardDiff.partials(p)
    return FD(binomlogpdf(n, val, x),  Δ * (x / val - (n - x) / (1 - val)))
end

## Negative binomial ##

# Note the definition of NegativeBinomial in Julia is not the same as Wikipedia's.
# Check the docstring of NegativeBinomial, r is the number of successes and
# k is the number of failures
_nbinomlogpdf_grad_1(r, p, k) = k == 0 ? log(p) : sum(1 / (k + r - i) for i in 1:k) + log(p)
_nbinomlogpdf_grad_2(r, p, k) = -k / (1 - p) + r / p

nbinomlogpdf(n::TrackedReal, p::TrackedReal, x::Int) = track(nbinomlogpdf, n, p, x)
nbinomlogpdf(n::Real, p::TrackedReal, x::Int) = track(nbinomlogpdf, n, p, x)
nbinomlogpdf(n::TrackedReal, p::Real, x::Int) = track(nbinomlogpdf, n, p, x)
@grad function nbinomlogpdf(r::TrackedReal, p::TrackedReal, k::Int)
    return nbinomlogpdf(data(r), data(p), k),
        Δ->(Δ * _nbinomlogpdf_grad_1(r, p, k), Δ * _nbinomlogpdf_grad_2(r, p, k), nothing)
end
@grad function nbinomlogpdf(r::Real, p::TrackedReal, k::Int)
    return nbinomlogpdf(data(r), data(p), k),
        Δ->(Tracker._zero(r), Δ * _nbinomlogpdf_grad_2(r, p, k), nothing)
end
@grad function nbinomlogpdf(r::TrackedReal, p::Real, k::Int)
    return nbinomlogpdf(data(r), data(p), k),
        Δ->(Δ * _nbinomlogpdf_grad_1(r, p, k), Tracker._zero(p), nothing)
end

function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    val_r = ForwardDiff.value(r)
    Δ_r = ForwardDiff.partials(r)
    dr = _nbinomlogpdf_grad_1(val_r, val_p, k)
    Δ_p = ForwardDiff.partials(p)
    dp = _nbinomlogpdf_grad_2(val_r, val_p, k)
    Δ = ForwardDiff._mul_partials(Δ_r, Δ_p, dr, dp)
    return FD(nbinomlogpdf(val_r, val_p, k),  Δ)
end
function nbinomlogpdf(r::Real, p::ForwardDiff.Dual{T}, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_p = ForwardDiff.value(p)
    Δ_p = ForwardDiff.partials(p) * _nbinomlogpdf_grad_2(r, val_p, k)
    return FD(nbinomlogpdf(r, val_p, k),  Δ_p)
end
function nbinomlogpdf(r::ForwardDiff.Dual{T}, p::Real, k::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val_r = ForwardDiff.value(r)
    Δ_r = ForwardDiff.partials(r) * _nbinomlogpdf_grad_1(val_r, p, k)
    return FD(nbinomlogpdf(val_r, p, k),  Δ_r)
end

## Integer dual ##

function BetaBinomial(n::ForwardDiff.Dual{<:Any, <:Integer}, α::Real, β::Real; check_args = true)
    return BetaBinomial(ForwardDiff.value(n), α, β; check_args = check_args)
end
Binomial(n::ForwardDiff.Dual{<:Any, <:Integer}, p::Real) = Binomial(ForwardDiff.value(n), p)
function Erlang(α::ForwardDiff.Dual{<:Any, <:Integer}, θ::Real; check_args = true)
    return Erlang(ForwardDiff.value(α), θ, check_args = check_args)
end

## Poisson ##

poislogpdf(v::TrackedReal, x::Int) = track(poislogpdf, v, x)
@grad function poislogpdf(v::TrackedReal, x::Int)
      return poislogpdf(data(v), x),
          Δ->(Δ * (x/v - 1), nothing)
end
ZygoteRules.@adjoint function poislogpdf(v::Real, x::Int)
    return poislogpdf(v, x),
        Δ->(Δ * (x/v - 1), nothing)
end

function poislogpdf(v::ForwardDiff.Dual{T}, x::Int) where {T}
    FD = ForwardDiff.Dual{T}
    val = ForwardDiff.value(v)
    Δ = ForwardDiff.partials(v)
    return FD(poislogpdf(val, x), Δ * (x/val - 1))
end

## PoissonBinomial ##

struct TuringPoissonBinomial{T<:Real, TV1<:AbstractVector{T}, TV2<:AbstractVector} <: DiscreteUnivariateDistribution
    p::TV1
    pmf::TV2
end
function TuringPoissonBinomial(p::AbstractArray{<:Real}; check_args = true)
    pb = Distributions.poissonbinomial_pdf_fft(p)
    ϵ = eps(eltype(pb))
    check_args && @assert all(x -> x >= -ϵ, pb) && isapprox(sum(pb), 1, atol = ϵ)
    return TuringPoissonBinomial(p, pb)
end
function logpdf(d::TuringPoissonBinomial{T}, k::Int) where T<:Real
    insupport(d, k) ? log(d.pmf[k + 1]) : -T(Inf)
end
quantile(d::TuringPoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1
PoissonBinomial(p::TrackedArray{<:Real}; check_args=true) = TuringPoissonBinomial(p; check_args = check_args)
Base.minimum(d::TuringPoissonBinomial) = 0
Base.maximum(d::TuringPoissonBinomial) = length(d.p)

poissonbinomial_pdf_fft(x::TrackedArray) = track(poissonbinomial_pdf_fft, x)
@grad function poissonbinomial_pdf_fft(x::TrackedArray)
    x_data = data(x)
    T = eltype(x_data)
    fft = poissonbinomial_pdf_fft(x_data)
    return  fft, Δ -> begin
        ((ForwardDiff.jacobian(x -> poissonbinomial_pdf_fft(x), x_data)::Matrix{T})' * Δ,)
    end
end
# FIXME: This is inefficient, replace with the commented code below once Zygote supports it.
ZygoteRules.@adjoint function poissonbinomial_pdf_fft(x::AbstractArray)
    T = eltype(x)
    fft = poissonbinomial_pdf_fft(x)
    return  fft, Δ -> begin
        ((ForwardDiff.jacobian(x -> poissonbinomial_pdf_fft(x), x)::Matrix{T})' * Δ,)
    end
end

# The code below doesn't work because of bugs in Zygote. The above is inefficient.
#=
ZygoteRules.@adjoint function poissonbinomial_pdf_fft(x::AbstractArray{<:Real})
    return ZygoteRules.pullback(poissonbinomial_pdf_fft_zygote, x)
end
function poissonbinomial_pdf_fft_zygote(p::AbstractArray{T}) where {T <: Real}
    n = length(p)
    ω = 2 * one(T) / (n + 1)

    lmax = ceil(Int, n/2)
    x1 = [one(T)/(n + 1)]
    x_lmaxp1 = map(1:lmax) do l
        logz = zero(T)
        argz = zero(T)
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        return dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
    end
    x_lmaxp2_end = [conj(x[l + 1]) for l in lmax:-1:1 if n + 1 - l > l]
    x = vcat(x1; x_lmaxp1, x_lmaxp2_end)
    y = [sum(x[j] * cis(-π * float(T)(2 * mod(j * k, n)) / n) for j in 1:n) for k in 1:n]
    return max.(0, real.(y))
end
function poissonbinomial_pdf_fft_zygote2(p::AbstractArray{T}) where {T <: Real}
    n = length(p)
    ω = 2 * one(T) / (n + 1)

    x = Vector{Complex{T}}(undef, n+1)
    lmax = ceil(Int, n/2)
    x[1] = one(T)/(n + 1)
    for l=1:lmax
        logz = zero(T)
        argz = zero(T)
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        x[l + 1] = dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
        if n + 1 - l > l
            x[n + 1 - l + 1] = conj(x[l + 1])
        end
    end
    max.(0, real.(_dft_zygote(copy(x))))
end
function _dft_zygote(x::Vector{T}) where T
    n = length(x)
    y = Zygote.Buffer(zeros(complex(float(T)), n))
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(-π * float(T)(2 * mod(j * k, n)) / n)
    end
    return copy(y)
end
=#

## Categorical ##

function Base.convert(
    ::Type{Distributions.DiscreteNonParametric{T,P,Ts,Ps}},
    d::Distributions.DiscreteNonParametric{T,P,Ts,Ps},
) where {T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:AbstractVector{P}}
    DiscreteNonParametric{T,P,Ts,Ps}(support(d), probs(d), check_args=false)
end

# Fix SubArray support
function Distributions.DiscreteNonParametric{T,P,Ts,Ps}(
    vs::Ts,
    ps::Ps;
    check_args=true,
) where {T<:Real, P<:Real, Ts<:AbstractVector{T}, Ps<:SubArray{P, 1}}
    cps = ps[:]
    return DiscreteNonParametric{T,P,Ts,typeof(cps)}(vs, cps; check_args = check_args)
end

function Distributions.DiscreteNonParametric{T,P,Ts,Ps}(
    vs::Ts,
    ps::Ps;
    check_args=true,
) where {T<:Real,P<:Real,Ts<:AbstractVector{T},Ps<:TrackedArray{P, 1, <:SubArray{P, 1}}}
    cps = ps[:]
    return DiscreteNonParametric{T,P,Ts,typeof(cps)}(vs, cps; check_args = check_args)
end

## ForwardDiff broadcasting support ##

function Distributions.logpdf(d::DiscreteUnivariateDistribution, k::ForwardDiff.Dual)
    return logpdf(d, convert(Integer, ForwardDiff.value(k)))
end


text:
## Landcover variables loading functions

"""
    Copernicus

Data from the Copernicus v2.0.1 land cover data from the Copernicus Glocal Land
Service. The original data by Buchhorn et al. is archived on Zenodo at
https://zenodo.org/record/3243509. See `SimpleSDMPredictor(Copernicus,
LandCover, ...)` for the list of variables.
"""
struct Copernicus end

import SimpleSDMLayers: SimpleSDMPredictor

"""
    SimpleSDMPredictor(::Type{Copernicus}, ::Type{LandCover}, layer::Integer=1; resolution::Float64=10.0, path::AbstractString=joinpath("assets", "landcover"), kwargs...)
    SimpleSDMPredictor(::Type{Copernicus}, ::Type{LandCover}, layers::AbstractArray; kwargs...)

Prepares the Copernicus v2.0.1 landcover data, and returns them as an array of
`SimpleSDMPredictor`s. Layers are called by their number, from 1 to 10. The list
of available layers is given in the table below.

This function loads the pre-downloaded and coarsened data from the
`assets/landcover/` folder (the `path` argument). These steps can be performed
with the scripts in `src/shell/` using Bash and GDAL. The available resolutions
are 10 and 5 arc-minutes, and must be specified as a floating point value with
the `resolution` keyword (with `10.0` as the default).

The original data by Buchhorn et al. is archived on Zenodo at
https://zenodo.org/record/3243509, and is available at a much finer resolution.

| Variable | Description        |
| ------   | ------             |
| 1        | Bare               |
| 2        | Crops              |
| 3        | Grass              |
| 4        | Moss               |
| 5        | Shrub              |
| 6        | Snow               |
| 7        | Tree               |
| 8        | Urban              |
| 9        | Water permanent    |
| 10       | Water seasonal     |

"""
function SimpleSDMPredictor(
    ::Type{Copernicus},
    ::Type{LandCover},
    layer::Integer=1;
    resolution::Float64=10.0,
    path::AbstractString=joinpath("assets", "landcover"),
    kwargs...,
)
    @assert resolution in [5.0, 10.0]
    1 ≤ layer ≤ 10 || throw(ArgumentError("The layer must be between 1 and 10"))

    # List files in path
    lc_files = readdir(path)
    # Filter for selected resolution
    filter!(contains("$(Int(resolution))m.tif"), lc_files)
    # Create path for selected layer only
    p = joinpath.(path, lc_files)[layer]

    # Read layer
    layer = geotiff(SimpleSDMPredictor, p; kwargs...)
    # Convert from UInt8 to Float32
    layer = convert(Float32, layer)

    return convert(SimpleSDMPredictor, layer)
end

function SimpleSDMPredictor(
    ::Type{Copernicus}, ::Type{LandCover}, layers::AbstractArray; kwargs...
)
    @assert eltype(layers) <: Integer
    return [SimpleSDMPredictor(Copernicus, LandCover, l; kwargs...) for l in layers]
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using Routing
using Test

include("debugging.jl")

test_start_time = time()

include("espprc_test.jl")
include("vrptw_test.jl")

@testset "Solomon Instance Test" begin
    solomon_dataset = Dict(
        "C101_025" => 191.30,
        "C102_025" => 190.30,
        "C201_025" => 214.70,
        "RC105_025" => 411.30,
        "R101_025" => 617.10,
        "R102_025" => 547.10,
    )

    for (name, val) in solomon_dataset
        @info("Solving the Solomon Instance $name")
        solomon = load_solomon(name)
        vrptw = generate_solomon_vrptw_instance(solomon)

        @testset "$name" begin
            println("-- Solomon Instance $name --")
            @time routes, objective_value = solveVRP(vrptw, pricing_method="pulse");
            @test isapprox(objective_value, val, atol=1e-7)
        end
    end
end

total_test_time = time() - test_start_time
@show total_test_time


text:
struct TapedTaskException
    exc::Exception
    backtrace::Vector{Any}
end

struct TapedTask{F}
    task::Task
    tf::TapedFunction{F}
    produce_ch::Channel{Any}
    consume_ch::Channel{Int}
    produced_val::Vector{Any}

    function TapedTask(
        t::Task,
        tf::TapedFunction{F},
        produce_ch::Channel{Any},
        consume_ch::Channel{Int}
    ) where {F}
        new{F}(t, tf, produce_ch, consume_ch, Any[])
    end
end

function producer()
    ttask = current_task().storage[:tapedtask]::TapedTask
    if length(ttask.produced_val) > 0
        val = pop!(ttask.produced_val)
        put!(ttask.produce_ch, val)
        take!(ttask.consume_ch) # wait for next consumer
    end
    return nothing
end

function wrap_task(tf, produce_ch, consume_ch, args...)
    try
        tf(args...; callback=producer, continuation=true)
    catch e
        bt = catch_backtrace()
        put!(produce_ch, TapedTaskException(e, bt))
        # @error "TapedTask Error: " exception=(e, bt)
        rethrow()
    finally
        @static if VERSION >= v"1.4"
            # we don't do this under Julia 1.3, because `isempty` always hangs on
            # an empty channel.
            while !isempty(produce_ch)
                yield()
            end
        end
        close(produce_ch)
        close(consume_ch)
    end
end

function TapedTask(tf::TapedFunction, args...)
    produce_ch = Channel()
    consume_ch = Channel{Int}()
    task = @task wrap_task(tf, produce_ch, consume_ch, args...)
    t = TapedTask(task, tf, produce_ch, consume_ch)
    task.storage === nothing && (task.storage = IdDict())
    task.storage[:tapedtask] = t
    return t
end

# NOTE: evaluating model without a trace, see
# https://github.com/TuringLang/Turing.jl/pull/1757#diff-8d16dd13c316055e55f300cd24294bb2f73f46cbcb5a481f8936ff56939da7ceR329
function TapedTask(f, args...)
    tf = TapedFunction(f, args...; cache=true)
    TapedTask(tf, args...)
end

TapedTask(t::TapedTask, args...) = TapedTask(func(t), args...)
func(t::TapedTask) = t.tf.func

#=
# ** Approach (A) to implement `produce`:
# Make`produce` a standalone instturction. This approach does NOT
# support `produce` in a nested call
function internal_produce(instr::Instruction, val)
    tf = instr.tape
    ttask = tf.owner
    put!(ttask.produce_ch, val)
    take!(ttask.consume_ch) # wait for next consumer
end

function produce(val)
    error("Libtask.produce can only be directly called in a task!")
end

function (instr::Instruction{typeof(produce)})()
    args = val(instr.input[1])
    internal_produce(instr, args)
end
=#

# ** Approach (B) to implement `produce`:
# This way has its caveat:
# `produce` may deeply hide in an instruction, but not be an instruction
# itself, and when we copy a task, the newly copied task will resume from
# the instruction after the one which contains this `produce` call. If the
# call to `produce` is not the last expression in the instuction, that
# instruction will not be whole executed in the copied task.
@inline function is_in_tapedtask()
    ct = current_task()
    ct.storage === nothing && return false
    haskey(ct.storage, :tapedtask) || return false
    # check if we are recording a tape
    ttask = ct.storage[:tapedtask]::TapedTask
    return !isempty(ttask.tf.tape)
end

function produce(val)
    is_in_tapedtask() || return nothing
    ttask = current_task().storage[:tapedtask]::TapedTask
    length(ttask.produced_val) > 1 &&
        error("There is a produced value which is not consumed.")
    push!(ttask.produced_val, val)
    return nothing
end

function consume(ttask::TapedTask)
    if istaskstarted(ttask.task)
        # tell producer that a consumer is coming
        put!(ttask.consume_ch, 0)
    else
        schedule(ttask.task)
    end

    val = try
        take!(ttask.produce_ch)
    catch e
        isa(e, InvalidStateException) || rethrow()
        istaskfailed(ttask.task) && throw(ttask.task.exception)
        # TODO: we return nothing to indicate the end of a task,
        #       remove this when AdvancedPS is udpated.
        istaskdone(ttask.task) && return nothing
    end

    # yield to let the task resume, this is necessary when there's
    # an exception is thrown in the task, it gives the task the chance
    # to rethow the exception and set its proper status:
    yield()
    isa(val, TapedTaskException) && throw(val.exc)
    return val
end

# Iteration interface.
function Base.iterate(t::TapedTask, state=nothing)
    try
        consume(t), nothing
    catch ex
        !isa(ex, InvalidStateException) && rethrow
        nothing
    end
end
Base.IteratorSize(::Type{<:TapedTask}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:TapedTask}) = Base.EltypeUnknown()

# copy the task

function Base.copy(t::TapedTask)
    tf = copy(t.tf)
    new_t = TapedTask(tf)
    storage = t.task.storage::IdDict{Any,Any}
    new_t.task.storage = copy(storage)
    new_t.task.storage[:tapedtask] = new_t
    return new_t
end

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
"""
    Materials

Module will weave all the Julia markdown files to HTML.
"""
module Materials
    using Weave: weave

    function weave_all()
        for (root, dirs, files) ∈ walkdir(@__DIR__)
            for file ∈ filter!(file -> endswith(file, ".jmd"), files)
                weave(joinpath(root, file)) # path to files
            end
        end
    end
    export weave_all
end

max_stars_repo_name: Nosferican/Introduction-to-Julia

max_issues_repo_name: Nosferican/Introduction-to-Julia

max_forks_repo_name: Nosferican/Introduction-to-Julia


text:
module TestCat
    using Base.Test
    using DataTables

    #
    # hcat
    #

    nvint = NullableArray(Nullable{Int}[1, 2, Nullable(), 4])
    nvstr = NullableArray(Nullable{String}["one", "two", Nullable(), "four"])

    dt2 = DataTable(Any[nvint, nvstr])
    dt3 = DataTable(Any[nvint])
    dt4 = convert(DataTable, [1:4 1:4])
    dt5 = DataTable(Any[NullableArray([1,2,3,4]), nvstr])

    dth = hcat(dt3, dt4)
    @test size(dth, 2) == 3
    @test names(dth) == [:x1, :x1_1, :x2]
    @test isequal(dth[:x1], dt3[:x1])
    @test isequal(dth, [dt3 dt4])
    @test isequal(dth, DataTables.hcat!(DataTable(), dt3, dt4))

    dth3 = hcat(dt3, dt4, dt5)
    @test names(dth3) == [:x1, :x1_1, :x2, :x1_2, :x2_1]
    @test isequal(dth3, hcat(dth, dt5))
    @test isequal(dth3, DataTables.hcat!(DataTable(), dt3, dt4, dt5))

    @test isequal(dt2, DataTables.hcat!(dt2))

    @testset "hcat ::AbstractDataTable" begin
        dt = DataTable(A = repeat('A':'C', inner=4), B = 1:12)
        gd = groupby(dt, :A)
        answer = DataTable(A = fill('A', 4), B = 1:4, A_1 = 'B', B_1 = 5:8, A_2 = 'C', B_2 = 9:12)
        @test hcat(gd...) == answer
        answer = answer[1:4]
        @test hcat(gd[1], gd[2]) == answer
    end

    @testset "hcat ::Vectors" begin
        dt = DataTable()
        DataTables.hcat!(dt, NullableCategoricalVector(1:10))
        @test isequal(dt[1], NullableCategoricalVector(1:10))
        DataTables.hcat!(dt, NullableArray(1:10))
        @test isequal(dt[2], NullableArray(1:10))
    end

    #
    # vcat
    #

    null_dt = DataTable(Int, 0, 0)
    dt = DataTable(Int, 4, 3)

    # Assignment of rows
    dt[1, :] = dt[1, :]
    dt[1:2, :] = dt[1:2, :]
    dt[[true,false,false,true], :] = dt[2:3, :]

    # Scalar broadcasting assignment of rows
    dt[1, :] = 1
    dt[1:2, :] = 1
    dt[[true,false,false,true], :] = 3

    # Vector broadcasting assignment of rows
    dt[1:2, :] = [2,3]
    dt[[true,false,false,true], :] = [2,3]

    # Assignment of columns
    dt[1] = zeros(4)
    dt[:, 2] = ones(4)

    # Broadcasting assignment of columns
    dt[:, 1] = 1
    dt[1] = 3
    dt[:x3] = 2

    # assignment of subtables
    dt[1, 1:2] = dt[2, 2:3]
    dt[1:2, 1:2] = dt[2:3, 2:3]
    dt[[true,false,false,true], 2:3] = dt[1:2,1:2]

    # scalar broadcasting assignment of subtables
    dt[1, 1:2] = 3
    dt[1:2, 1:2] = 3
    dt[[true,false,false,true], 2:3] = 3

    # vector broadcasting assignment of subtables
    dt[1:2, 1:2] = [3,2]
    dt[[true,false,false,true], 2:3] = [2,3]

    @test vcat(null_dt) == DataTable()
    @test vcat(null_dt, null_dt) == DataTable()
    @test_throws ArgumentError vcat(null_dt, dt)
    @test_throws ArgumentError vcat(dt, null_dt)
    @test eltypes(vcat(dt, dt)) == Type[Float64, Float64, Int]
    @test size(vcat(dt, dt)) == (size(dt,1)*2, size(dt,2))
    @test eltypes(vcat(dt, dt, dt)) == Type[Float64,Float64,Int]
    @test size(vcat(dt, dt, dt)) == (size(dt,1)*3, size(dt,2))

    alt_dt = deepcopy(dt)
    vcat(dt, alt_dt)

    # Don't fail on non-matching types
    dt[1] = zeros(Int, nrow(dt))
    vcat(dt, alt_dt)

    dtr = vcat(dt4, dt4)
    @test size(dtr, 1) == 8
    @test names(dt4) == names(dtr)
    @test isequal(dtr, [dt4; dt4])

    @test eltypes(vcat(DataTable(a = [1]), DataTable(a = [2.1]))) == Type[Float64]
    @test eltypes(vcat(DataTable(a = NullableArray(Int, 1)), DataTable(a = [2.1]))) == Type[Nullable{Float64}]

    # Minimal container type promotion
    dta = DataTable(a = NullableCategoricalArray([1, 2, 2]))
    dtb = DataTable(a = NullableCategoricalArray([2, 3, 4]))
    dtc = DataTable(a = NullableArray([2, 3, 4]))
    dtd = DataTable(Any[2:4], [:a])
    dtab = vcat(dta, dtb)
    dtac = vcat(dta, dtc)
    @test isequal(dtab[:a], Nullable{Int}[1, 2, 2, 2, 3, 4])
    @test isequal(dtac[:a], Nullable{Int}[1, 2, 2, 2, 3, 4])
    @test isa(dtab[:a], NullableCategoricalVector{Int})
    @test isa(dtac[:a], NullableCategoricalVector{Int})
    # ^^ container may flip if container promotion happens in Base/DataArrays
    dc = vcat(dtd, dtc)
    @test isequal(vcat(dtc, dtd), dc)

    # Zero-row DataTables
    dtc0 = similar(dtc, 0)
    @test isequal(vcat(dtd, dtc0, dtc), dc)
    @test eltypes(vcat(dtd, dtc0)) == eltypes(dc)

    # vcat should be able to concatenate different implementations of AbstractDataTable (PR #944)
    @test isequal(vcat(view(DataTable(A=1:3),2),DataTable(A=4:5)), DataTable(A=[2,4,5]))

    @testset "vcat >2 args" begin
        @test vcat(DataTable(), DataTable(), DataTable()) == DataTable()
        dt = DataTable(x = trues(1), y = falses(1))
        @test vcat(dt, dt, dt) == DataTable(x = trues(3), y = falses(3))
    end

    @testset "vcat mixed coltypes" begin
        drf = CategoricalArrays.DefaultRefType
        dt = vcat(DataTable([[1]], [:x]), DataTable([[1.0]], [:x]))
        @test dt == DataTable([[1.0, 1.0]], [:x])
        @test typeof.(dt.columns) == [Vector{Float64}]
        dt = vcat(DataTable([[1]], [:x]), DataTable([["1"]], [:x]))
        @test dt == DataTable([[1, "1"]], [:x])
        @test typeof.(dt.columns) == [Vector{Any}]
        dt = vcat(DataTable([NullableArray([1])], [:x]), DataTable([[1]], [:x]))
        @test dt == DataTable([NullableArray([1, 1])], [:x])
        @test typeof.(dt.columns) == [NullableVector{Int}]
        dt = vcat(DataTable([CategoricalArray([1])], [:x]), DataTable([[1]], [:x]))
        @test dt == DataTable([CategoricalArray([1, 1])], [:x])
        @test typeof.(dt.columns) == [CategoricalVector{Int, drf}]
        dt = vcat(DataTable([CategoricalArray([1])], [:x]),
                  DataTable([NullableArray([1])], [:x]))
        @test dt == DataTable([NullableCategoricalArray([1, 1])], [:x])
        @test typeof.(dt.columns) == [NullableCategoricalVector{Int, drf}]
        dt = vcat(DataTable([CategoricalArray([1])], [:x]),
                  DataTable([NullableCategoricalArray([1])], [:x]))
        @test dt == DataTable([NullableCategoricalArray([1, 1])], [:x])
        @test typeof.(dt.columns) == [NullableCategoricalVector{Int, drf}]
        dt = vcat(DataTable([NullableArray([1])], [:x]),
                  DataTable([NullableArray(["1"])], [:x]))
        @test dt == DataTable([NullableArray([1, "1"])], [:x])
        @test typeof.(dt.columns) == [NullableVector{Any}]
        dt = vcat(DataTable([CategoricalArray([1])], [:x]),
                  DataTable([CategoricalArray(["1"])], [:x]))
        @test dt == DataTable([CategoricalArray([1, "1"])], [:x])
        @test typeof.(dt.columns) == [CategoricalVector{Any, drf}]
        dt = vcat(DataTable([trues(1)], [:x]), DataTable([[false]], [:x]))
        @test dt == DataTable([[true, false]], [:x])
        @test typeof.(dt.columns) == [Vector{Bool}]
    end

    @testset "vcat errors" begin
        err = @test_throws ArgumentError vcat(DataTable(), DataTable(), DataTable(x=[]))
        @test err.value.msg == "column(s) x are missing from argument(s) 1 and 2"
        err = @test_throws ArgumentError vcat(DataTable(), DataTable(), DataTable(x=[1]))
        @test err.value.msg == "column(s) x are missing from argument(s) 1 and 2"
        dt1 = DataTable(A = 1:3, B = 1:3)
        dt2 = DataTable(A = 1:3)
        # right missing 1 column
        err = @test_throws ArgumentError vcat(dt1, dt2)
        @test err.value.msg == "column(s) B are missing from argument(s) 2"
        # left missing 1 column
        err = @test_throws ArgumentError vcat(dt2, dt1)
        @test err.value.msg == "column(s) B are missing from argument(s) 1"
        # multiple missing 1 column
        err = @test_throws ArgumentError vcat(dt1, dt2, dt2, dt2, dt2, dt2)
        @test err.value.msg == "column(s) B are missing from argument(s) 2, 3, 4, 5 and 6"
        # argument missing >1 columns
        dt1 = DataTable(A = 1:3, B = 1:3, C = 1:3, D = 1:3, E = 1:3)
        err = @test_throws ArgumentError vcat(dt1, dt2)
        @test err.value.msg == "column(s) B, C, D and E are missing from argument(s) 2"
        # >1 arguments missing >1 columns
        err = @test_throws ArgumentError vcat(dt1, dt2, dt2, dt2, dt2)
        @test err.value.msg == "column(s) B, C, D and E are missing from argument(s) 2, 3, 4 and 5"
        # out of order
        dt2 = dt1[reverse(names(dt1))]
        err = @test_throws ArgumentError vcat(dt1, dt2)
        @test err.value.msg == "column order of argument(s) 1 != column order of argument(s) 2"
        # first group >1 arguments
        err = @test_throws ArgumentError vcat(dt1, dt1, dt2)
        @test err.value.msg == "column order of argument(s) 1 and 2 != column order of argument(s) 3"
        # second group >1 arguments
        err = @test_throws ArgumentError vcat(dt1, dt2, dt2)
        @test err.value.msg == "column order of argument(s) 1 != column order of argument(s) 2 and 3"
        # first and second groups >1 argument
        err = @test_throws ArgumentError vcat(dt1, dt1, dt1, dt2, dt2, dt2)
        @test err.value.msg == "column order of argument(s) 1, 2 and 3 != column order of argument(s) 4, 5 and 6"
        # >2 groups out of order
        srand(1)
        dt3 = dt1[shuffle(names(dt1))]
        err = @test_throws ArgumentError vcat(dt1, dt1, dt1, dt2, dt2, dt2, dt3, dt3, dt3, dt3)
        @test err.value.msg == "column order of argument(s) 1, 2 and 3 != column order of argument(s) 4, 5 and 6 != column order of argument(s) 7, 8, 9 and 10"
        # missing columns throws error before out of order columns
        dt1 = DataTable(A = 1, B = 1)
        dt2 = DataTable(A = 1)
        dt3 = DataTable(B = 1, A = 1)
        err = @test_throws ArgumentError vcat(dt1, dt2, dt3)
        @test err.value.msg == "column(s) B are missing from argument(s) 2"
        # unique columns for both sides
        dt1 = DataTable(A = 1, B = 1, C = 1, D = 1)
        dt2 = DataTable(A = 1, C = 1, D = 1, E = 1, F = 1)
        err = @test_throws ArgumentError vcat(dt1, dt2)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, and column(s) B are missing from argument(s) 2"
        err = @test_throws ArgumentError vcat(dt1, dt1, dt2, dt2)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1 and 2, and column(s) B are missing from argument(s) 3 and 4"
        dt3 = DataTable(A = 1, B = 1, C = 1, D = 1, E = 1)
        err = @test_throws ArgumentError vcat(dt1, dt2, dt3)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, column(s) B are missing from argument(s) 2, and column(s) F are missing from argument(s) 3"
        err = @test_throws ArgumentError vcat(dt1, dt1, dt2, dt2, dt3, dt3)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1 and 2, column(s) B are missing from argument(s) 3 and 4, and column(s) F are missing from argument(s) 5 and 6"
        err = @test_throws ArgumentError vcat(dt1, dt1, dt1, dt2, dt2, dt2, dt3, dt3, dt3)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, 2 and 3, column(s) B are missing from argument(s) 4, 5 and 6, and column(s) F are missing from argument(s) 7, 8 and 9"
        # dt4 is a superset of names found in all other datatables and won't be shown in error
        dt4 = DataTable(A = 1, B = 1, C = 1, D = 1, E = 1, F = 1)
        err = @test_throws ArgumentError vcat(dt1, dt2, dt3, dt4)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, column(s) B are missing from argument(s) 2, and column(s) F are missing from argument(s) 3"
        err = @test_throws ArgumentError vcat(dt1, dt1, dt2, dt2, dt3, dt3, dt4, dt4)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1 and 2, column(s) B are missing from argument(s) 3 and 4, and column(s) F are missing from argument(s) 5 and 6"
        err = @test_throws ArgumentError vcat(dt1, dt1, dt1, dt2, dt2, dt2, dt3, dt3, dt3, dt4, dt4, dt4)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, 2 and 3, column(s) B are missing from argument(s) 4, 5 and 6, and column(s) F are missing from argument(s) 7, 8 and 9"
        err = @test_throws ArgumentError vcat(dt1, dt2, dt3, dt4, dt1, dt2, dt3, dt4, dt1, dt2, dt3, dt4)
        @test err.value.msg == "column(s) E and F are missing from argument(s) 1, 5 and 9, column(s) B are missing from argument(s) 2, 6 and 10, and column(s) F are missing from argument(s) 3, 7 and 11"
    end
    x = view(DataTable(A = NullableArray(1:3)), 2)
    y = DataTable(A = NullableArray(4:5))
    @test isequal(vcat(x, y), DataTable(A = NullableArray([2, 4, 5])))
end


text:
function fib20()
#
#  test polynomial suggested by Goedecker
#
    n = 20;
    p = fib(n)
    z = [-.9473885474109854;
         -.9029230998844866-.2884421994216345*im; -.9029230998844866+.2884421994216345*im;
         -.7734751441751908-.5504669225440185*im; -.7734751441751908+.5504669225440185*im;
         -.5705411148090332-.7619498890625472*im; -.5705411148090332+.7619498890625472*im;
         -.3121525651549263-.9031354642584495*im; -.3121525651549263+.9031354642584495*im;
         -.02129876087300105-.9602891496489760*im; -.02129876087300105+.9602891496489760*im;
         .2760438985065760-.9267417870492216*im; .2760438985065760+.9267417870492216*im;
         .5529956567245064-.8031600819414702*im; .5529956567245064+.8031600819414702*im;
         .7834063973613306-.5969423289081712*im; .7834063973613306+.5969423289081712*im;
         .9416394828514235-.3213776783384489*im; .9416394828514235+.3213776783384489*im;
         1.999999046316589];
    z = [z ones(n)]
    p, PolyZeros(z)
end


text:
# ---
# title: 891. Sum of Subsequence Widths
# id: problem891
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Math
# link: <https://leetcode.com/problems/sum-of-subsequence-widths/description/>
# hidden: true
# ---
#
# Given an array of integers `A`, consider all non-empty subsequences of `A`.
#
# For any sequence S, let the  _width_  of S be the difference between the
# maximum and minimum element of S.
#
# Return the sum of the widths of all subsequences of A.
#
# As the answer may be very large, **return the answer modulo 10^9 + 7**.
#
#
#
# **Example 1:**
#
#
#
#     Input: [2,1,3]
#     Output: 6
#     Explanation: Subsequences are [1], [2], [3], [2,1], [2,3], [1,3], [2,1,3].
#     The corresponding widths are 0, 0, 0, 1, 1, 2, 2.
#     The sum of these widths is 6.
#
#
#
#
# **Note:**
#
#   * `1 <= A.length <= 20000`
#   * `1 <= A[i] <= 20000`
#
#
## @lc code=start
using LeetCode

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


text:
module Schedulers

include("../docstrings.jl")

using Metatheory.Rules
using Metatheory.EGraphs
using Metatheory.Patterns
using DocStringExtensions

export AbstractScheduler
export SimpleScheduler
export BackoffScheduler
export ScoredScheduler
export cansaturate
export cansearch
export inform!
export setiter!

"""
Represents a rule scheduler for the equality saturation process

"""
abstract type AbstractScheduler end

"""
Should return `true` if the e-graph can be said to be saturated
```
cansaturate(s::AbstractScheduler)
```
"""
function cansaturate end

"""
Should return `false` if the rule `r` should be skipped
```
cansearch(s::AbstractScheduler, r::Rule)
```
"""
function cansearch end

"""
This function is called **after** pattern matching on the e-graph,
informs the scheduler about the yielded matches.
Returns `false` if the matches should not be yielded and ignored.
```
inform!(s::AbstractScheduler, r::AbstractRule, n_matches)
```
"""
function inform! end

function setiter! end

# ===========================================================================
# SimpleScheduler
# ===========================================================================

"""
A simple Rewrite Scheduler that applies every rule every time
"""
struct SimpleScheduler <: AbstractScheduler end

cansaturate(s::SimpleScheduler) = true
cansearch(s::SimpleScheduler, r::AbstractRule) = true
function SimpleScheduler(G::EGraph, theory::Vector{<:AbstractRule})
    SimpleScheduler()
end
inform!(s::SimpleScheduler, r, n_matches) = true
setiter!(s::SimpleScheduler, iteration) = nothing

# ===========================================================================
# BackoffScheduler
# ===========================================================================

mutable struct BackoffSchedulerEntry
    match_limit::Int
    ban_length::Int
    times_banned::Int
    banned_until::Int
end

"""
A Rewrite Scheduler that implements exponential rule backoff.
For each rewrite, there exists a configurable initial match limit.
If a rewrite search yield more than this limit, then we ban this rule
for number of iterations, double its limit, and double the time it
will be banned next time.

This seems effective at preventing explosive rules like
associativity from taking an unfair amount of resources.
"""
mutable struct BackoffScheduler <: AbstractScheduler
    data::IdDict{AbstractRule, BackoffSchedulerEntry}
    G::EGraph
    theory::Vector{<:AbstractRule}
    curr_iter::Int
end

cansearch(s::BackoffScheduler, r::AbstractRule)::Bool = s.curr_iter > s.data[r].banned_until

function BackoffScheduler(g::EGraph, theory::Vector{<:AbstractRule})
    # BackoffScheduler(g, theory, 128, 4)
    BackoffScheduler(g, theory, 1000, 5)
end

function BackoffScheduler(G::EGraph, theory::Vector{<:AbstractRule}, match_limit::Int, ban_length::Int)
    gsize = length(G.uf)
    data = IdDict{AbstractRule, BackoffSchedulerEntry}()

    for rule ∈ theory
        data[rule] = BackoffSchedulerEntry(match_limit, ban_length, 0, 0)
    end

    return BackoffScheduler(data, G, theory, 1)
end

# can saturate if there's no banned rule
cansaturate(s::BackoffScheduler)::Bool = all(kv -> s.curr_iter > last(kv).banned_until, s.data)

function inform!(s::BackoffScheduler, rule::AbstractRule, n_matches)
    # println(s.data[rule])

    rd = s.data[rule]
    treshold = rd.match_limit << rd.times_banned
    if n_matches > treshold
        ban_length = rd.ban_length << rd.times_banned
        rd.times_banned += 1
        rd.banned_until = s.curr_iter + ban_length
        # @info "banning rule $rule until $(rd.banned_until)!"
        return false
    end
    return true
end

function setiter!(s::BackoffScheduler, curr_iter)
    s.curr_iter = curr_iter
end

# ===========================================================================
# ScoredScheduler
# ===========================================================================

mutable struct ScoredSchedulerEntry
    match_limit::Int
    ban_length::Int
    times_banned::Int
    banned_until::Int
    weight::Int
end

"""
A Rewrite Scheduler that implements exponential rule backoff.
For each rewrite, there exists a configurable initial match limit.
If a rewrite search yield more than this limit, then we ban this rule
for number of iterations, double its limit, and double the time it
will be banned next time.

This seems effective at preventing explosive rules like
associativity from taking an unfair amount of resources.
"""
mutable struct ScoredScheduler <: AbstractScheduler
    data::IdDict{AbstractRule, ScoredSchedulerEntry}
    G::EGraph
    theory::Vector{<:AbstractRule}
    curr_iter::Int
end

cansearch(s::ScoredScheduler, r::AbstractRule)::Bool = s.curr_iter > s.data[r].banned_until

exprsize(a) = 1

function exprsize(e::PatTerm)
    c = 1 + length(e.args)
    for a ∈ e.args
        c += exprsize(a)
    end
    return c
end

function exprsize(e::Expr)
    start = Meta.isexpr(e, :call) ? 2 : 1

    c = 1 + length(e.args[start:end])
    for a ∈ e.args[start:end]
        c += exprsize(a)
    end

    return c
end

function ScoredScheduler(g::EGraph, theory::Vector{<:AbstractRule})
    # BackoffScheduler(g, theory, 128, 4)
    ScoredScheduler(g, theory, 1000, 5, exprsize)
end

function ScoredScheduler(G::EGraph, theory::Vector{<:AbstractRule}, match_limit::Int, ban_length::Int, complexity::Function)
    gsize = length(G.uf)
    data = IdDict{AbstractRule, ScoredSchedulerEntry}()

    for rule ∈ theory
        if rule isa DynamicRule
            w = 2
            data[rule] = ScoredSchedulerEntry(match_limit, ban_length, 0, 0, w)
            continue
        end
        (l, r) = rule.left, rule.right

        cl = complexity(l)
        cr = complexity(r)
        # println("$rule HAS SCORE $((cl, cr))")
        if cl > cr
            w = 1   # reduces complexity
        elseif cr > cl
            w = 3   # augments complexity
        else
            w = 2   # complexity is equal
        end
        # println(w)
        data[rule] = ScoredSchedulerEntry(match_limit, ban_length, 0, 0, w)
    end

    return ScoredScheduler(data, G, theory, 1)
end

# can saturate if there's no banned rule
cansaturate(s::ScoredScheduler)::Bool = all(kv -> s.curr_iter > last(kv).banned_until, s.data)

function inform!(s::ScoredScheduler, rule::AbstractRule, n_matches)
    # println(s.data[rule])

    rd = s.data[rule]
    treshold = rd.match_limit * (rd.weight^rd.times_banned)
    if n_matches > treshold
        ban_length = rd.ban_length * (rd.weight^rd.times_banned)
        rd.times_banned += 1
        rd.banned_until = s.curr_iter + ban_length
        # @info "banning rule $rule until $(rd.banned_until)!"
        return false
    end
    return true
end

function setiter!(s::ScoredScheduler, curr_iter)
    s.curr_iter = curr_iter
end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TypeTrees
##
# Generate a text graphic of Julia modules type tree
##

struct Binding
    mod::Module
    sym::Symbol
end
Binding(tn::TypeName) = Binding(tn.module, tn.name)
Base.isless(a::Binding, b::Binding) = isless(a.sym, b.sym)

# The node type holds the type of the current node and a dict of subtypes
struct TTNode
    typ # ::Type
    subtypes::Dict{Binding, TTNode}

    TTNode(@nospecialize t) = new(t, Dict{Binding, TTNode}())
end

# Add a node to a dict if not added
function add_ttnode(subtypes::Dict{Binding, TTNode}, sname::Binding, tnode::TTNode)
    return get!(subtypes, sname, tnode)
end

function add_ttnode(subtypes::Dict{Binding, TTNode}, sname::Binding, t::Type)
    return get!(subtypes, sname, TTNode(t))
end

# Store a type and its type hierarchy chain
# Recurse till we reach the top level type
function store_type(sname::Binding, t::Union)
    suptype = Union
    tnode = TTNode(t)

    # store unions under Union type
    subtypes = store_type(Binding(suptype.name), suptype)
    add_ttnode(subtypes, sname, tnode)
    store_union(sname, tnode, t)

    # unions are also in a sense related to the types of their components

    return tnode.subtypes
end
function store_union(sname::Binding, tnode::TTNode, @nospecialize t)
    t = Base.unwrap_unionall(t)
    if isa(t, Union)
        store_union(sname, tnode, t.a)
        store_union(sname, tnode, t.b)
    elseif isa(t, DataType)
        binding = Binding(t.name)
        subtypes = store_type(binding, t)
        add_ttnode(subtypes, sname, tnode)
    end
    nothing
end

function store_type(sname::Binding, t::UnionAll)
    suptype = Base.unwrap_unionall(t)
    binding = isa(suptype, DataType) ? Binding(suptype.name) : Binding(Main, string(suptype::Union))
    subtypes = store_type(binding, suptype)
    tnode = add_ttnode(subtypes, sname, t)
    return tnode.subtypes
end

function store_type(sname::Binding, t::DataType)
    suptype = supertype(t)
    subtypes = (suptype != t) ? store_type(Binding(suptype.name), suptype) : types_tree
    tnode = add_ttnode(subtypes, sname, t)
    return tnode.subtypes
end

# examine all symbols in module and store those that are types
function store_all_from(m::Module)
    for s in names(m, true)
        if isdefined(m, s) && !Base.isdeprecated(m, s)
            t = getfield(m, s)
            if isa(t, Type) && t !== Union{}
                store_type(Binding(m, s), t)
            elseif isa(t, Module) && module_name(t) === s && module_parent(t) === m && t !== m
                store_all_from(t)
            end
        end
    end
end

type_props(typ) = ""
type_props(typ::DataType) = string("<<",
                                 typ.abstract    ? " abstract"    : " concrete",
                                 typ.mutable     ? " mutable"     : " immutable",
                                 typ.layout != C_NULL ? string(
                                     Base.datatype_pointerfree(typ) ? " pointerfree" : "",
                                     Base.datatype_haspadding(typ) ? " haspadding" : "",
                                     " nfields:", Core.nfields(typ),
                                     " size:", typ.size,
                                     ", align:", Base.datatype_alignment(typ)) : "",
                                 " >>")

function print_tree(subtypes::Dict{Binding, TTNode}, pfx::String="")
    for b in sort!(collect(keys(subtypes)))
        v = subtypes[b]
        ishidden = unsafe_load(Base.unsafe_convert(Ptr{UInt8}, b.sym)) == UInt8('#')
        if ishidden && supertype(v.typ) === Function
            continue
        end
        if b.mod === Main
            n = string(b.sym)
        elseif !isa(v.typ, DataType) || v.typ.name.module != b.mod || v.typ.name.name != b.sym
            n_io = IOBuffer()
            print(n_io, b.mod, '.', b.sym)
            ua = v.typ
            if isa(ua, UnionAll)
                print(n_io, "{")
                while true
                    print(n_io, ua.var)
                    ua = ua.body
                    if isa(ua, UnionAll)
                        print(n_io, ", ")
                    else
                        break
                    end
                end
                print(n_io, "}")
            end
            n = String(take!(n_io))
        else
            n = string(v.typ)
        end
        if n == string(v.typ)
            println(pfx, "+- ", n, " ", type_props(v.typ))
        else
            println(pfx, "+- ", n, " = ", v.typ, " ", type_props(v.typ))
        end
        v.typ === Function && println(pfx, ".     ## hiding implicit Function subtypes ##")
        print_tree(v.subtypes, pfx * ".  ")
    end
end

# TODO: optionally take module names in command line
# TODO: option to list subtrees of type tree, or other symbol types
const types_tree = Dict{Binding, TTNode}()

store_all_from(Main)

# print_tree(types_tree)

end # module

if !isinteractive()
    TypeTrees.print_tree(TypeTrees.types_tree)
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
## this generates data from the simple supply-demand model
#  q = α1 + α2*p + α3*m + e1
#  q = β1 + β2*p +        e2
## setup
using Econometrics, Plots, DelimitedFiles
# number of obsn
n = 500
# model parameters
α1 = 100.0
α2 = -1.0
α3 = 1.0
β1 = 20.0
β2 = 1.0

## make variables and shocks
# exog var income
m = randn(n,5)
m = sum(m.*m,dims=2) # chi square df=5
# structural shocks
ϵ1 = 1.0*randn(n)
ϵ2 = 1.0*randn(n)
# rf shocks
ν1 = (β2*ϵ1-α2*ϵ2)/(β2-α2)
ν2 = (ϵ1-ϵ2)/(β2-α2)
# rf coefs
π11 = (β2*α1-α2*β1)/(β2-α2)
π21 = β2*α3/(β2-α2)
π12 = (α1-β1)/(β2-α2)
π22 = α3/(β2-α2)

##  generate endog variables computed from rf
q = π11 .+ π21*m + ν1
p = π12 .+ π22*m + ν2
## write and plot data
data = [q p m]
data = sortbyc(data,3)
nn = Int(n/4)
p = scatter(data[1:nn,2], data[1:nn,1], markersize=8, xlabel="price", ylabel="quantity", legend=false)
for i = 2:4
    scatter!(data[(i-1)*nn+1:i*nn,2], data[(i-1)*nn:i*nn,1], markersize=8, show=true)
end
display(p)
#writedlm("data.txt", data)

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
module bcrypt_b64

export bcrypt_b64_decode, bcrypt_b64_encode

import ..base64

const codec = [UInt8(x) for x in ['.'; '/'; 'A':'Z'; 'a':'z'; '0':'9']]

const BcryptEncoding = base64.NewEncoding(codec)

@inline function bcrypt_b64_encode(src::Array{UInt8, 1})
	n = base64.EncodedLen(BcryptEncoding, length(src))
	dst = zeros(UInt8, n)
	base64.Encode!(BcryptEncoding, dst, src)
	while dst[n] == 0x3d
		n -= 1
	end
	return dst[1:n]
end

@inline function bcrypt_b64_decode(src::Array{UInt8, 1}) :: Array{UInt8, 1}
	j = 4-(length(src) % 4)
	for _ = 0:j-1
		push!(src, 0x3d)
	end
	n = base64.DecodedLen(BcryptEncoding, length(src))
	dst = zeros(UInt8, n)
	i = base64.Decode!(BcryptEncoding, dst, src)
	return dst[1:i]
end

end


text:
# import some useful things

import Base: +, -, *, <
import Base: copy, show
# provides polynomial ring and polynomials interface
using AbstractAlgebra
# provides dot product operation
import LinearAlgebra: dot

#-----------------------------------------------------------------------------

struct Signature
    index::Integer
    monom::Any
end

# our module element will store one thing: a tuple of polynomials
struct ModuleElement
    ev::Any
    sgn::Signature
end

# some magic to allow copying ModuleElement
function Base.copy(m::ModuleElement)
    return ModuleElement(deepcopy(m.ev), deepcopy(m.sgn))
end

# some magic to allow copying ModuleElement
function Base.copy(m::Signature)
    return Signature(deepcopy(m.index), deepcopy(m.monom))
end

# some magic to prin2t ModuleElement nicely
Base.show(io::IO, ::MIME"text/plain", m::ModuleElement) = println(io, "($(m.ev), $(m.sgn))")
Base.show(io::IO, m::ModuleElement) = print(io, "($(m.ev), $(m.sgn))")

#-----------------------------------------------------------------------------
# Orderings stuff

function <ₚₒₜ(a::Signature, b::Signature)
    b.index < a.index || (b.index == a.index && a.monom < b.monom)
end

>ₚₒₜ(a::Signature, b::Signature) = b <ₚₒₜ a
<(a::Signature, b::Signature) = a <ₚₒₜ b
Base.isless(a::Signature, b::Signature) = a < b

# compare two module elements w.r.t. pot strategy,
# a nice thing about julia is that we can use some latex capabilities in naming variables
function <ₚₒₜ(a::ModuleElement, b::ModuleElement)
    a.sgn <ₚₒₜ b.sgn
end

>ₚₒₜ(a::ModuleElement, b::ModuleElement) = b <ₚₒₜ a
<(a::ModuleElement, b::ModuleElement) = a <ₚₒₜ b
Base.isless(a::ModuleElement, b::ModuleElement) = a < b

#-----------------------------------------------------------------------------
# Arithmetic operations

*(t, s::Signature) = Signature(s.index, leading_monomial(s.monom*t))
*(s::Signature, t) = t * s

# returns t*m where t is something, and m is a module element
*(t, m::ModuleElement) = m * t
*(m::ModuleElement, t) = ModuleElement(m.ev*t, m.sgn*t)

# group operations
-(a::ModuleElement) = ModuleElement(-a.ev, a.sgn)
+(a::ModuleElement, b::ModuleElement) = ModuleElement(a.ev+b.ev, max(a.sgn, b.sgn))
-(a::ModuleElement, b::ModuleElement) = a + (-b)

function (≃)(a::ModuleElement, b::ModuleElement)
    a.ev == b.ev && a.sgn == b.sgn
end

#-----------------------------------------------------------------------------
# Signatures and evaluation

dot(F, m::ModuleElement) = dot(m, F)
dot(m::ModuleElement, F) = sum(m.mtuple .* F)

# evaluation hom
v(F, m::ModuleElement) = dot(F, m)
v(F, M) = [v(F, m) for m in M]

issyzygy(h::ModuleElement) = iszero(h.ev)


text:
# Autogenerated wrapper script for SeqGen_jll for x86_64-linux-musl
export seqgen

JLLWrappers.@generate_wrapper_header("SeqGen")
JLLWrappers.@declare_executable_product(seqgen)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        seqgen,
        "bin/seq-gen",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()