text:
"""
```
euler_equation_error(m, cₜ, logSDFxR, 𝔼_quadrature, zₜ = m.z;
    c_init = NaN, return_soln = false, kwargs...)
euler_equation_error(m, cₜ, logSDFxR, 𝔼_quadrature, shock_matrix, zₜ = m.z;
    c_init = NaN, summary_statistic = x -> norm(x, Inf), burnin = 0,
    return_soln = false, kwargs...)
```
calculates standard Euler equation errors, as recommended by Judd (1992).
The first method calculates the error at some state `zₜ`, which defaults
to the stochastic steady state. The second method simulates the state
vector from an initial state `zₜ` (defaults to stochastic steady state)
given a sequence of drawn shocks, evaluates the Euler equation errors,
and returns some summary statistic of the errors specified by the keyword
`summary_statistic`.

The Euler equation is

```math
\\begin{aligned}
0 = \\log \\mathbb{E}_t \\exp(m_{t + 1} + r_{t + 1}) = \\log \\mathbb{E}_t[M_{t + 1} R_{t + 1}],
\\end{aligned}
```

where ``m_{t + 1} = \\log(M_{t + 1})`` is the log stochastic discount factor and ``r_{t + 1} = \\log(R_{t + 1})``
is the risk free rate.

### Inputs
- `m::RiskAdjustedLinearization`: A solved instance of a risk-adjusted linearization
- `cₜ::Function`: a function of `(m, zₜ)` that calculates consumption at state `zₜ`, given the
    state-space representation implied by `m`.
- `logSDFxR::Function`: a `Function` evaluating ``m_{t + 1} + r_{t + 1}``. The `Function` must
    take as input `(m, zₜ, εₜ₊₁, c)`, where `m` is a `RiskAdjustedLinearization`,
    `zₜ` is a state vector at which to evaluate, `εₜ₊₁` is a draw from the distribution
    of exogenous shocks, and `c` is a guess for consumption at `zₜ` implied by
    the conditional expectation in the Euler equation when calculated with a quadrature rule.
    Note that `c` can be either the consumption level or some transformation (e.g. log consumption),
    but the user should be consistent in the definition of the `cₜ` function with the guess `c`,
    i.e. both should return the same transformation of consumption (e.g. both should return the level).
- `𝔼_quadrature::Function`: a quadrature rule whose single input is a `Function` with a single
    input, which is a shock `εₜ₊₁`.
- `zₜ::AbstractVector`: a state at which to evaluate the Euler equation error
- `shock_matrix::Abstractmatrix`: a `Nε × T` matrix of shocks drawn from the distribution of exogenous shocks.

### Keywords
- `c_init::Number`: an initial guess to be used when solving the "true" consumption policy using
    quadrature. The default is the consumption policy according to the `RiskAdjustedLinearization`
- `summary_statistic::Function`: a `Function` used to compute a summary statistic from the
    ergodic set of Euler equation errors. The default is the maximum absolute error.
- `burnin::Int`: number of periods to drop as burn-in
- `return_soln::Bool`: if true, return the solution to the nonlinear equation isntead of the error
- `kwargs`: Any keyword arguments for `nlsolve` can be passed, too, e.g. `ftol` or `autodiff`
    since `nlsolve` is used to calculate the "true" consumption policy.
"""
function euler_equation_error(m::RiskAdjustedLinearization, cₜ::Function, logSDFxR::Function, 𝔼_quadrature::Function,
                              zₜ::AbstractVector = m.z; c_init::Number = NaN,
                              return_soln::Bool = false, kwargs...)

    # Compute expected consumption according to RAL
    c_ral = cₜ(m, zₜ)

    # Compute implied consumption according to the quadrature rule
    out = nlsolve(c -> [log(𝔼_quadrature(εₜ₊₁ -> exp(logSDFxR(m, zₜ, εₜ₊₁, c[1]))))], [isnan(c_init) ? c_ral : c_init];
                  kwargs...)

    if out.f_converged
        c_impl = out.zero[1]
    else
        error("Failed to solve implied consumption.")
    end

    # Return error in unit-free terms
    if return_soln
        return c_impl
    else
        return (c_ral - c_impl) / c_ral
    end
end

function euler_equation_error(m::RiskAdjustedLinearization, cₜ::Function, logSDFxR::Function, 𝔼_quadrature::Function,
                              shock_matrix::AbstractMatrix, zₜ::AbstractVector = m.z; c_init::Number = NaN,
                              summary_statistic::Function = x -> norm(x, Inf), burnin::Int = 0,
                              return_soln::Bool = false, kwargs...)

    # Set up
    T = size(shock_matrix, 2)

    # Simulate states
    states, _ = simulate(m, T, shock_matrix, zₜ)

    # Compute implied consumption according to the quadrature rule for each state
    # and expected consumption according to RAL
    err = [euler_equation_error(m, cₜ, logSDFxR, 𝔼_quadrature, (@view states[:, t]); c_init = c_init,
                                return_soln = return_soln, kwargs...) for t in (burnin + 1):T]

    # Return error in unit-free terms
    if return_soln
        return err
    else
        return summary_statistic(err)
    end
end

# n-period ahead euler equation error obtained by the
# Compute cₜ⁽¹⁾ by solving
# 0 = log Eₜ[exp(mₜ₊₁(cₜ₊₁⁽⁰⁾, cₜ⁽¹⁾) + rₜ)],
# where cₜ₊₁⁽⁰⁾ is just the consumption function obtained from
# an RAL evaluated at a state in t + 1 (i.e., given zₜ and a
# a draw of shocks εₜ₊₁, one can calculate zₜ₊₁ and then
# evaluate cₜ₊₁⁽⁰⁾ = c + Ψ (zₜ₊₁ - z)).
# This approach implicitly defines a mapping cₜ⁽¹⁾ at every
# state zₜ since the expectation must be calculated conditional
# on the state zₜ. It follows that if cₜ⁽²⁾ solves
# 0 = log Eₜ[exp(log(β) - γ * (cₜ₊₁⁽¹⁾ - cₜ⁽²⁾) + rₜ)],
# then we calculate cₜ₊₁⁽¹⁾ by solving the first equation
# given an initial state zₜ₊₁. Since solving the first equation
# requires pushing the initial state forward by one period,
# calculating cₜ₊₁⁽¹⁾ means solving
# 0 = log Eₜ₊₁[exp(mₜ₊₂(cₜ₊₂⁽⁰⁾, cₜ₊₁⁽¹⁾) + rₜ₊₁)],
# hence cₜ⁽²⁾ is indeed a 2-period ahead Euler equation error
# b/c the recursion implicitly involves calculating what
# expected consumption is 2 periods ahead. Since
# we approximate the expectation with quadrature,
# this approach is recursive in that, we select some
# initial state zₜ, compute a set of {zₜ₊₁} implied
# by zₜ and draws εₜ₊₁, compute a new set {zₜ₊₂} for each
# zₜ₊₁ and draws εₜ₊₂ (for each zₜ₊₁), and so on until
# we have computed zₜ₊ₙ

"""
```
dynamic_euler_equation_error(m, cₜ, logSDFxR, 𝔼_quadrature, endo_states, n_aug,
    shock_matrix, zₜ = m.z; c_init = NaN, summary_statistic = x -> norm(x, Inf),
    burnin = 0, return_soln = false, kwargs...)
```
calculates dynamic Euler equation errors, as proposed in Den Haan (2009).
The Euler equation is

``math
\\begin{aligned}
0 = \\log \\mathbb{E}_t \\exp(m_{t + 1} + r_{t + 1}) = \\log \\mathbb{E}_t[M_{t + 1} R_{t + 1}],
\\end{aligned}
``

where ``m_{t + 1} = \\log(M_{t + 1})`` is the log stochastic discount factor and ``r_{t + 1} = \\log(R_{t + 1})``
is the risk free rate.

The dynamic errors are computed according the following algorithm.

1. Simulate according to the risk-adjusted linearization time series for the state variables
2. Using the time series from 1, compute time series for consumption and
    some state variable (usually capital) that can ensure budget constraints hold and markets
    clear when computing consumption by applying quadrature.
3. Generate a second "implied" time series for consumption and the "capital" state variable,
    starting from the same initial state as 2. Repeat the following steps at each time period.
    (i)  Compute the conditional expectation in the Euler equation using quadrature to
         obtain implied consumption.
    (ii) Use budget constraint/market-clearing to compute implied capital.

By default, `dynamic_euler_equation_error` returns some summary statistic of the errors
specified by the keyword `summary_statistic`.

### Inputs
- `m::RiskAdjustedLinearization`: A solved instance of a risk-adjusted linearization
- `cₜ::Function`: a function of `(m, zₜ)` that calculates consumption at state `zₜ`, given the
    state-space representation implied by `m`.
- `logSDFxR::Function`: a `Function` evaluating ``m_{t + 1} + r_{t + 1}``. The `Function` must
    take as input `(m, zₜ, εₜ₊₁, cₜ)`, where `m` is a `RiskAdjustedLinearization`,
    `zₜ` is a state vector at which to evaluate, `εₜ₊₁` is a draw from the distribution
    of exogenous shocks, and `cₜ` is the a guess for consumption at `zₜ` implied by
    the conditional expectation in the Euler equation when calculated with a quadrature rule.
- `𝔼_quadrature::Function`: a quadrature rule whose single input is a `Function` with a single
    input, which is a shock `εₜ₊₁`.
- `endo_states::Function`: augments the state variables in the risk-adjusted linearization,
    usually with one additional variable, which represents capital or assets.
- `n_aug::Int`: number of extra state variables added by `endo_states` (usually 1).
- `zₜ::AbstractVector`: a state at which to evaluate the Euler equation error
- `shock_matrix::Abstractmatrix`: a `Nε × T` matrix of shocks drawn from the distribution of exogenous shocks.

### Keywords
- `c_init::Number`: an initial guess to be used when solving the true consumption policy using
    quadrature. The default is the consumption policy according to the `RiskAdjustedLinearization`
- `summary_statistic::Function`: a `Function` used to compute a summary statistic from the
    ergodic set of Euler equation errors. The default is the maximum absolute error.
- `burnin::Int`: number of periods to drop as burn-in
- `kwargs`: Any keyword arguments for `nlsolve` can be passed, too, e.g. `ftol` or `autodiff`
    since `nlsolve` is used to calculate the "true" consumption policy.
"""
function dynamic_euler_equation_error(m::RiskAdjustedLinearization, cₜ::Function, logSDFxR::Function,
                                      𝔼_quadrature::Function, endo_states::Function, n_aug::Int,
                                      shock_matrix::AbstractMatrix, z₀::AbstractVector = m.z;
                                      c_init::Number = NaN, summary_statistic::Function = x -> norm(x, Inf),
                                      burnin::Int = 0, return_soln::Bool = false, kwargs...)

    # Set up
    T = size(shock_matrix, 2)
    c_impl = Vector{eltype(shock_matrix)}(undef, T)

    # Simulate states and calculate consumption according to RAL
    states, _  = simulate(m, T, shock_matrix, z₀)
    c_ral      = [cₜ(m, (@view states[:, t])) for t in 1:T]
    orig_i     = 1:size(states, 1)

    # Additional set up
    endo_states_impl = similar(states, length(orig_i) + n_aug, T)
    endo_states_ral  = similar(endo_states_impl)

    # For each state, calculate conditional expectation using quadrature rule
    # and compute the implied states
    out = nlsolve(c -> [log(𝔼_quadrature(εₜ₊₁ -> exp(logSDFxR(m, (@view states[:, 1]), εₜ₊₁, c[1]))))], [isnan(c_init) ? c_ral[1] : c_init];
                  kwargs...) # Do period 1 separately b/c needed to initialize endo_states_impl
    if out.f_converged
        c_impl[1] = out.zero[1]
    else
        error("Failed to solve implied consumption in period 1 of $T.")
    end
    endo_states_impl[:, 1] = endo_states(m, (@view states[:, 1]), z₀, c_impl[1])
    endo_states_ral[:, 1]  = endo_states(m, (@view states[:, 1]), z₀, c_ral[1])

    for t in 2:T
        out = nlsolve(c -> [log(𝔼_quadrature(εₜ₊₁ -> exp(logSDFxR(m, (@view states[:, t]), εₜ₊₁, c[1]))))], [isnan(c_init) ? c_ral[t] : c_init];
                      kwargs...)
        if out.f_converged
            c_impl[t] = out.zero[1]
        else
            error("Failed to solve implied consumption in period $t of $T.")
        end
        endo_states_impl[:, t] = endo_states(m, (@view states[:, t]), (@view endo_states_impl[orig_i, t - 1]), c_impl[t])
        endo_states_ral[:, t]  = endo_states(m, (@view states[:, t]), (@view endo_states_ral[orig_i, t - 1]),  c_ral[t])
    end

    # Calculate the errors
    if return_soln
        return c_ral[(burnin + 1):end], c_impl[(burnin + 1):end], endo_states_ral[(burnin + 1):end], endo_states_impl[(burnin + 1):end]
    else
        return summary_statistic(((@view c_ral[(burnin + 1):end]) - (@view c_impl[(burnin + 1):end])) ./ (@view c_ral[(burnin + 1):end])), summary_statistic(vec((@view endo_states_ral[:, (burnin + 1):end]) - (@view endo_states_impl[:, (burnin + 1):end])) ./ vec((@view endo_states_ral[:, (burnin + 1):end])))
    end
end


text:
"""
Computes the residuals of the arbitrage equations. The general form of the arbitrage equation is

    `0 = E_t [f(m, s, x, M, S, X; p)]`

where `m` are current exogenous variables, `s` are current states,
`x` are current controls, `M` are next period's exogenous variables, `S` are next period's states, `X` are next period's controls, and `p` are the model parameters. This function evaluates the right hand side of the arbitrage equation for the given inputs.

If the list of current controls `x` is provided as a two-dimensional array (`ListOfPoints`), it is transformed to a one-dimensional array (`ListOfListOfPoints`).

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `dprocess`: Discretized exogenous process.
* `s::ListOfPoints`: List of state variable values.
* `x::ListOfListOfPoints`: List of control variable values associated with each exogenous shock.
* `p::Vector{Float64}`: Model parameters.
* `dr`: Current guess for the decision rule.
# Returns
* `res`: Residuals of the arbitrage equation associated with each exogenous shock.
"""
function euler_residuals_ti(model, dprocess::AbstractDiscretizedProcess,s::ListOfPoints{d}, x::Vector{ListOfPoints{n_x}}, p::SVector, dr) where n_x where d

    N = length(s)
    # TODO: allocate properly...
    res = deepcopy(x)
    for i_m=1:length(res)
        res[i_m][:] *= 0.0
    end
    for i in 1:size(res, 1)
        m = node(Point, dprocess, i)
        for (w, M, j) in get_integration_nodes(Point, dprocess,i)
            # Update the states
            S = Dolo.transition(model, m, s, x[i], M, p)
            X = dr(i, j, S)
            res[i] += w*Dolo.arbitrage(model, m, s, x[i], M, S, X, p)
        end
    end
    return res
end

struct TimeIterationLog
    header::Array{String, 1}
    keywords::Array{Symbol, 1}
    entries::Array{Any, 1}
end

function TimeIterationLog()
    header = [ "It", "ϵₙ", "ηₙ=|xₙ-xₙ₋₁|", "λₙ=ηₙ/ηₙ₋₁", "Time", "Newton steps"]
    keywords = [:it, :err, :gain, :time, :nit]
    TimeIterationLog(header, keywords, [])
end

function append!(log::TimeIterationLog; verbose=true, entry...)
    d = Dict(entry)
    push!(log.entries, d)
    verbose && show_entry(log, d)
end

function initialize(log::TimeIterationLog; verbose=true)
    verbose && show_start(log)
end

function finalize(log::TimeIterationLog; verbose=true)
    verbose && show_end(log)
end

function show(log::TimeIterationLog)
    show_start(log)
    for entry in log.entries
        show_entry(log,entry)
    end
    show_end(log)
end

function show_start(log::TimeIterationLog)
    println(repeat("-", 66))
    @printf "%-6s%-16s%-16s%-16s%-16s%-5s\n" "It" "ϵₙ" "ηₙ=|xₙ-xₙ₋₁|" "λₙ=ηₙ/ηₙ₋₁" "Time" "Newton steps"
    println(repeat("-", 66))
end

function show_entry(log::TimeIterationLog, entry)
    it = entry[:it]
    epsilon = entry[:epsilon]
    err = entry[:err]
    gain = entry[:gain]
    time = entry[:time]
    nit = entry[:nit]
    @printf "%-6i%-16.2e%-16.2e%-16.2e%-16.2e%-5i\n" it epsilon err gain time nit
end

function show_end(log::TimeIterationLog)
    println(repeat("-", 66))
end

###

mutable struct IterationTrace
    trace::Array{Any,1}
end

mutable struct TimeIterationResult
    dr::AbstractDecisionRule
    iterations::Int
    complementarities::Bool
    dprocess::AbstractDiscretizedProcess
    x_converged::Bool
    x_tol::Float64
    err::Float64
    log::TimeIterationLog
    trace::Union{Nothing,IterationTrace}
end

converged(r::TimeIterationResult) = r.x_converged
function Base.show(io::IO, r::TimeIterationResult)
    @printf io "Results of Time Iteration Algorithm\n"
    @printf io " * Complementarities: %s\n" string(r.complementarities)
    @printf io " * Discretized Process type: %s\n" string(typeof(r.dprocess))
    @printf io " * Decision Rule type: %s\n" string(typeof(r.dr))
    @printf io " * Number of iterations: %s\n" string(r.iterations)
    @printf io " * Convergence: %s\n" converged(r)
    @printf io "   * |x - x'| < %.1e: %s\n" r.x_tol r.x_converged
end

"""
Computes a global solution for a model via backward time iteration. The time iteration is applied to the residuals of the arbitrage equations.

If the initial guess for the decision rule is not explicitly provided, the initial guess is provided by `ConstantDecisionRule`.
If the stochastic process for the model is not explicitly provided, the process is taken from the default provided by the model object, `model.exogenous`

# Arguments
* `model::NumericModel`: Model object that describes the current model environment.
* `process`: The stochastic process associated with the exogenous variables in the model.
* `init_dr`: Initial guess for the decision rule.
# Returns
* `dr`: Solved decision rule.
"""
function time_iteration(model::Model, dprocess::AbstractDiscretizedProcess,
                        grid, init_dr;
                        verbose::Bool=true,
                        maxit::Int=500, tol_η::Float64=1e-7, trace::Bool=false, maxbsteps=5, dampen=1.0,
                        solver=Dict(), complementarities::Bool=true)

    if get(solver, :type, :__missing__) == :direct
        return time_iteration_direct(
            model, dprocess, grid, init_dr; verbose=verbose,
            maxit=maxit, tol_η=tol_η
        )
    end

    endo_nodes = nodes(ListOfPoints, grid)
    N = n_nodes(grid)
    n_s_endo = size(endo_nodes, 2)
    n_s_exo = n_nodes(dprocess)

    # initial guess
    # number of smooth decision rules
    n_m = max(n_s_exo, 1)
    p = SVector(model.calibration[:parameters]...)

    x0 = [init_dr(i, endo_nodes) for i=1:n_m]

    ti_trace = trace ? IterationTrace([x0]) : nothing

    n_x = length(model.calibration[:controls])
    if complementarities == true
        x_lb = [controls_lb(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
        x_ub = [controls_ub(model, node(Point,dprocess,i),endo_nodes,p) for i=1:n_m]
        BIG = 100000
        for i=1:n_m
          for n=1:N
            x_lb[i][n] = max.(x_lb[i][n],-BIG)
            x_ub[i][n] = min.(x_ub[i][n], BIG)
          end
        end
    end

    # create decision rule (which interpolates x0)
    dr = CachedDecisionRule(dprocess, grid, x0)

    steps = 0.5.^collect(0:maxbsteps)

    err_0 = NaN
    err = 1.0

    log = TimeIterationLog()
    initialize(log, verbose=verbose)

    it = 0
    while it<maxit && err>tol_η

        it += 1

        t1 = time_ns()

        set_values!(dr, x0)
        fobj(u) = euler_residuals_ti(model, dprocess, endo_nodes, u, p, dr)

        tt = euler_residuals_ti(model, dprocess, endo_nodes, x0, p, dr)

        if complementarities
            res = newton(fobj, x0, x_lb, x_ub; solver...)
        else
            res = newton(fobj, x0; solver...)
        end

        x1 = res.solution
        nit = res.iterations
        epsil = res.errors[1]

        trace && push!(ti_trace.trace, x1)

        err = maxabs(x1-x0)
        x0 = (1-dampen)*x0 + dampen*x1

        gain = err / err_0
        err_0 = err

        elapsed = time_ns() - t1

        append!(log; verbose=verbose, it=it, epsilon=epsil, err=err, gain=gain, time=elapsed, nit=nit)
    end

    finalize(log, verbose=verbose)

    # TODO: somehow after defining `fobj` the `dr` object gets `Core.Box`ed
    #       making the return type right here non-inferrable.

    converged = err < tol_η
    res = TimeIterationResult(dr.dr, it, true, dprocess, converged, tol_η, err, log, ti_trace)

end

# get grid for endogenous
function time_iteration(model, dprocess, dr; grid=Dict(), kwargs...)
    grid = get_grid(model, options=grid)
    return time_iteration(model, dprocess, grid, dr;  kwargs...)
end

# get stupid initial rule
function time_iteration(model, dprocess::AbstractDiscretizedProcess; grid=Dict(), kwargs...)

    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return time_iteration(model, dprocess, init_dr;  grid=grid, kwargs...)
end

function time_iteration(model, init_dr; grid=Dict(), kwargs...)
    dprocess = discretize( model.exogenous )
    return time_iteration(model, dprocess, init_dr; grid=grid, kwargs...)
end

function time_iteration(model; grid=Dict(), kwargs...)
    dprocess = discretize( model.exogenous )
    init_dr = ConstantDecisionRule(model.calibration[:controls])
    return time_iteration(model, dprocess, init_dr; grid=grid, kwargs...)
end

ext: jl


text:
module MusicFadeTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Music Fade" => Ahorn.EntityPlacement(
        Maple.MusicFadeTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.MusicFadeTrigger)
    return Dict{String, Any}(
        "direction" => Maple.music_fade_trigger_directions
    )
end

end


text:
"""
    buildW(G,hind,hcovers,pr)

Precomputes
"""
function buildW{D, BF}(G::NGrid{D, BF}, hind, hcovers, pr = 1:length(G))
    bf = (BF == Linear ? cc_bf_l : cc_bf_q)
    nc 		= size(G.covers, 1)

    coversid = zeros(Int32, length(G))
    for i = 1:length(hind)
        coversid[i] = findfirst(hcovers, hind[i])
    end
    mL      = maximum(G.L)
    J 		= ones(Int, mL + 1, D)
    B 		= zeros(mL + 1, D)
    for i = pr
        G.IDs[i] = Int[]
        G.Bs[i] = Float64[]
        for d = 1:D
            for l = 1:mL + 1
                j 	= clamp(round(Int, G.grid[i, d] * (dM(l)) + 1 / 2), 1, dM(l))
                B[l, d] 	= bf(G.grid[i, d], cc_dg(l, j), Int16(M(l)))
                J[l, d]  = j
            end
        end
        for ii = 1:coversid[i]
            b  = B[G.covers[ii, D], D] * B[G.covers[ii, 1], 1]
            id1 = J[G.covers[ii, D], D] - 1
            for d = D - 1:-1:2
                b *= B[G.covers[ii, d], d]
                id1 = id1 * G.covers_dM[ii, d] + (J[G.covers[ii, d], d] - 1)
            end
            id1 = (J[G.covers[ii, 1], 1] - 1) + G.covers_dM[ii, 1] * id1 + 1 + G.covers_loc[ii] - 1
            if b > 0
                push!(G.IDs[i], id1)
                push!(G.Bs[i], id1 == i ? b : -b)
            end
        end
    end
    return
end

# This section constructs specialised functions for grids of dimensions 2-12
Lbj = quote
@inbounds for d in dr1
    J[2, d] = (x[i, d] > 0.5)
    B[2, d] = (1.0 - 2 * abs(x[i, d] - J[2, d]))
end
@inbounds for l = 3:mL
    dm = 2^(l - 2)
    m  = dm * 2
    for d in dr1
        j 	= clamp(round(Int, x[i, d] * dm + 1 / 2), 1, dm)
        xij = (2j - 1) / m
        dx = (1.0 - (m) * abs(x[i, d] - xij))
        B[l, d] 	= dx
        J[l, d]  = j - 1
    end
end
end

Qbj = quote
@inbounds for d in dr1
    J[2, d] = (x[i, d] > 0.5)
    B[2, d] = 1.0 - (2 * (x[i, d] - J[2, d]))^2
end
@inbounds for l = 3:mL
    dm = 2^(l - 2)
    m  = dm * 2
    for d in dr1
        j 	= clamp(round(Int, x[i, d] * dm + 1 / 2), 1, dm)
        xij = (2j - 1) / m
        dx = 1.0 - (m * (x[i, d] - xij))^2
        B[l, d]  = dx
        J[l, d]  = j - 1
    end
end
end

for b in [(Linear, Lbj), (Quadratic, Qbj)]
    for D = 2:15
        coverloop = :(@inbounds for ii in nc
            b   = B[G.covers[ii, $D], $D] * B[G.covers[ii, 1], 1]
            id1 = J[G.covers[ii, $D], $D]
        end)
        for d in D - 1:-1:2
            push!(coverloop.args[2].args[2].args, :(b  *= B[G.covers[ii, $d], $d]))
            push!(coverloop.args[2].args[2].args, :(id1 = id1 * G.covers_dM[ii, $d] + J[G.covers[ii, $d], $d]))
        end

        push!(coverloop.args[2].args[2].args, :(id1 = id1 * G.covers_dM[ii, 1] + G.covers_loc[ii] + J[G.covers[ii, 1], 1];yi += b * w[id1]))

        f = :(function interpolate(G::NGrid{$D, $(b[1])}, A::Vector{Float64}, xi::Array{Float64, 2}, y = zeros(Float64, size(xi, 1)))
            w         = getW(G, A)
            x         = nXtoU(xi, G.bounds)
            nx        = size(x, 1)
            cr         = 1:size(G.covers, 1)
            dr1 = 1:$D
            mL      = maximum(G.L) + 1
            J         = zeros(Int, mL, $D)
            B         = ones(mL, $D)

            @threadsfixed [J, B] for i = 1:nx
                $(b[2])
                yi = 0.0
                @inbounds for ii in cr
                    $(Expr(:block,
                        :(b   = B[G.covers[ii, $D], $D] * B[G.covers[ii, 1], 1]),
                        (:(b  *= B[G.covers[ii, $d], $d]) for d in D - 1:-1:2)...,
                        :(id1 = J[G.covers[ii, $D], $D]),
                        (:(id1 = id1 * G.covers_dM[ii, $d] + J[G.covers[ii, $d], $d]) for d in D - 1:-1:2)...,
                        :(id1 = id1 * G.covers_dM[ii, 1] + G.covers_loc[ii] + J[G.covers[ii, 1], 1]),
                        :(yi += b * w[id1])
                        ))
                end
                y[i] = yi
            end
            y
        end)

        subs!(f, :coverloop => coverloop)
        eval(current_module(), f)
    end
end

mutable struct dimdef{n} end

for b in [(Linear, Lbj), (Quadratic, Qbj)]
    for D = 2:15
        for adim = 1:8
            coverloop = quote
                fill!(yi, 0.0)
                @inbounds for ii in nc
                    b   = B[G.covers[ii, $D], $D] * B[G.covers[ii, 1], 1]
                    id1 = J[G.covers[ii, $D], $D]
                end
            end
            for d in D - 1:-1:2
                push!(coverloop.args[4].args[2].args[2].args, :(b  *= B[G.covers[ii, $d], $d]))
                push!(coverloop.args[4].args[2].args[2].args, :(id1 = id1 * G.covers_dM[ii, $d] + J[G.covers[ii, $d], $d]))
            end
            push!(coverloop.args[4].args[2].args[2].args, :(id1 = id1 * G.covers_dM[ii, 1] + G.covers_loc[ii] + J[G.covers[ii, 1], 1]))
            for d = 1:adim
                push!(coverloop.args[4].args[2].args[2].args, :(yi[$d] += b * w[id1, $d]))
                push!(coverloop.args, :(@inbounds y[i, $d] = yi[$d]))
            end

            f = :(function interpolate(G::NGrid{$D, $(b[1])}, A::Array{Float64, 2}, xi::Array{Float64, 2}, nA::dimdef{$adim}, y::Array{Float64, 2})
                w         = getW(G, A)
                x         = nXtoU(xi, G.bounds)
                nx        = size(x, 1)

                nc         = 1:size(G.covers, 1)
                dr1 = 1:$D
                mL      = maximum(G.L) + 1
                J         = zeros(Int, mL, $D)
                B         = ones(mL, $D)
                yi        = zeros($adim)

                @threadsfixed [J, B, yi] for i = 1:nx
                    $(b[2])
                    coverloop
                end
                y
            end)

            subs!(f, :coverloop => coverloop)
            eval(current_module(), f)
        end
    end
end

interpolate(G::NGrid, A::Array{Float64, 2}, x::Array{Float64, 2}, y = zeros(size(x, 1), size(A, 2))) = interpolate(G, A, x, dimdef{size(A, 2)}(), y)

function getW(G::NGrid, A::Vector{Float64})
    w = copy(A)
    for i = 1:length(G)
        for l ∈ 1:length(G.Bs[i]) - 1
            @inbounds w[i] +=  G.Bs[i][l] * w[G.IDs[i][l]]
        end
    end
    return w
end

function getW(G::NGrid, A::Array{Float64, 2})
    nA = size(A, 2)
    w = copy(A)
    for i = 1:length(G)
        for l ∈ 1:length(G.Bs[i]) - 1
            for a = 1:nA
                @inbounds w[i, a] +=  G.Bs[i][l] * w[G.IDs[i][l], a]
            end
        end
    end
    return w
end


text:
using Random: AbstractRNG

export invert,
       primdist,
       distapply,
       StdNormal,
       StdUniform,
       UniformInt,
       PrimClass,
       PrimRandVar

"Primitive random variable class"
abstract type PrimClass <: AbstractClass end
Var.isclass(class::PrimClass) = true
Var.traitvartype(class::PrimClass) = Var.TraitIsClass()
Var.traitvartype(class::Type{<:PrimClass}) = Var.TraitIsClass()

"An primitive random variable"
const PrimRandVar{F <: PrimClass, ID} = Member{F, ID}

(primclass::PrimClass)(id, ω) = Member(id, primclass)(ω)

struct StdNormal{T} <: PrimClass end
Base.eltype(::Type{StdNormal{T}}) where T = T
Base.rand(rng::AbstractRNG, ::StdNormal{T}) where {T} = randn(rng, T)

"StdUniform Uniformly Distributed over [0,1]"
struct StdUniform{T} <: PrimClass end
Base.eltype(::Type{StdUniform{T}}) where T = T
Base.rand(rng::AbstractRNG, ::StdUniform{T}) where T = rand(rng, T)

"Uniformly Distributed ove Integers"
struct UniformInt{T} <: PrimClass end
Base.eltype(::Type{UniformInt{T}}) where T = T
Base.rand(rng::AbstractRNG, ::UniformInt{T}) where T = rand(rng, T)

"""`primdist(d::Distribution)``
Primitive (parameterless) distribution that `d` is defined in terms of"""
function primdist end

"""`invert(d::Distribution, val)`
If output of `val` is `val` what must its primitives have been?`"""
function invert end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
export fasta2matrix, letter2num, num2letter, matrix2fasta

##############################################################
##############################################################

"""
  function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)

  Reads an aligned fasta file and returns a matrix with amino acids converted
  to numbers following the convention A --> 1, C --> 2, ... '-' --> 21.

    "filename": full path of the fasta file.
    "max_gap_fraction": to be removed.

"""

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

function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)
    f = FastaReader(filename)

    max_gap_fraction = Float64(max_gap_fraction)

    # pass 1

    seqs = Int[]
    inds = Int[]
    fseqlen = 0

    for (name, seq) in f
        ngaps = 0
        if f.num_parsed == 1
            ls = length(seq)
            resize!(inds, ls)
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    fseqlen += 1
                    inds[fseqlen] = i
                    c == '-' && (ngaps += 1)
                end
            end
        else
            ls = length(seq)
            ls == length(inds) || error("inputs are not aligned")
            tstfseqlen = 0
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    tstfseqlen += 1
                    inds[tstfseqlen] == i || error("inconsistent inputs")
                    c == '-' && (ngaps += 1)
                end
            end
            tstfseqlen == fseqlen || error("inconsistent inputs")
        end
        ngaps / fseqlen <= max_gap_fraction && push!(seqs, f.num_parsed)
    end

    length(seqs) > 0 || error("Out of $(f.num_parsed) sequences, none passed the filter (max_gap_fraction=$max_gap_fraction)")

    # pass 2
    Z = Array{Int16}(undef, fseqlen, length(seqs))
    seqid = 1
    for (name, seq) in f
        seqs[end] < f.num_parsed && break
        seqs[seqid] == f.num_parsed || continue
        for i = 1:fseqlen
            c = seq[inds[i]]
            Z[i, seqid] = letter2num(c)
        end
        seqid += 1
    end
    @assert seqid == length(seqs) + 1

    close(f)

    size(Z)[2] == 1 && return Z[:, 1]
    return Int16.(Z')
end

##############################################################
"""
    letter2num(c::Union{Char,UInt8})

    Takes as input a char (representing an amino acid)
    and returns its conventional number representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

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

let alphabet = [ 1, 21, 2, 3, 4, 5, 6, 7, 8, 21, 9, 10, 11, 12, 21, 13, 14, 15, 16, 17, 21, 18, 19, 21, 20]
               # A, B,  C, D, E, F, G, H, I, J,  K, L,  M,  N,  O,  P,  Q,  R,  S,  T,  U,  V,  W,  X,  Y
    global letter2num
    function letter2num(c::Union{Char,UInt8})
        i = UInt8(c) - 0x40
        1 <= i <= 25 && return alphabet[i]
        return 21
     end
end

##############################################################
"""
    num2letter(i::Integer)

    Takes as input an integer (representing an amino acid)
    and returns its conventional character representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

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

let alphabet = ["A", "C", "D", "E", "F", "G", "H", "I",  "K", "L",  "M",  "N", "P",  "Q",  "R",
"S",  "T", "V",  "W",  "Y"]
    global num2letter
    function num2letter(i :: Integer)
        1 <= i <= 20 && return alphabet[i]
        return "-"
    end
end

##############################################################
"""
    vec2string(v::Array{:<Integer, 1})

    Takes as input a vector of integers (representing an amino acids)
    and returns a list of characters, the corresponding amino acids.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

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

function vec2string(v::Array{<:Integer, 1})
    s = ""
    for i in v
        s = s*num2letter(i)
    end
    return s
end

##############################################################
"""
	string2vec(s::AbstractString)

    Takes as input a string (representing amino acids)
    and returns a vector of numbers, the corresponding number representation.
    In this case with the convention "1 2 .. 21" <==> "A C .. -"
"""

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

function string2vec(s::AbstractString)
    v = Vector{Int8}(undef, length(s))
    for (i, l) in enumerate(s)
        v[i] = letter2num(l)
    end
    return v
end

##############################################################
"""
	write_single_muts_MSA(path_wt::AbstractString, fitness_wt::Array{Union{Missing, Float64}},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")
	write_single_muts_MSA(wt::Array{<:Integer, 1}, fitness_wt::Array{Union{Missing, Float64}},
			dir_out::AbstractString; wt_name::AbstractString = "WT")
	write_single_muts_MSA(path_wt::AbstractString, fitness_wt::Array{Any, 2},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")
	write_single_muts_MSA(wt::Array{<:Integer, 1}, fitness_wt::Array{Any, 2},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")

    USE:
    Given a sequence, it returns all single mutations in a MSA fashon.

    INPUT:
    "fitness_wt": if it is a 1d vector is must contain missings for missing mutations
    			  and should cover only the non gapped positions

    			  if it is 2d vector the first column must show the position, the second the mutation
    			  and the third the fitness value

    OUTPUT: writes a .fasta and a .fit file in the "dir_out" folder

"""

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

function write_single_muts_MSA(path_wt::AbstractString, fitness_wt::Array{Union{Missing, Float64}},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")

	wt = fasta2matrix(path_wt)
	dir_out == "" && (  dir_out = (*).(split(path_wt, "/")[1:end-1])  )
	path_MSA_out = joinpath(dir_out, "single_muts_MSA_$(wt_name).fasta")
	path_fit_out = joinpath(dir_out, "single_muts_fitness_$(wt_name).fit")
	vec_fit = []
	k = 0
	FastaWriter(path_MSA_out, "w") do file
	    writeentry(file, "$k | $(wt_name)", vec2string(wt))
	    for i in [pos for (pos, amino) in enumerate(wt) if amino!= 21]
	        mutant = copy(wt)
	        for amino in 1:20
	            mutant[i] = amino
	            new_amino = num2letter(mutant[i])
	            old_amino = num2letter(wt[i])
	            if !ismissing(fitness_wt[(i-1)*20 + amino]) && (amino != wt[i])
	                k+=1
	                append!(vec_fit, fitness_wt[(i-1)*20 + amino])
	                writeentry(file, "$k | $(old_amino)$(i)$(new_amino)", vec2string(mutant))
	            end
	        end
	    end
	end

	pushfirst!(vec_fit, 0)
	writedlm(path_fit_out, vec_fit)
end

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

function write_single_muts_MSA(wt::Array{<:Integer, 1}, fitness_wt::Array{Union{Missing, Float64}},
			dir_out::AbstractString; wt_name::AbstractString = "WT")

	path_MSA_out = joinpath(dir_out, "single_muts_MSA_$(wt_name).fasta")
	path_fit_out = joinpath(dir_out, "single_muts_fitness_$(wt_name).fit")
	vec_fit = []
	k = 0
	FastaWriter(path_MSA_out, "w") do file
	    writeentry(file, "$k | $(wt_name)", vec2string(wt))
	    for i in [pos for (pos, amino) in enumerate(wt) if amino!= 21]
	        mutant = copy(wt)
	        for amino in 1:20
	            mutant[i] = amino
	            new_amino = num2letter(mutant[i])
	            old_amino = num2letter(wt[i])
	            if !ismissing(fitness_wt[(i-1)*20 + amino]) && (amino != wt[i])
	                k+=1
	                append!(vec_fit, fitness_wt[(i-1)*20 + amino])
	                writeentry(file, "$k | $(old_amino)$(i)$(new_amino)", vec2string(mutant))
	            end
	        end
	    end
	end

	pushfirst!(vec_fit, 0)
	writedlm(path_fit_out, vec_fit)
end

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

function write_single_muts_MSA(path_wt::AbstractString, fitness_wt::Array{Any, 2},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")

	wt = fasta2matrix(path_wt)
	dir_out == "" && (  dir_out = (*).(split(path_wt, "/")[1:end-1])  )
	path_MSA_out = joinpath(dir_out, "single_muts_MSA_$(wt_name).fasta")
	path_fit_out = joinpath(dir_out, "single_muts_fitness_$(wt_name).fit")
	vec_fit = []

    l = 0
	FastaWriter(path_MSA_out, "w") do file
	    writeentry(file, "0 | $(wt_name)", vec2string(wt))
        for k in 1:size(fitness_wt, 1)

	        mutant = copy(wt)
            pos = fitness_wt[k, 1]

            old_amino_lettr = num2letter(wt[pos])
            new_amino_lettr = fitness_wt[k, 2][1]
            new_amino_num = letter2num(new_amino_lettr)

            mutant[pos] = new_amino_num

            if (new_amino_num != wt[pos])
                l +=1
                append!(vec_fit, fitness_wt[k, 3])
                writeentry(file, "$l | $(old_amino_lettr)$(pos)$(new_amino_lettr)", vec2string(mutant))
            end
	    end
	end

	pushfirst!(vec_fit, 0)
	writedlm(path_fit_out, vec_fit)
end

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

function write_single_muts_MSA(wt::Array{<:Integer, 1}, fitness_wt::Array{Any, 2},
			dir_out::AbstractString = ""; wt_name::AbstractString = "WT")

	dir_out == "" && (  dir_out = (*).(split(path_wt, "/")[1:end-1])  )
	path_MSA_out = joinpath(dir_out, "single_muts_MSA_$(wt_name).fasta")
	path_fit_out = joinpath(dir_out, "single_muts_fitness_$(wt_name).fit")
	vec_fit = []

    l = 0
	FastaWriter(path_MSA_out, "w") do file
	    writeentry(file, "0 | $(wt_name)", vec2string(wt))
        for k in 1:size(fitness_wt, 1)

	        mutant = copy(wt)
            pos = fitness_wt[k, 1]

            old_amino_lettr = num2letter(wt[pos])
            new_amino_lettr = fitness_wt[k, 2][1]
            new_amino_num = letter2num(new_amino_lettr)

            mutant[pos] = new_amino_num

            if (new_amino_num != wt[pos])
                l +=1
                append!(vec_fit, fitness_wt[k, 3])
                writeentry(file, "$l | $(old_amino_lettr)$(pos)$(new_amino_lettr)", vec2string(mutant))
            end
	    end
	end

	pushfirst!(vec_fit, 0)
	writedlm(path_fit_out, vec_fit)
end

##############################################################
"""
    function matrix2fasta(path::AbstractString, MSA, desc_name::AbstractString = "")

    (path, MSA) --> write MSA to path

    USE:
    Given a vector or matrix of numbers representing an amino acid, write it in letter format to "path".

    INPUT:
    "path": path of output
    "MSA": matrix or vector of numbers (1 to 21 in the convention used here)
    "desc_name": name of the sequence (works only if NOT an MSA)

    OUTPUT:

"""

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

function matrix2fasta(path::AbstractString, MSA, desc_name::AbstractString = "", start_index = 1)

    M = size(MSA, 1)

    single_seq_flag = (length(size(MSA)) == 1)

    if single_seq_flag == true
        FastaWriter(path, "w") do file
            desc_name == "" && (desc_name = "1")
            writeentry(file, desc_name, vec2string(MSA) )
        end
    else
	seq_index = start_index
        FastaWriter(path, "w") do file
            for i in 1:M
                seq = MSA[i, :]
                writeentry(file, "$(seq_index) $(desc_name)", vec2string(seq) )
		seq_index +=1
            end
        end
    end
end

max_stars_repo_name: matteobisardi/KitMSA

max_issues_repo_name: matteobisardi/KitMSA

max_forks_repo_name: matteobisardi/KitMSA


text:
using Distributions, Random, CSV, DataFrames

# num_tracts^length(percentage_buckets) possible states
num_tracts = 10
percentage_buckets = collect(0:0.2:1)

Random.seed!(238)

# education = ["high_school", "college_grad"]
econ = ["low_income", "avg_income", "high_income"]
# race = ["white", "black", "native", "alaska_native", "asian", "hawaiian", "pacific_islander"] # official census buckets
demo_types = econ

struct ResponseInc
  self_resp_inc::Distribution
  followup_inc::Distribution
end

response_inc_demo_file = "response_inc_demo.csv"
response_inc_demo = Dict()
for row in CSV.File(read(joinpath(@__DIR__, response_inc_demo_file)))
  response_inc_demo[row.demo] = ResponseInc(
    Normal(row.self_resp_inc_mean, row.self_resp_inc_std),
    Normal(row.followup_inc_mean, row.followup_inc_std),
  )
end

function genStates(num_states, demo, file)
  df = DataFrame()
  for d in demo
    df[d] = zeros(num_states)
  end
  for row in eachrow(df)
    r = rand(1, length(demo))
    r = r / sum(r)
    row[:] = r
  end
  CSV.write(joinpath(@__DIR__, file), df)
end

tract_file = "tract_demo.csv"
genStates(num_tracts, demo, tract_file)

tract_demo = CSV.File(joinpath(@__DIR__, tract_file)) |> DataFrame

struct Tract
  id::Int
  response_rate::Float64
  response_rate_bucket::Float64
end

function T(s, a)
  sp = copy(s)
  map(update_tracts!(a), sp)
  return sp
end

# findnearest gets the index of the item that is closest to the given value.
findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]

function update_tracts!(a)
  return function(tract::Tract)
    demo = tract_demo[tract.id]
    total_increase = 0
    if tract.id == a
      map_func = get_increase(demo, "self_resp_inc")
      avg(map(map_func, demo_types))
    else
      map_func = get_increase(demo, "followup_inc")
      avg(map(map_func, demo_types))
    end
    sp.response_rate += total_increase
    sp.response_rate_bucket = findnearest(percentage_buckets, response_rate)# round response rate to bucket
  end
end

function get_increase(demo::DataFrameRow, inc_type)
  return function(demo_type::String)
    sym = Symbol(inc_type)
    return rand(getfield(response_inc_demo[demo_type], sym)) * demo[demo_type]
  end
end

max_stars_repo_name: Census-MDP/Census-MDP

max_issues_repo_name: Census-MDP/Census-MDP

max_forks_repo_name: Census-MDP/Census-MDP


text:
module CF

import ..d

#cross entropy cost function

logloss{F <: AbstractFloat}(x::F) = -log(x)
dlogloss{F <: AbstractFloat}(x::F) = -one(F)/x

d(::typeof(logloss)) = dlogloss

crossentropy{F <: AbstractFloat}(expected::Bool, result::F) = (expected) * logloss(result) + (!expected) * logloss(one(F)-result)
crossentropy{F <: AbstractFloat}(expected::Array{Bool}, result::Array{F}) = sum(crossentropy.(expected, result))

dcrossentropy{F <: AbstractFloat}(expected::Bool, result::F) = (expected) * d(logloss)(result) - (!expected) * d(logloss)(one(F) - result)
dcrossentropy{F <: AbstractFloat}(expected::Array{Bool}, result::Array{F}) = dcrossentropy.(expected, result)

d(::typeof(crossentropy)) = dcrossentropy

#mean square cost function

meansquare{F <: AbstractFloat}(result::F, expected::F) = sqr(expected - result)
dmeansquare{F <: AbstractFloat}(result::F, expected::F) = F(2.0) * (result + expected)

d(::typeof(meansquare)) = dmeansquare
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
Base.isapprox(a::Array,b::CuArray;o...)=isapprox(a,Array(b);o...)
Base.isapprox(a::CuArray,b::Array;o...)=isapprox(Array(a),b;o...)


text:
function part1()
	crabs = [parse(Int, i) for i in split(readline("inputD7.txt"), ",")]
	highest = maximum(crabs)
	minfuel = minimum(n -> sum(abs(n - c) for c in crabs), 1:highest)

	return minfuel
end

function part2()
	crabs = [parse(Int, i) for i in split(readline("inputD7.txt"), ",")]
	highest = maximum(crabs)
	minfuel = minimum((n -> sum(sum(1:abs(n - c)) for c in crabs)), 1:highest)

	return minfuel
end

function part1_broadcasted()
	# An alternative solution to part1 that uses broadcasting instead of for-loops
	crabs = parse.(Int, split(readline("inputD7.txt"), ","))
	highest = maximum(crabs)
	minfuel = minimum(n -> sum(abs.(n .- crabs)), 1:highest)

	return minfuel
end

println(part1())
println(part2())


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export logbackward!

"""
`logbackward!(logbeta,Nt,A,y)`

In-place version of [`logbackward`](@ref).
"""
function logbackward!(beta::AbstractMatrix{T},
                      Nt::Integer,
                      A::AbstractMatrix{T},
                      y::AbstractMatrix{T}) where {T<:AbstractFloat}
  Ns = size(A,1)
  fill!(beta, -T(Inf))
  for i = 1:Ns
    beta[Nt,i] = zero(T)
  end
  for t = Nt-1:-1:1
    for j = 1:Ns
      for k = 1:Ns
        ajk = A[j,k]
        if ajk != -T(Inf)
          beta[t,j] = logadd(beta[t,j], beta[t+1,k]+ajk+y[t+1,k])
        end
      end
    end
  end
  return beta
end

# backward with constrained path
function logbackward!(beta::Matrix{T},
                      Nt::Integer,
                      A::AbstractMatrix{T},
                      t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                      y::AbstractMatrix{T}) where {D<:Integer,T<:AbstractFloat}
  Ns = size(A,1)
  fill!(beta, -T(Inf))
  for i = 1:Ns
    beta[Nt,i] = zero(T)
  end

  for t = Nt-1:-1:1
    I,J = t2IJ[t]
    for k in eachindex(I)
      i,j = I[k],J[k]
      beta[t,i] = logadd(beta[t,i], beta[t+1,j]+A[i,j]+y[t+1,j])
    end
  end
  return beta
end

export logbackward

"""
`logbackward(Nt,A,y)`

Computes the backward log-probabilities (``\\hat{\\beta}``) of an Hidden Markov Model (HMM) of `Ns` states with log-transition matrix `A` (must be a `Ns` × `Ns` matrix), scale coefficients and observation log-likelihoods `y` (must be of size `Nt2` × `Ns` with `Nt2` ≥ `Nt`).
"""
function logbackward(Nt::Integer,A::AbstractMatrix{T},y::AbstractMatrix{T}) where {T <: AbstractFloat}
  Ns = size(A,1)
  beta = zeros(T,Nt,Ns)
  logbackward!(beta,Nt,A,y)
end

"""
`logbackward(Nt,A,t2tr,y)`

Computes the backward log-probabilities with constrained path indicated by the `t2tr` vector. See [`forward`](@ref) for a description of the structure and meaning of `t2tr`.
"""
function logbackward(Nt::Integer,
                     A::AbstractMatrix{T},
                     t2tr,
                     y::AbstractMatrix{T}) where {T <: AbstractFloat, D <: Integer}
  Ns = size(A,1)
  beta = zeros(T,Nt,Ns)
  logbackward!(beta,Nt,A,t2tr,y)
end


text:
using ADCME
using Plots
using PyCall
@pyimport scipy.interpolate as si
pyplot()

################# --- functions --- #################

# flux function non-convex
function f_w(sw, M=1)
    # u := saturation of water
    # M := ratio of viscosities
    return sw^2 / (sw^2 + ((1 - sw)^2 / M))
end

# gradient of non-convex flux
function f_w_grad(sw, M=1)
    return 2 * M * sw * (1-sw) / ((M+1) * sw^2 - 2 * sw + 1)^2
end

# analytical Buckley-Leverett solution to sw (sw := Saturation of water)
function sw(x, t, M=1)
    # x := spacial coord
    # t := time nodes
    # M := viscosity ratio μ_oil/μ_water
    s = range(0, 1, length=length(x))[end:-1:1]
    x_s = @. f_w_grad(s)

    # compute sasturation before shock
    dif = x_s[2] - x_s[1]
    for i in 2:length(x)
        if x_s[i+1] - x_s[i] < dif
            x_s = x_s[1:i]
            @goto stop
        else
            dif = x_s[i+1] - x_s[i]
        end
    end
    @label stop

    x_s_times = []
    [push!(x_s_times, x_s * k) for k in t ]

    sw_ = zeros(Nx*Nt)

    # Perform the interpolation
    for k in 1:length(t)
        bench = x_s_times[k][end]
        tmp = si.griddata(Array(x_s_times[k]), Array(s), Array(x))
        sw_[(k-1)*Nx+1:k*Nx] = @. ifelse(x <= bench, tmp, 0.0)
        convert(Array{Float64,1}, x_s_times[k])
    end

    return sw_
end

### --- numerical parameters --- ###

X = 1 # spacial domain from 0 to 1
T = 1 # time domain from 0 to 1

dx = 0.01
dt = 0.01
Nx = round(Int, X/dx + 1) # spacial nodes
Nt = round(Int, T/dt) # time nodes

# simple coordinates
x = range(0, X, step=dx) # x coords
t = range(dt, T, step=dt) # t steps

# 1D array of all possible coordinte pairs
# will also be all colocation points for PDE regularization
x_star = repeat(x, Nt)
t_star = collect(Iterators.flatten(repeat(t', Nx)))

################# --- anal solution --- #################

sw_anal = sw(x, t)
sw_anal_tf = constant(sw_anal)

################# --- neural net --- #################

# training input for neural network
x_star_tf = constant(x_star)
t_star_tf = constant(t_star)
#train_input = constant(hcat(train_XX, train_TT))

# layers
config = [20, 20, 20, 20, 20, 20, 20, 20, 1]
# xavier initialization
θ = Variable(fc_init([2; config]))

# neural network function
sw_nn = squeeze(fc(hcat(x_star_tf, t_star_tf), config, θ)) + 1

# derivatives for PDE
sw_t_nn = tf.gradients(sw_nn, t_star_tf)[1]
sw_x_nn = tf.gradients(sw_nn, x_star_tf)[1]
sw_xx_nn = tf.gradients(sw_x_nn, x_star_tf)[1]
f_w_nn = f_w(sw_nn)
f_w_x_nn = tf.gradients(f_w_nn, x_star_tf)[1]

### --- loss --- ###

# loss funciton: (IC/BC - analitical_solution) + (sw_t + f_w(sw) - ϵ * sw_xx)
# sw_t       := time derivative of water saturation
# f_w(sw)    := flux funciton dependend on water saturation
# ϵ         := 1/Pe
# sw_xx      := second spatial derivative of water saturation

ϵ = constant(2.5*10e-3)

# Loss with diffusion
loss = (sum((sw_nn[1:Nx] - sw_anal_tf[1:Nx])^2) / Nx +
        sum((sw_nn[Nx+1:Nx:end] - sw_anal_tf[Nx+1:Nx:end])^2)/ (Nt - 1) +
        sum((sw_t_nn + f_w_x_nn - ϵ * sw_xx_nn)^2) / (Nx * Nt))
# Loss without diffusion
#= loss = (sum((u_nn[1:Nx] - u_anal_train[1:Nx])^2) / Nx +
        sum((u_nn[Nx+1:Nx:end] - u_anal_train[Nx+1:Nx:end])^2)/ (Nt - 1) +
        sum((u_nn[Nx*2-1:Nx:end] - u_anal_train[Nx*2-1:Nx:end])^2)/ (Nt - 1) +
        sum((u_t_nn + f_w_x_nn)^2) / (Nx * Nt)) =#

################# --- training --- #################

### optimising

losses = []
# callback function for optimisation
i = 0
function cb()
    global i
    if i % 10 == 0
        println("Epoch: $i")
    end
    #push!(losses, loss)
    i += 1
end

sess = Session(); init(sess)

# Scipy optimizer
opt = ScipyOptimizerInterface(loss; method="L-BFGS-B", options=Dict("maxiter"=> 3000, "ftol"=>1e-12, "gtol"=>1e-12))
ScipyOptimizerMinimize(sess, opt; loss_callback=cb)

#BFGS!(sess, loss)

### --- predict --- ###

sw_test_nn = run(sess, squeeze(fc(hcat(x_star_tf, t_star_tf), config, θ))+1)

error = sum(@. sqrt((sw_test_nn - sw_anal).^2)) / sum(@. sqrt(sw_anal.^2))
println("\nFinal Loss:\t", run(sess, loss))
println("sw error:\t", error)

### --- Plots --- ###

println("\nPrinting images...")
close("all")

# times at which to plot
t_plot = [0.25, 0.5, 0.75]
sw_anal = reshape(sw_anal, (Nx, Nt))
sw_test_nn = reshape(sw_test_nn, (Nx, Nt))

for i in 1:length(t_plot)
    idx = round(Int, t_plot[i] / dt * Nt*dt)
    println("0.$idx")
    plot(x, sw_anal[:, idx], label="anal")
    plot!(x, sw_test_nn[:, idx], label="NN", line=:dash)
    name = "./results/0.$idx.png"
    savefig(name)
end

### --- Save NN --- ###

println("saving NN...")
# note if you run script multiple times -> number epochs adds up
name = "./model.mat"
save(sess, name)

sess.close()

max_stars_repo_path: 1. Buckley-Leverett/nn_prediciton.jl

max_issues_repo_path: 1. Buckley-Leverett/nn_prediciton.jl

max_forks_repo_path: 1. Buckley-Leverett/nn_prediciton.jl


text:
using YearMonths, Test, Dates

@testset "Constructors" begin
    ym = YearMonth(2018, 12)
    @test year(ym) == 2018
    @test month(ym) == 12
    @test yearmonth(ym) == (2018, 12)
    show(ym)
    println("")

    ym = YearMonth(2018, 10)
    ym2 = YearMonth(2018, 11)

    @test ym == YearMonth("201810")
    @test ym == YearMonth("2018-10")
    @test string(ym) == "2018-10"
    @test ym == YearMonth(201810)
    @test ym != YearMonth(2018, 11)
    @test Date(ym) == Date(2018, 10)
    @test Date(ym) == Date(2018, 10, 1)

    @test hash(ym) == hash(YearMonth(2018, 10))
    @test hash(ym) != hash(ym2)

    @test string(YearMonth(0xff, 1)) == "255-01"
    @test string(YearMonth("20000-12")) == "20000-12"
    @test string(YearMonth("2000012")) == "20000-12"
    @test string(YearMonth("2-01")) == "2-01"
    @test string(YearMonth("201")) == "2-01"
    @test_throws ErrorException YearMonth("10")
end

@testset "Math" begin
    ym = YearMonth(2018, 12)
    @test ym - Month(1) == YearMonth(2018, 11)
    @test ym - Year(1) == YearMonth(2017, 12)
    @test ym + Month(1) == YearMonth(2019, 1)
    @test Month(1) + ym == YearMonth(2019, 1)
    @test ym + Year(1) == YearMonth(2019, 12)
    @test Year(1) + ym == YearMonth(2019, 12)
    @test ym + Month(12) == YearMonth(2019, 12)

    ym1 = YearMonth(2018, 12)
    ym2 = YearMonth(2019, 1)
    @test ym1 <= YearMonth(2018, 12)
    @test ym1 < ym2
    @test ym2 > ym1
    @test isequal(ym1, YearMonth(2018, 12))
    @test isequal(YearMonth(2018, 12), ym1)
    @test !isequal(ym1, ym2)
    @test isless(ym1, ym2)
    @test !isless(ym2, ym1)
end

@testset "Conversions" begin
    dt = Date(2018, 12, 10)
    ym1 = YearMonth(2018, 12)
    @test YearMonth(dt) == ym1
    @test_throws ErrorException YearMonth("2018-122")

    ym = YearMonth(2018, 5)
    @test firstdayofmonth(ym) == Date(2018, 5, 1)
    @test lastdayofmonth(ym) == Date(2018, 5, 31)
end

@testset "Tutorial" begin
    include("tutorial.jl")
end


text:
using LinearAlgebra
using StaticArrays
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}
const VecFrame = SVector{18, Float64}

using PyCall
using Conda
using Pkg

ENV["PYTHON"] = "D:\\Anaconda\\python.exe"
Pkg.build("PyCall")
println(PyCall.conda)
println(PyCall.libpython)
pyimport_conda("torch", "pytorch", "pytorch")
pyimport_conda("torchvision", "torchvision", "https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/")

py"""
import pandas as pd
import numpy as np
import keras
import matplotlib.pyplot as plt
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation, Flatten, LSTM, TimeDistributed, RepeatVector
from keras.models import load_model
import random
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import train_test_split
from keras.callbacks import ModelCheckpoint
from keras.models import load_model
import time

def get_frames(db,Vehicle_IDs,Frame_ID):#Vehicle_IDs in order, for better performance
    pred_n=10
    cnt=0
    vid=db['Vehicle_ID'].values
    f=db['Frame_ID'].values

    lx=db['Local_X'].values
    ly=db['Local_Y'].values
    v=db['v_Vel'].values
    a=db['v_Acc'].values
    l=db['Lane_ID'].values
    t=db['theta'].values
    frames_for_all_vehicles=[]
    case_list=[]
    distance=[] #if second model, we record the dx between line and current point to distinguish online case

    #Requirement: must know its previous 20 frames
    for i in range(len(db)):
        # print(vid[i])
        if(vid[i]==Vehicle_IDs[cnt] and Frame_ID==f[i]):
            record=[]

            #which model to predict, check previous 20 frames(not includes itself)
            flag=False
            for j in range(20):
                if(l[i-j]!=l[i-j-1]):
                  flag=True
                  distance.append(lx[i]-lx[i-j-1])
                  break
            if(flag==False):#B->before lane change: first model, no need to record distance dx
                case_list.append(0)
                distance.append(None)
            if(flag==True):#A->after lane change: second model
                case_list.append(1)

            #necessary information
            for j in range(pred_n):
                qid=i-pred_n+1+j
                record.append([lx[qid],ly[qid],v[qid],a[qid],l[qid],t[qid]])

            frames_for_all_vehicles.append(record)
            cnt+=1
            if(cnt==len(Vehicle_IDs)):
                break
    # print(len(frames_for_all_vehicles))
    return np.array(frames_for_all_vehicles),np.array(case_list),np.array(distance)

def predict_action(Vehicle_IDs,Frame_ID,single_model=False,noise_threshold=0.5,pred_threshold=0.8):
    pred_n=10

    #test case
    time_stamp="0820am-0835am"
    db=pd.read_csv("pre_data/"+time_stamp+"/add_label.csv")
    # print(len(db))
    single_model=False
    model_dir='./model/'
    model=None
    model0=None#0 means first model before lane change
    model1=None#1 means second model after lane change
    if(single_model):
        model_cate="_np"
        model_name="best_weights"+model_cate+".h5"
        model=load_model(model_dir+model_name)
    else:
        model_cate0="_B"
        model_cate1="_An"
        model_name0="best_weights_npnew"+model_cate0+".h5"
        model_name1="best_weights_npnew"+model_cate1+".h5"
        model0=load_model(model_dir+model_name0)
        model1=load_model(model_dir+model_name1)
    print("Model loaded")
    st1=time.time()
    # print(Vehicle_IDs)
    frames_list,case_list,distance=get_frames(db,Vehicle_IDs,Frame_ID)
    st2=time.time()
    if(single_model==False):
        # print(frames_list.shape)
        result0=model0.predict(frames_list)
        result1=model1.predict(frames_list)
        pred_category0=np.argmax(result0,axis=1)
        pred_category1=np.argmax(result1,axis=1)
        pred_category=(case_list==0).astype(int)*pred_category0+(case_list==1).astype(int)*pred_category1

        second=np.array([-1]*len(pred_category))#find second possible direction, if possibility of first possible direction >0.8, second = -1
        sorted_result0=np.sort(result0,axis=1)
        sorted_result1=np.sort(result1,axis=1)
        for i in range(len(Vehicle_IDs)):
            # print(sorted_result1)
            if(case_list[i]==0 and sorted_result0[i][-1]<pred_threshold):
                for j in range(3):
                    if(result0[i][j]==sorted_result0[i][-2]):
                        second[i]=j
                        break
            if(case_list[i]==1 and sorted_result1[i][-1]<pred_threshold):
                for j in range(3):
                    if(result1[i][j]==sorted_result1[i][-2]):
                        second[i]=j
                        break

    else:
        result=model.predict(frames_list)
        pred_category=np.argmax(result,axis=1)

    decision=[]
    for i in range(len(pred_category)):
        if(case_list[i]==1):#After lane change model

            #right lane change occurs in previous 20 frames
            if(distance[i]>0):
                if(pred_category[i]==2):
                    decision.append(pred_category[i])#turn right case
                elif(distance<=noise_threshold):
                    decision.append(0)#on line case
                else:
                    decision.append(pred_category[i])#keep line case or turn left case

            #left lane change occurs in previous 20 frames
            if(distance[i]<=0):
                if(pred_category[i]==1):
                    decision.append(pred_category[i])#turn left case
                elif(distance>=-noise_threshold):
                    decision.append(0)#on line case
                else:
                    decision.append(pred_category[i])#keep line case or turn right case

        else:#Before lane change model
            decision.append(pred_category[i])
    et=time.time()
    # print("et-st1:",et-st1)
    # print("et-st2:",et-st2)
    # print(decision,second,pred_category)
    # print(case_list,distance,result1,result0)
    return decision,second
"""

function get_Vehicle_IDs(frame::Vector{VecFrame})
    res = Int64[]
    for v in frame
        push!(res, v[1])
    end
    return res
end

function predict_action(prev_frame::Vector{VecFrame}, frame::Vector{VecFrame}, Frame_ID::Int64)
    Vehicle_IDs = get_Vehicle_IDs(frame)
    # println(Vehicle_IDs)
    decisions, seconds = py"predict_action"(Vehicle_IDs, Frame_ID)
    threshold_keep = [-0.075, 0.075]
    threshold_left = [-0.2, 0.0]
    threshold_right = [0.0, 0.2]
    keep_cal = 0.0
    left_cal = -0.1
    right_cal = 0.075
    frame1 = copy(frame)
    frame2 = copy(frame)
    frame3 = copy(frame)
    frame4 = copy(frame)
    frame5 = copy(frame)
    frame6 = copy(frame)
    frame7 = copy(frame)
    frame8 = copy(frame)
    frame9 = copy(frame)
    for i = 1:length(decisions)
        if seconds[i] != -1
            println("lalala")
        end
        v_id = Vehicle_IDs[i]
        decision = decisions[i]
        theta = 0.0
        for v in prev_frame
            if v[1] == v_id
                theta = atan((frame[i][5] - v[5])/(frame[i][6] - v[6]))
                break
            end
        end
        if decision == 0
            if !(threshold_keep[1] <= theta <= threshold_keep[2])
                theta = (theta + keep_cal)/2
            end
        elseif decision == 1 # left
            if !(threshold_left[1] <= theta <= threshold_left[2])
                theta = (theta + left_cal)/2
            end
        elseif decision == 2 # right
            if !(threshold_right[1] <= theta <= threshold_right[2])
                theta = (theta + right_cal)/2
            end
        end
        vx = frame[i][12] * sin(theta)
        vy = frame[i][12] * cos(theta)
        frame1[i] = setindex(frame1[i], frame[i][5] + 0.1 * vx, 5)
        frame2[i] = setindex(frame2[i], frame[i][5] + 0.2 * vx, 5)
        frame3[i] = setindex(frame3[i], frame[i][5] + 0.3 * vx, 5)
        frame4[i] = setindex(frame4[i], frame[i][5] + 0.4 * vx, 5)
        frame5[i] = setindex(frame5[i], frame[i][5] + 0.5 * vx, 5)
        frame6[i] = setindex(frame6[i], frame[i][5] + 0.6 * vx, 5)
        frame7[i] = setindex(frame7[i], frame[i][5] + 0.7 * vx, 5)
        frame8[i] = setindex(frame8[i], frame[i][5] + 0.8 * vx, 5)
        frame9[i] = setindex(frame9[i], frame[i][5] + 0.9 * vx, 5)
        frame1[i] = setindex(frame1[i], frame[i][6] + 0.1 * vy, 6)
        frame2[i] = setindex(frame2[i], frame[i][6] + 0.2 * vy, 6)
        frame3[i] = setindex(frame3[i], frame[i][6] + 0.3 * vy, 6)
        frame4[i] = setindex(frame4[i], frame[i][6] + 0.4 * vy, 6)
        frame5[i] = setindex(frame5[i], frame[i][6] + 0.5 * vy, 6)
        frame6[i] = setindex(frame6[i], frame[i][6] + 0.6 * vy, 6)
        frame7[i] = setindex(frame7[i], frame[i][6] + 0.7 * vy, 6)
        frame8[i] = setindex(frame8[i], frame[i][6] + 0.8 * vy, 6)
        frame9[i] = setindex(frame9[i], frame[i][6] + 0.9 * vy, 6)
    end
    return frame1, frame2, frame3, frame4, frame5, frame6, frame7, frame8, frame9
end

# Vehicle_IDs=[1,2,3]
# Frame_ID=300
# predict_action(Vehicle_IDs,Frame_ID)


text:
using ObjectDetectors
using GeometryBasics

@testset "Intersection over Union" begin
    core_args = ("class", 0.7, 256, 256, MACHINE, "detector")

    box1 = BoundingBoxDetection(Point2(0.0, 0.0), Point2(3.0, 4.0), core_args...)
    @test iou(box1, box1) == 1

    box2 = BoundingBoxDetection(Point2(1.0, 2.0), Point2(3.0, 4.0), core_args...)
    @test iou(box1, box2) == 1 / 3

    box3 = BoundingBoxDetection(Point2(1.0, 2.0), Point2(4.0, 6.0), core_args...)
    @test iou(box1, box3) == 1 / 5

    box4 = BoundingBoxDetection(Point2(5.0, 5.0), Point2(6.0, 6.0), core_args...)
    @test iou(box1, box4) == 0
end


text:
function secret_handshake(code::Integer)
    # NOTE: Handle wrong codes
    # code < 0 && error("Invalid secret handshake code")
    # NOTE: Handle common actions
    masks = [(1, "wink"), (2, "double blink"),
             (4, "close your eyes"), (8, "jump")]
    actions = map(x->code&x[1]==x[1]&&x[2], masks)
    filter!(x->x!=false, actions)
    # NOTE: Handle reversing bitmask
    code&16==16 && reverse!(actions)
    actions
end


text:
# Objective function & sense.
#
# SCIP only supports affine objectives. For quadratic or nonlinear objectives,
# the solver will depend on bridges with auxiliary variables. Single variable
# objectives are also accepted, but the type is not correctly remembered.

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SAF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SVF}) = true

function MOI.set(o::Optimizer, ::MOI.ObjectiveFunction{SAF}, obj::SAF)
    allow_modification(o)

    # reset objective coefficient of all variables first
    for v in values(o.mscip.vars)
        @SC SCIPchgVarObj(o, v[], 0.0)
    end

    # set new objective coefficients, summing coefficients
    for t in obj.terms
        v = var(o, t.variable_index)
        oldcoef = SCIPvarGetObj(v)
        newcoef = oldcoef + t.coefficient
        @SC SCIPchgVarObj(o, v, newcoef)
    end

    @SC SCIPaddOrigObjoffset(o, obj.constant - SCIPgetOrigObjoffset(o))

    return nothing
end

# Note that SCIP always uses a scalar affine function internally!
function MOI.set(o::Optimizer, ::MOI.ObjectiveFunction{SVF}, obj::SVF)
    aff_obj = SAF([AFF_TERM(1.0, obj.variable)], 0.0)
    return MOI.set(o, MOI.ObjectiveFunction{SAF}(), aff_obj)
end

function MOI.get(o::Optimizer, ::MOI.ObjectiveFunction{SAF})
    terms = AFF_TERM[]
    for vr = keys(o.mscip.vars)
        vi = VI(vr.val)
        coef = SCIPvarGetObj(var(o, vi))
        coef == 0.0 || push!(terms, AFF_TERM(coef, vi))
    end
    constant = SCIPgetOrigObjoffset(o)
    return SAF(terms, constant)
end

# Note that SCIP always uses a scalar affine function internally!
function MOI.get(o::Optimizer, ::MOI.ObjectiveFunction{SVF})
    aff_obj = MOI.get(o, MOI.ObjectiveFunction{SAF}())
    if (length(aff_obj.terms) != 1
        || aff_obj.terms[1].coefficient != 1.0
        || aff_obj.constant != 0.0)
        error("Objective is not single variable: $aff_obj !")
    end
    return SVF(aff_obj.terms[1].variable_index)
end

function MOI.set(o::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    allow_modification(o)
    if sense == MOI.MIN_SENSE
        @SC SCIPsetObjsense(o, SCIP_OBJSENSE_MINIMIZE)
    elseif sense == MOI.MAX_SENSE
        @SC SCIPsetObjsense(o, SCIP_OBJSENSE_MAXIMIZE)
    elseif sense == MOI.FEASIBLITY_SENSE
        @warn "FEASIBLITY_SENSE not supported by SCIP.jl" maxlog=1
    end
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ObjectiveSense)
    return SCIPgetObjsense(o) == SCIP_OBJSENSE_MAXIMIZE ? MOI.MAX_SENSE : MOI.MIN_SENSE
end

function MOI.modify(o::Optimizer, ::MOI.ObjectiveFunction{SAF},
                    change::MOI.ScalarCoefficientChange{Float64})
    allow_modification(o)
    @SC SCIPchgVarObj(o, var(o, change.variable), change.new_coefficient)
    return nothing
end

MOI.get(o::Optimizer, ::MOI.ObjectiveFunctionType) = SAF


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 2, Zadanie 5.2.
#

fl1 = Float32
fl2 = Float64

# starting points
p = fl1(0.01)
pp = fl2(0.01)

r = 3

# calculate the next sequence member
function calc(x, fl)
    return fl(x + fl(r*fl(x*fl(1 - x))))
end

# print out the consecutive members of the sequences
for n in 1:40

    global p = calc(p, fl1)
    global pp = calc(pp, fl2)
    println(n, " & ", p, " & ", pp, "\\\\")

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# AUTO GENERATED FILE - DO NOT EDIT

export textinput

"""
    textinput(;kwargs...)

A TextInput component.
Custom input with label and description. For more information, see: https://mantine.dev/core/text-input/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `description` (String; optional): Input description, displayed after label
- `disabled` (Bool; optional): The component can show it is currently unable to be interacted with
- `error` (String; optional): Displays error message after input
- `label` (String; optional): Input label, displayed before input
- `multiline` (Bool; optional): Will input have multiple lines?
- `placeholder` (String; optional): Placeholder, displayed when date is not selected
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Input border-radius from theme or number to set border-radius in px
- `required` (Bool; optional): Adds red asterisk on the right side of label
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Input size
- `style` (Dict; optional): Inline style override
- `type` (a value equal to: "number", "text", "password"; optional): The type of control to render.
- `value` (a list of or a singular dash component, string or number; optional): Input value
- `variant` (a value equal to: "default", "filled", "unstyled", "headless"; optional): Defines input appearance, defaults to default in light color scheme and filled in dark
"""
function textinput(; kwargs...)
        available_props = Symbol[:id, :className, :description, :disabled, :error, :label, :multiline, :placeholder, :radius, :required, :size, :style, :type, :value, :variant]
        wild_props = Symbol[]
        return Component("textinput", "TextInput", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
# ---
# title: 973. K Closest Points to Origin
# id: problem973
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Divide and Conquer, Heap, Sort
# link: <https://leetcode.com/problems/k-closest-points-to-origin/description/>
# hidden: true
# ---
#
# We have a list of `points` on the plane.  Find the `K` closest points to the
# origin `(0, 0)`.
#
# (Here, the distance between two points on a plane is the Euclidean distance.)
#
# You may return the answer in any order.  The answer is guaranteed to be unique
# (except for the order that it is in.)
#
#
#
# **Example 1:**
#
#
#
#     Input: points = [[1,3],[-2,2]], K = 1
#     Output: [[-2,2]]
#     Explanation:
#     The distance between (1, 3) and the origin is sqrt(10).
#     The distance between (-2, 2) and the origin is sqrt(8).
#     Since sqrt(8) < sqrt(10), (-2, 2) is closer to the origin.
#     We only want the closest K = 1 points from the origin, so the answer is just [[-2,2]].
#
#
# **Example 2:**
#
#
#
#     Input: points = [[3,3],[5,-1],[-2,4]], K = 2
#     Output: [[3,3],[-2,4]]
#     (The answer [[-2,4],[3,3]] would also be accepted.)
#
#
#
#
# **Note:**
#
#   1. `1 <= K <= points.length <= 10000`
#   2. `-10000 < points[i][0] < 10000`
#   3. `-10000 < points[i][1] < 10000`
#
#
## @lc code=start
using LeetCode

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


text:
# Definitions of functions to compile
twox(x) = 2x

const aa = [4, 5]
arrayfun(x) = x + aa[1] + aa[2]

jsin(x) = sin(x)

function arridx(i)
    a = collect(1.0:0.1:10.0)
    @inbounds i > 3 ? a[1] : a[5]
end

fsimple() = [0:.001:2;][end]

include("ode.jl")
fode() = ode23s((t,y)->2.0t^2, 0.0, [0:.001:2;], initstep = 1e-4)[2][end]

# Functions to compile and arguments to pass
funcalls = [
    (twox, Tuple{Int}, 4),
    (arrayfun, Tuple{Int}, 4),
    (jsin, Tuple{Float64}, 0.5),
    (arridx, Tuple{Int}, 4),
    (fsimple, Tuple{}, ()),
    (fode, Tuple{}, ()),         # Broken on Julia v1.2.0; works on Julia v1.3.0-rc3
]

StaticCompiler.exegen(funcalls)

using Formatting
@testset "exegen" begin
    cd("standalone") do
        for (func, tt, val) in funcalls
            fname = nameof(func)
            rettype = Base.return_types(func, tt)[1]
            fmt = StaticCompiler.Cformatmap[rettype]
            @test Formatting.sprintf1(fmt, func(val...)) == read(`./$fname`, String)
        end
    end
end


text:
## Lexer types
abstract type AbstractLex end

struct Keyword <: AbstractLex
    text::String
end

struct GText <: AbstractLex
    text::String
end

struct Slashed <: AbstractLex
    text::String
end

struct GNumber{T<:Union{Int,Float64}} <: AbstractLex
    val::T
end

struct GArray{N} <: AbstractLex
    name::String
    indices::NTuple{N,Any}
end
GArray(name::AbstractString, indices::NTuple{N,Any}) where N =
    GArray{N}(name, indices)

struct Table{N} <: AbstractLex
    name::GArray{N}
    body::String
end

struct Dots <: AbstractLex
    text::String
end

# perhaps Operator, Relation?

struct GCall <: AbstractLex
    name::String
    args::Vector{Any}
end

struct Parens <: AbstractLex
    args::Vector{Any}
end

struct StatementEnd <: AbstractLex end

# For comparisons and hashing, we deliberately exclude the indices of GArray,
# since dictionary lookup wants to return the variable independently of the
# specific indices being accessed.
Base.:(==)(::AbstractLex, ::AbstractLex) = false
Base.:(==)(a::T, b::T) where T<:Union{Keyword,GText,Slashed,Dots} = a.text == b.text
Base.:(==)(a::GNumber, b::GNumber) = a.val == b.val
Base.:(==)(a::GArray{N}, b::GArray{N}) where N = a.name == b.name # && a.indices == b.indices
Base.:(==)(a::Table{N}, b::Table{N}) where N = a.name == b.name && a.body == b.body
Base.:(==)(a::GCall, b::GCall) = a.name == b.name && a.args == b.args
Base.:(==)(a::Parens, b::Parens) = a.args == b.args
Base.:(==)(::StatementEnd, ::StatementEnd) = true

const hash_lex = Sys.WORD_SIZE == 64 ? 0xc4e3eb2da1eebf4a : 0x834afa53

function Base.hash(x::AbstractLex, h::UInt)
    h ⊻= hash_lex
    _hash(x, h)
end
_hash(x::Keyword, h::UInt) = hash(x.text, hash(1, h))
_hash(x::GText, h::UInt)   = hash(x.text, hash(2, h))
_hash(x::GArray, h::UInt)  = hash(x.name, hash(3, h)) # hash(x.indices, )
_hash(x::Table, h::UInt)   = hash(x.body, hash(x.name, hash(4, h)))
_hash(x::Slashed, h::UInt) = hash(x.text, hash(5, h))
_hash(x::GNumber, h::UInt) = hash(x.val, hash(6, h))
_hash(x::Dots, h::UInt)    = hash(x.text, hash(7, h))
_hash(x::GCall, h::UInt)   = hash(x.args, hash(x.name, hash(8, h)))
_hash(x::Parens, h::UInt)  = hash(x.args, hash(9, h))
_hash(x::StatementEnd, h::UInt) = hash(10, h)

Base.convert(::Type{String}, l::Union{Keyword,GText,Slashed}) = l.text
Base.string(l::Union{Keyword,GText,Slashed}) = l.text
Base.string(l::GNumber) = string(l.val)
function Base.string(part::GArray)
    str = part.name * '(' * part.indices[1]
    for j = 2:length(part.indices)
        str *= ',' * part.indices[j]
    end
    str *= ')'
    return str
end

Base.convert(::Type{Expr}, l::GNumber) = l.val
Base.convert(::Type{Expr}, l::GText) = Symbol(string(l))
Base.convert(::Type{Expr}, l::GArray) =
    Expr(:ref, Symbol(getname(l)), map(x->convert(Expr, x), l.indices)...)
Base.convert(::Type{Expr}, l::GCall) =
    Expr(:call, Symbol(getname(l)), map(x->convert(Expr, x), l.args)...)
function Base.convert(::Type{Expr}, l::Parens)
    @assert length(l.args) == 1
    ex = convert(Expr, l.args[1])
    return :( ($ex) )
end

getname(a::GArray) = a.name
getname(a::Table) = a.name
getname(a::GCall) = a.name
getname(a::GNumber) = nothing
getname(a::AbstractLex) = a.text

sexpr(l::AbstractLex) = sexpr(stdout, l)
sexprlf(l::AbstractLex) = (sexpr(l); println())
sexpr(io::IO, l::Union{Keyword,Slashed,Table,Dots}) = print(io, l)
sexpr(io::IO, l::GText)   = print(io, l.text)
sexpr(io::IO, l::GNumber) = print(io, l.val)
sexpr(io::IO, l::GArray)  = print_delim(io, l.name, l.indices, '[', ']')
sexpr(io::IO, l::GCall)   = print_delim(io, l.name, l.args, '(', ')')
sexpr(io::IO, l::Parens)  = print_delim(io, "", l.args, '(', ')')
sexpr(io::IO, l::StatementEnd) = print(io, ';')
function print_delim(io, name, args, opener, closer)
    print(io, name, opener)
    sexpr(io, args[1])
    for i = 2:length(args)
        print(io, ',')
        sexpr(io, args[i])
    end
    print(io, closer)
end

## Parser types

# Attributes of a Variable
struct VarInfo{Axs<:Tuple}
    typ::String  # free, positive, negative, integer, binary
    axs::Axs
    assignments::Vector{Pair{Any,Any}}   # assignments that set lo, up, ...

    function VarInfo{Axs}(typ::AbstractString, axs) where Axs
        typ ∈ vartypes || error(typ, " must be one of ", vartypes)
        new{Axs}(typ, axs, Pair{Any,Any}[])
    end
end
VarInfo(typ::AbstractString, axs::Tuple{Vararg{<:AbstractUnitRange}}) =
    VarInfo{typeof(axs)}(typ, axs)

struct ModelInfo
    equations::Vector{String}
    assignments::Dict{String,Any}

    ModelInfo(equations) = new(equations, Dict{String,Any}())
end


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

name = "papi"
version = v"0.2.0"

# Collection of sources required to build papi
sources = [
     "https://sourceforge.net/projects/perfmon2/files/libpfm4/libpfm-4.10.1.tar.gz" =>
     "c61c575378b5c17ccfc5806761e4038828610de76e2e34fac9f7fa73ba844b49",

     "https://bitbucket.org/icl/papi.git" =>
     "0fdac4fc7f95f0ac8039e431419a5133088911af",

     # Patch to apply
     joinpath(@__DIR__, "patch"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libpfm-4.10.1/
make -j${nprocs}
make PREFIX=$prefix install
mkdir $prefix/bin

cp examples/showevtinfo $prefix/bin/showevtinfo
cp examples/check_events $prefix/bin/check_events
cd ../papi/src
git apply ../../inherit.patch

./configure --prefix=$prefix --host=$host --enable-perfevent-rdpmc=no
make -j${nprocs}
make install
exit

"""
# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpapi", :libpapi),
    LibraryProduct(prefix, "libpfm", :libpfm),
    ExecutableProduct(prefix, "check_events", :check_events),
    ExecutableProduct(prefix, "showevtinfo", :showevtinfo),
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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


text:
@testset "Proportion distances" begin

    @testset "Biological Sequences" begin
        dnaA = dna"ATCGCCA-"
        dnaB = dna"ATCCCCTA"
        rnaA = rna"ATCGCCA-"
        rnaB = rna"ATCCCCTA"

        @test_approx_eq_eps pdistance(Mismatch, dnaA, dnaB) 0.375 1e-3
        @test_approx_eq_eps pdistance(Mismatch, rnaA, rnaB) 0.375 1e-3

        @test_approx_eq_eps pdistance(Mutated, dnaA, dnaB) 0.2857142857142857 1e-3
        @test_approx_eq_eps pdistance(Mutated, rnaA, rnaB) 0.2857142857142857 1e-3
    end

end


text:
immutable compare_distributions
    data::AbstractArray
    log_likehoods_ratio::AbstractArray
    sig_level::Float64
    xmin::Number
    V_test_stat::Float64
    V_p_val::Float64
    V_preff_distr::Int64
    C_b::Int64
    C_p_val::Float64
    C_preff_distr::Int64
end

function show(io::IO, x::compare_distributions)
    println(io, "Data:  $(x.data)")
    println(io, "Log likehood ratio: $(x.log_likehoods_ratio)")
    println(io, "Significance level: $(x.sig_level)")
    println(io, "Xmin: $(x.xmin)")
    println(io, "Vuong - test statictic: $(x.V_test_stat)")
    println(io, "Vuong - p-value: $(x.V_p_val)")
    println(io, "Vuong - preffered distribution: $(x.V_preff_distr)")
    println(io, "Clarke test - number of possitive values in log likehood ratio: $(x.C_b)")
    println(io, "Clarke test - p-value: $(x.C_p_val)")
    println(io, "Clarke test - preffered distribution: $(x.C_preff_distr)")
end

function compare_distributions(d1::con_powerlaw, d2::DataType, data::AbstractArray; sig_level = 0.05)
    if !(d2 <: ContinuousUnivariateDistribution)
        println("Both distributions should be continuos.")
        return Union{}
    end
    xmin = d1.θ
    data = sort(data)
    data = data[findfirst(x -> x >= xmin, data): end]
    d2 = fit(d2, data)

    _compare_distributions(d1,d2,data,xmin,sig_level)

end

function compare_distributions(d1::dis_powerlaw, d2::DataType, data::AbstractArray; sig_level = 0.05)
    if !(d2 <: DiscreteUnivariateDistribution)
        println("Both distributions should be discrete.")
        return Union{}
    end
    xmin = d1.θ
    data = sort(data)
    data = data[findfirst(x -> x >= xmin, data): end]
    d2 = fit(d2, data)

    _compare_distributions(d1,d2,data,xmin,sig_level)
end

function compare_distributions(d1::DataType, d2::DataType, data::AbstractArray, xmin::Number = 0; sig_level = 0.05)
    if !(((d1 <: ContinuousUnivariateDistribution) && (d2 <: ContinuousUnivariateDistribution)) || ((d1 <: DiscreteUnivariateDistribution) && (d2 <: DiscreteUnivariateDistribution)))
        println("Both distributions should be of same type (continuous or discrete).")
        return Union{}
    end

    if (xmin == 0)
        if (d1 <: DiscreteUnivariateDistribution) xmin = 1 end
        xmin = minimum(data)
    else
        data = sort(data)
        data = data[findfirst(x -> x >= xmin, data): end]
    end
    d1 = fit(d1, data)
    d2 = fit(d2, data)

    _compare_distributions(d1,d2,data,xmin,sig_level)
end

function compare_distributions(d1::DiscreteUnivariateDistribution, d2::DiscreteUnivariateDistribution, data::AbstractArray, xmin::Number = 1; sig_level = 0.05)
   _compare_distributions(d1,d2,data,xmin,sig_level)
end

function compare_distributions(d1::ContinuousUnivariateDistribution, d2::ContinuousUnivariateDistribution, data::AbstractArray, xmin::Number = 0; sig_level = 0.05)
   _compare_distributions(d1,d2,data,xmin,sig_level)
end

function _compare_distributions(d1::UnivariateDistribution, d2::UnivariateDistribution, data::AbstractArray, xmin::Number, sig_level)
    data = sort(data)
    xmin == 0 ? xmin = minimum(data) : data = data[findfirst(x -> x >= xmin, data): end]
    # Vuong's test
    log_likehoods_ratio = logpdf(d1,data) - logpdf(d2,data)
    n = length(log_likehoods_ratio)
    m = mean(log_likehoods_ratio)
    standard_deviation = std(log_likehoods_ratio)
    test_stat = sqrt(n) * m / standard_deviation
    v_p_val = cdf(Normal(), test_stat)
    v_preff_distr = 0
    if (test_stat > cdf(Normal(), 1 - sig_level/2))
        v_preff_distr = 1
    elseif (test_stat < -cdf(Normal(), 1 - sig_level/2))
        v_preff_distr = 2
    end

    # Clarke's test
    # b = number of possitive values in log_likehoods_ratio
    b = sum(x->x > 0 ? 1 : 0, log_likehoods_ratio)
    preff_distr = 0
    if (b >= n/2)
        pval = 2 * (1 - cdf(Binomial(n,0.5),b-1))
        if (pval <= sig_level) preff_distr = 1 end
    else
        pval = 2 * (cdf(Binomial(n,0.5),b))
        if (pval <= sig_level) preff_distr = 2 end
    end

    return compare_distributions(data,log_likehoods_ratio,sig_level,xmin,test_stat,v_p_val,v_preff_distr,b,pval,preff_distr)
end


text:
# This file was generated, do not modify it. # hide
linear = @load LinearRegressor pkg=MLJLinearModels
knn = @load KNNRegressor; knn.K = 4
tree_booster = @load EvoTreeRegressor; tree_booster.nrounds = 100
forest = @load RandomForestRegressor pkg=DecisionTree; forest.n_trees = 500
svm = @load SVMRegressor;


text:
#!/usr/bin/env julia

using Random

# mod must always be positive

function mod(a, n)
    if (t = a % n) < 0
        n + t
    else
        t
    end
end

#=
a^b (mod n) using the method of repeated squares.

The key here is that every integer can be written as a sum of powers of 2 (binary numbers)
and that includes the exponent. By repeated squaring we get a raised to a power of 2. Also
recall that a^b * a^c = a^(b + c), so rather than adding we multiply since we are dealing
with the exponent.
=#

function powerMod(a, d, n)
    v = 1 # Value
    p = a # Powers of a
    while d > 0
        if isodd(d) # 1 bit in the exponent
           v = mod(v * p, n)
        end
        p = mod(p^2, n) # Next power of two
        d >>>= 1
    end
    v
end

#=
Greatest common divisor, Euclidean version.
=#

function gcd(a, b)
    while b ≠ 0
        a, b = b, mod(a, b)
    end
    a
end

#=
Witness loop of the Miller-Rabin probabilistic primality test.
=#

function witness(a, n)
    u, t = n - 1, 0
    while iseven(u) # n = u * 2^t + 1
        t += 1   # Increase exponent
        u >>>= 1 # Decrease the multiplier
    end
    x = powerMod(a, u, n)
    for i in 1:t
        y = powerMod(x, 2, n)
        if y == 1 && x ≠ 1 && x ≠ n - 1
            return true
        end
        x = y
    end
    x ≠ 1
end

#=
Miller-Rabin probabilistic primality test: the chance of being wrong is ≈ 1/4 each pass through
the loop.
=#

function isPrime(n, k)
    if n < 2 || (n ≠ 2 && iseven(n)) # 0, 1, and even except for 2 are not prime.
        return false
    elseif n < 4 # 3 is prime
        return true
    end # We must test all others
    for j in 1:k
        a = rand(2:n - 2) # Choose a random witness
        if witness(a, n)
            return false
        end
    end
    true
end

#=
We need a random prime number in [low, high] and for now a 4^–100 chance of a composite is
good enough.
=#

function randomPrime(low, high)
    guess = 0 # Certainly not prime!
    while !isPrime(guess, 100)
        guess = rand(low:high) # Half will be even, the rest have Pr[prime] ≈ 1/log(N).
    end
    guess
end

#=
A safe prime is the one following a Sophie German prime. If prime(p) and prime(2p + 1) then
2p + 1 is a safe prime.
=#

function safePrime(low, high)
    p = randomPrime(low, high)
    while !isPrime(2 * p + 1,100)
        p = randomPrime(low, high)
    end
    return 2 * p + 1
end

#=
A generator must not be congruent to 1 for any of its powers that are
proper divisors of p – 1.  Since p is safe prime, there are only two:
2 and (p – 1) / 2. The number of such generators is 𝜑(p – 1).
=#

function generator(n, p)
    g = n
    q = (p - 1) ÷ 2
    while powerMod(g, 2, p) == 1 && powerMod(g, q, p) == 1
        g = g + 1
    end
    g
end

#=
Generate an efficient description of a cyclic group G of order p, with generator r.

Choose a random integer a ∊ {(p – 1)/2, ..., p − 1}

Compute b = r^a

The public key consists of the values (p, r, b)

The private key consists of the values (p, a)
=#

function keys(safe, k)
    f = safe ? safePrime : randomPrime
    p = f(big"2"^(k - 1), big"2"^k - 1)
    r = generator(big"2"^16 + 1, p)
    a = rand((p - 1) ÷ 2:p - 1)
    b = powerMod(r, a, p)
    ((p, a), (p, r, b))
end

#=
Choose a random k ∊ {1, ..., p – 2}

Compute 𝛾 = r^k (mod p)

Compute 𝛿 = m b^k (mod p)

The encrypted message is (𝛾, 𝛿)
=#

function encrypt(m, key)
    (p, r, b) = key
    k = rand(1:p - 2)
    𝛾 = powerMod(r, k, p)
    𝛿 = mod(m * powerMod(b, k, p), p)
    (𝛾, 𝛿)
end

#=
The decrypted message is 𝛿 𝛾^(p – 1 – a) (mod p)
=#

function decrypt(m, key)
    (p, a) = key
    (𝛾, 𝛿) = m
    mod((powerMod(𝛾, p - 1 - a, p) * 𝛿), p)
end

function encode(s)
    sum::BigInt = 0
    pow::BigInt = 1
    for c in s
        sum += pow * BigInt(c)
        pow *= 256
    end
    sum
end

#=
Transform a BigInt back into a string, subtracting off the 0xAA. We treat it as a base-256
integer and just pull off the digits.
=#

function decode(n)
    s = ""
    while n > 0
        s = s * Char(mod(n, 256))
        n ÷= 256
    end
    s
end

print("How many bits? ")

bits = parse(Int64, readline())

(prv, pub) = keys(false, bits)

println("pub = $pub")
println("prv = $prv")

print(">> ")
for m in eachline()
    c = encrypt(encode(m), pub); println("En[$m] = $c")
    t = decode(decrypt(c, prv)); println("De[$c] = $t")
    print(">> ")
end

max_stars_repo_name: darrelllong/El-Gamal-in-Julia

max_issues_repo_name: darrelllong/El-Gamal-in-Julia

max_forks_repo_name: darrelllong/El-Gamal-in-Julia


text:
module testMoments
using Test
using LinearAlgebra ; const la = LinearAlgebra

srcDir = dirname(dirname(@__FILE__))*"\\src\\"
include(srcDir*"Moments.jl")
import .Moments ; const mom = Moments

@testset "eᵢ" begin
    n = 7
    for k in 1:n
        e = mom.eᵢ(n,k)
        @test length(e) == n
        @test maximum(e) == 1
        @test minimum(e) == 0
        @test sum(e) == 1
    end
end

@testset "var_kron" begin
    A = repeat([[1,0,1]],2,2)
    B = repeat([[0,1,0]],2,2)
    A_tens_B =  mom.var_kron(A,B)
    @test A_tens_B == repeat([[1,1,1]],4,4)

    n = rand(2:4)
    m = rand(2:5)
    k = n*m
    A = [rand(0:2,5) for i in 1:n, j in 1:n]
    B = [rand(0:2,5) for i in 1:m, j in 1:m]
    A_tens_B =  mom.var_kron(A,B)
    @test size(A_tens_B) == (k,k)

    A =  mom.make_mon_expo(4,(1,1))
    B =  mom.make_mon_expo(4,(1,1);isle = false)
    C = mom.var_kron(A,B)
    for i in 1:4, j in 1:4
        rc(k) = (1+(k-1)*4):(4+(k-1)*4)
        @assert C[rc(i),rc(j)] == [A[i,j]] .+ B
    end

    # Need more
end

@testset "make_mom_expo_mat" begin
    mon_expo_mat =  mom.make_mon_expo(3,(1,1);isle = true)
    @test la.diag(mon_expo_mat) == [ [0, 0, 0],
                                  [2, 0, 0],
                                  [0, 2, 0],
                                  [0, 0, 2]]

    n = 4
    t = 2

    M_vec =  mom.make_mon_expo(n,t)
    M_mat =  mom.make_mon_expo(n,(t,t))

    ## Complex moments

    t = 2
    d = (2,2)
    n = sum(2 .* d)
    M_vec =  mom.make_mon_expo(d,t)
    @assert length(M_vec) == binomial(n+t,t)
    M_mat =  mom.make_mon_expo(d,(t,t))
    @assert M_vec == M_mat[1,:]
end

@testset "make_mom_expo_keys" begin
    n = 3
    t = (1,1)
    mek = mom.make_mon_expo(n,2)
    @test mek == [[0, 0, 0],
                            [1, 0, 0],
                            [0, 1, 0],
                            [0, 0, 1],
                            [2, 0, 0],
                            [1, 1, 0],
                            [1, 0, 1],
                            [0, 2, 0],
                            [0, 1, 1],
                            [0, 0, 2]]

    for i in 1:3
        n = rand(3:6,1)[1]
        t = rand(2:4,1)[1]
        mek = mom.make_mon_expo(n,t*2)
        @test size(unique(mek))[1]  == size(mek)[1]
    end
end

## Real

@testset "make_xxᵀ_tens_yyᵀ" begin
    d = (3,3)
    xxᵀtensyyᵀ = mom.make_xxᵀ_tens_yyᵀ(d)
    @test size(xxᵀtensyyᵀ) == (prod(d),prod(d))

    @test diag(xxᵀtensyyᵀ) ==  [[2, 0, 0, 2, 0, 0],
                                [2, 0, 0, 0, 2, 0],
                                [2, 0, 0, 0, 0, 2],
                                [0, 2, 0, 2, 0, 0],
                                [0, 2, 0, 0, 2, 0],
                                [0, 2, 0, 0, 0, 2],
                                [0, 0, 2, 2, 0, 0],
                                [0, 0, 2, 0, 2, 0],
                                [0, 0, 2, 0, 0, 2]]

    @test sum.(xxᵀtensyyᵀ) == repeat([4],prod(d),prod(d))

    for i in 1:5
        d = rand(2:5,2)
        xxᵀtensyyᵀ = mom.make_xxᵀ_tens_yyᵀ(d)
        halfsum(arr) = [sum(arr[1:d[1]]),sum(arr[d[1]+1:end])]
        @test halfsum.(xxᵀtensyyᵀ) == repeat([[2,2]],prod(d),prod(d))
    end
end

@testset "get_ℝ_block_diagᵀ" begin
    d = (3,3)
    t = (4,4)
    ℝ_block_diag = mom.get_ℝ_block_diag(d,t)

    @test sum([size(ℝ_block_diag[b])[1] for b ∈ keys(ℝ_block_diag)]) == binomial(sum(d)+t[1],t[1])
end

## Complex

@testset "make_xxᵀ⨂yyᵀ" begin
    d = rand(2:5,2)
    xx̄ᵀ_tens_yȳᵀ = mom.make_xxᵀ⨂yyᵀ(d)
end

@test mom.split_expo([1,1,1,2,2,2,3,3,3,4,4,4],3,3) == ([1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4])

@testset "get_xx̄yȳMM_blocks" begin
    d = (2,2) ; t = (2,2)
    Iᵗ = make_mon_expo(d,t[1])
    Iᵗ_temp = map(x ->  sum.(split_expo(x,d...)),Iᵗ)
    r_list  = map(v -> v[1]-v[2],Iᵗ_temp) ; s_list  = map(v -> v[3]-v[4],Iᵗ_temp)

    Iᵗ_d = Dict()
    for r ∈ -t[1]:t[1], s ∈ -t[1]:t[1]
        T = Iᵗ[(r_list .== r) .& (s_list .== s)]
        isempty(T) ? continue : Iᵗ_d[r,s] = T
    end
    nar = [i  for k in keys(Iᵗ_d) for i in Iᵗ_d[k]]
    @test sort(nar) == sort(Iᵗ)

    d₁ = rand(2:4); t₁ = rand(2:4)
    d = (d₁,d₁) ; t = (t₁,t₁)
    xx̄yȳMM_blocks = mom.get_xx̄yȳMM_blocks(d,t)
    for b in keys(xx̄yȳMM_blocks)
        REF = map(x-> mom.split_expo(x,d...),xx̄yȳMM_blocks[b])
        rp =  sum.(map(v -> v[1]-v[2],REF)) ; sp  = sum.(map(v -> v[3]-v[4],REF))

        @test length(unique(rp)) == 1
        @test length(unique(sp)) == 1
        @test (rp[1,1],sp[1,1]) == (0,0)
    end
end

@testset "get_γδ_dict" begin
    d₁ = rand(2:3); t₁ = rand(2:3)
    d = (d₁,d₁) ; t = (t₁,t₁)
    γδ_dict = Moments.get_γδ_dict(d,t)
    for α ∈ keys(γδ_dict)
        γ,δ = γδ_dict[α]
        @test length(unique(γ+δ)) == 1
        @test unique(γ+δ)[1] == α
    end
end

@testset "get_coef" begin

end

@testset "get_expo" begin
end

@testset "get_ℜℑααᶥββᶥ" begin
end

@testset "get_ℝℂ_coefexpo" begin
end

@testset "get_ℂ_block_diagᵀ" begin

end

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
@doc raw"""
    PolynomialKernel([a=1 [,c=1 [,d=3]]])

The polynomial kernel is a Mercer kernel given by:

```math
\kappa(\mathbf{x},\mathbf{y}) =
(a \mathbf{x}^\intercal \mathbf{y} + c)^d
\qquad \alpha > 0, \; c \geq 0, \; d \in \mathbb{Z}_{+}
```
"""
struct PolynomialKernel{T<:AbstractFloat,U<:Integer} <: MercerKernel{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    d::HyperParameter{U}
    function PolynomialKernel{T}(a::Real, c::Real, d::U) where {T<:AbstractFloat,U<:Integer}
        new{T,U}(HyperParameter(convert(T,a), interval(OpenBound(zero(T)), nothing)),
                 HyperParameter(convert(T,c), interval(ClosedBound(zero(T)), nothing)),
                 HyperParameter(d, interval(ClosedBound(one(U)), nothing)))
    end
end
function PolynomialKernel(a::T1=1.0, c::T2=one(T1), d::Integer=3) where {T1<:Real,T2<:Real}
    PolynomialKernel{floattype(T1,T2)}(a, c, d)
end

@inline eltypes(::Type{<:PolynomialKernel{T,U}}) where {T,U} = (T,U)
@inline thetafieldnames(κ::PolynomialKernel) = Symbol[:a, :c]

@inline polynomialkernel(z::T, a::T, c::T, d::U) where {T<:AbstractFloat,U<:Integer} = (a*z + c)^d

@inline basefunction(::PolynomialKernel) = ScalarProduct()
@inline function kappa(κ::PolynomialKernel{T}, z::T) where {T}
    polynomialkernel(z, getvalue(κ.a), getvalue(κ.c), getvalue(κ.d))
end