text:
module alphaPredictableComponents

using LinearAlgebra, Printf, Statistics
using DataFrames, CoupledFields, FFTW, StatsBase

export modelf, prca, pca, CVfn
export Rsq1, Rsq2, Rsq3
export Rsq, omega, iota
export LCCA
export CoDa

include("CoDa.jl")

"""
    parzen(tau::OrdinalRange{Int64}, M::Int64)

Parzen filter
"""
function parzen(tau::OrdinalRange{Int64}, M::Int64)
    # lags =-50:50; f=  parzen(lags, 50)
    # p=plot(x=lags, y=f, Geom.line); draw(PNG(6.6inch, 3.3inch), p)

    τ = abs.(tau)
    i = τ .>= M/2
    z =   (i).*(2(1 .- τ./M).^3) +  (.!i).*(1 .- 6(τ./M).^2 + 6(τ./M).^3)
    z *= 1.0/sum(z)
    return z
end

"""
    modelf(par::Vector, X::Matrix, modelfn::Function)

Find the components of ``X`` using `modelfn` and parameters `par`.
par[1] = ``α``, par[2] = ``τ``
"""
function modelf(par::Vector{Float64}, X::Matrix{Float64}, modelfn::Function)
    Z0 = CoDa.alphaTransform(X, par[1])
    Z = Z0 .- mean(Z0,dims=1)
    lagmax = floor(Int, size(Z,1)/8)
    return modelfn(par, Z, lagmax)
end

"""
    prca(par::Vector, X::Matrix, lagmax::Int)

Predictable Components Analysis.
"""
function prca(par::Vector{Float64}, X::Matrix{Float64}, lagmax::Int)
    local n::Int64
    local p::Int64
    local Cxx_fac_inv::Matrix{Float64}

    n,p = size(X)
    Cxx_fac_inv = cov(X)^(-0.5)
    Zx = X * Cxx_fac_inv
#    lagmax = floor(Int, n/4)
    filt = parzen(0:lagmax, lagmax)
    mc = StatsBase.crosscov(Zx, Zx, 0:lagmax, demean=false)
    M = zeros(Float64, p, p)
    for k in 1:lagmax
        m1 = mc[1+k,:,:]
        M += filt[1+k] * m1 * m1'
    end

    d,U = eigen(Symmetric(M))
    W = Cxx_fac_inv * reverse(U,dims=2)
    Rx = X * W
    A1 = zeros(Float64, p, 4)
    for j in 1:4
        A1[:,j] = LCCA(Rx[:,j:j], Zx, Int(par[2]))
    end
    A = Cxx_fac_inv * A1
    Ry = X * A
    return ModelObj(W, Rx, A, Ry, reverse(d), par, "αPrC")
end

"""
    pca(par::Vector, X::Matrix, lagmax::Int)

Principal Components Analysis
"""
function pca(par::Vector{Float64}, X::Matrix{Float64}, lagmax::Int)

    Csc = diagm(0=>vec(1 ./ std(X, dims=1)))
    Z = X * Csc
    F, d, V = svd(Z)
#    W = Csc * V
    Rx = F./std(F,dims=1)
    W = (X'X) \ (X'Rx[:,1:4])
    Ry = X*W
    return ModelObj(W, Rx, W, Ry, d.^2, par, "αPC")
end

"""
    CVfn(parm::Matrix, X::Matrix, modelfn::Function, cvmetric::Function)

Cross-validation
"""
function CVfn(parm::Matrix{Float64}, X::Matrix{Float64}, modelfn::Function, cvmetric::Function)

    # CV metric parameters
    # par[2] = τ lags

    # Model parameters
    # par[1] = alpha

    nrg = size(parm, 1)
    n, p = size(X)

    # Setup
    ta = 1:floor(Int, n/2)
    tb = setdiff(1:n,ta)

    lagmax = floor(Int, n/8)

    CVmetric = zeros(Float64, nrg)
    for i in 1:nrg
        percent = round(100*i/nrg)
        @printf " %03d%%" percent

        par = parm[i,:]
        τ = Int(par[2])
        Z = CoDa.alphaTransform(X, par[1])
        zm1 = mean(Z[ta,:], dims=1); # zm2 = mean(Z[tb,:],1)

        Z1train = Z[ta,:].-zm1
        Z1test = Z[tb,:].-zm1

        model1 = modelfn(par, Z1train, lagmax)
#        modelfn==pca && println("αPCA")
#        modelfn==pca && (model1 = Afn(model1, Z1train, τ))
        CVmetric[i] = cvmetric(model1, Z1test, τ)
        println([parm[i:i,:] CVmetric[i] ])
    end
    parm = [parm CVmetric]
    imax = argmax(parm[:,end])
    return parm[imax,:]

end

"""
    Rsq1(model::ModelObj, Xtest::Matrix, τ::Int)

Predictive R²
"""
function Rsq1(model::ModelObj, Xtest::Matrix{Float64}, τ::Int)
    Rcv = Xtest * model.W[:,1]
    Tcv = Xtest * model.A[:,1]
    return cor(Rcv[1:(end-τ),1], Tcv[(1+τ):end,1])^2
end

"""
    Rsq2(model::ModelObj, Xtest::Matrix, τ::Int)

Explanation of variance R²
"""
function Rsq2(model::ModelObj, Xtest::Matrix{Float64}, τ::Int)
    Rcv = Xtest * model.W[:,1]
    numer = (Xtest'Rcv) / (Rcv'Rcv) * (Rcv'Xtest)
    denom = Xtest'Xtest
    return tr(numer)/tr(denom)
end

"""
    Rsq3(model::ModelObj, Xtest::Matrix, τ::Int)

"""
function Rsq3(model::ModelObj, Xtest::Matrix{Float64}, τ::Int)
    Rcv = Xtest * model.W[:,1]
    numer  = sum(abs2, Rcv / (Rcv'Rcv) * (Rcv'Xtest))
    denom = sum(abs2, Xtest)
    return numer/denom
end

function Afn(model::ModelObj, X::Matrix{Float64}, τ::Int)

    p = size(X,2)
    Cxx_fac_inv = cov(X)^(-0.5)
    Zx = X * Cxx_fac_inv
    Rx = X * model.W
    A1 = zeros(Float64, p, 4)
    for j in 1:4
        A1[:,j] = LCCA(Rx[:,j:j], Zx, τ)
    end
    A = Cxx_fac_inv * A1
    Ry = X * A
    return ModelObj(model.W, Rx, A, Ry, model.evals, model.pars, "αPCA")
end

"""
    LCCA(X::Matrix, Y::Matrix, τ::Int)

Linear Canonical Correlation Analysis between ``X`` and ``Y`` at lead time ``τ``
"""
function LCCA(X::T, Y::T, τ::Int) where T<:Matrix{Float64}
    nr = size(X,1)-τ
    Xa = X[1:nr,:]
    Ya = Y[(1+τ):end,:]
    M = cor(Xa, Ya)
    U, d, V = svd(M)
#    sv = cor(Xa*U[:,1], Ya*V[:,1])
    return V[:,1]
end

"""
    biplot(model::ModelObj, jv::Vector{Int}, Y::Matrix, label::Vector{Symbol})

Biplot
"""
function biplot(model::ModelObj, jv::Vector{Int}, Y::Matrix{Float64}, label::Vector{Symbol})

    Z = CoDa.alpha(Y, model.pars[1])
    R = model.R[:,jv]
    B = cov(R, zscore(Z,1))

    j = string(model.method,jv[2])
    return DataFrame(label=string.(label), B1=B[1,:], B2=B[2,:], x0=0.0, y0=0.0, a=string("α=",model.pars[1]), j=j)
end

"""
    Rsq(R::Matrix, X0::Matrix)

Explanation of variance R²
"""
function Rsq(R::Matrix{Float64}, X0::Matrix{Float64})
    numer = mapslices(Rj->tr((X0'Rj)/(Rj'Rj)*(Rj'X0)), R, dims=1)
    denom = tr(X0'X0)
    return vec(numer)./denom
end

"""
    omega(X::Matrix)

Spectral entropy
"""
function omega(X::Matrix{Float64})
    Xfft = fft(X,1)[2:end,:]
    sdf = real(Xfft .* conj(Xfft))
    n = size(Xfft,1)
    sdf[sdf.<0.0] .= 0.0
    sdf ./= sum(sdf,dims=1)
    o =  sdf .* log.(n, sdf)
    o[isinf.(o)] .= 0.0
    ιfca = 1.0 .+ sum(o,dims=1)
    return vec(ιfca)
end

omega(X::Matrix{Float64}, X0::Matrix{Float64}) = omega(X)

"""
    iota(model::ModelObj, iotafn::Function, jv::Vector, X0::Matrix, pos::Vector)

Iota
"""
function iota(model::ModelObj, iotafn::Function, jv::Vector{Int}, X0::Matrix{Float64}, pos::Vector)
    R = model.R[:,jv]
    ι = round.(Int, 100 .*iotafn(R, X0))
    istr = iotafn==omega ? string.("Ω = ",ι, "%") : string.(ι, "%")
    j = model.method .*string.(jv)
    return DataFrame(j=j, label=istr, x=pos[1], y=pos[2])
end

end # module


text:
function make_slist(suppZ)
    # one way to generalize is allow for more general slists
    # for these simulations we are only using fully saturated
    slist = [((d,z) -> Int((d == d̄) * (z == z̄)))
             for d̄ in 0:1, z̄ in eachrow(suppZ)][:]
end

function compute_βₛ(dgp::DGP)
    Γₛ = compute_Γₛ([(dgp.mtrs[1].basis, dgp.mtrs[2].basis)], dgp)
    # Only one model here, so Γₛ is a 1 x 2 array of matrices
    βₛ = fill(NaN, size(Γₛ[1,1])[1])
    for s in 1:length(βₛ)
        βₛ[s] = sum(Γₛ[1,1][s,:,:] .* dgp.mtrs[1].θ +
                    Γₛ[1,2][s,:,:] .* dgp.mtrs[2].θ)
    end
    return βₛ
end
export compute_βₛ

# return a 2-dimensional array Γₛ[ℓ, d], where each component is itself
# a 3-dimensional array with elements [s, j, k]
function compute_Γₛ(bases::Array{Tuple{MTRBasis, MTRBasis}, 1}, dgp::DGP)
    [compute_Γₛ(basis[d + 1], d, dgp) for basis in bases, d in 0:1]
end

function compute_Γₛ(basis::MTRBasis, d::Integer, dgp::DGP)
    @assert d in [0,1]
    slist = make_slist(dgp.suppZ)
    Γₛ = zeros(length(slist), length(basis.a), length(basis.b))
    for (i,z) in enumerate(eachrow(dgp.suppZ))
        intlb = (1 - d) * dgp.pscore[i]
        intub = d * dgp.pscore[i] + (1 - d) * 1
        Γₛ += [(aj(z) * ibk(intlb, intub) * s(d, z))
                for s in slist, aj in basis.a, ibk in basis.ib] .* dgp.densZ[i]
    end
    return Γₛ
end

ext: jl


text:
module QuantumBehavior

# """
# QuantumBehavior Module
#
# Constructs bell scenario behaviors using the laws of quantum mechanics.
# """

using ..Behavior
using LinearAlgebra

# """
# bipartite_scenario(α_povms, β_povms, ψ)
#
#     Constructs the quantum behavior corresponding to the system described by the
#     provided POVM operators and wavefunction.
#
# Assumption:
#     The number of outputs is uniform across inputs.
#
# Input:
#     α/β_povms: Array, contains the number of POVMs as there are inputs. Each POVM
#                is an array of matrices containing the num-outputs matrices.
#     ψ: Vector, the total wavefunction specified for the bipartite system as a whole.
#
# Output:
#     behavior: Vector, outcome probabilities for all possible events in bell scenario.
#               by default, the non-signaling representation is used.
# """
function bipartite_scenario(α_povms, β_povms, ψ)

    α_expt = (size(α_povms)[1], size(α_povms[1])[1])
    β_expt = (size(β_povms)[1], size(β_povms[1])[1])

    gen_behavior = zeros((Behavior.dimension(α_expt,β_expt,"generalized"),1))
    gen_behavior[1] = 1

    xc = 0
    yc = 0
    ac = 0
    bc = 0
    for x_povm in α_povms
        xc += 1
        for y_povm in β_povms
            yc += 1
            for a in x_povm
                ac += 1
                for b in y_povm
                    bc += 1
                    p_abxy = ψ'*kron(a,b)*ψ

                    index = Behavior.index((ac,bc,xc,yc),α_expt,β_expt,"generalized")

                    gen_behavior[index] = p_abxy[1]
                end
                bc = 0
            end
            ac = 0
        end
        yc = 0
    end
    xc = 0

    behavior = nothing
    if (α_expt[2] < 2) & (β_expt[2] < 2)
        behavior = gen_behavior
    elseif (α_expt[2] < 2) | (β_expt[2] < 2)
        behavior = Behavior.gen_to_norm_proj(α_expt,β_expt)*gen_behavior
    else
        behavior = Behavior.gen_to_ns_proj(α_expt,β_expt)*gen_behavior
    end

    behavior
end

# """
# prepare_and_measure(ρ_set, povms):
#
#     Computes the quantum behavior for a prepare and measure scenario where Alice
#     encodes their input into a quantum state (ρ) and sends it to Bob who then
#     measures the quantum state with the specified POVM.
#
# Input:
#     ρ_set: Array, contains valid density matrices
#     povms: Array, contains a valid povm arrays. One povm for each of Bob's inputs
#
# Output:
#     behavior: Col Vector, a behavior represented in the normalized subspace.
# """
function prepare_and_measure(ρ_set, povms; rep="normalized")

    α_in = length(ρ_set)
    β_in = length(povms)
    num_out = length(povms[1])

    α_expt = (α_in,1)
    β_expt = (β_in,num_out)

    gen_behavior = zeros((Behavior.dimension(α_expt,β_expt,"generalized"),1))
    gen_behavior[1] = 1

    xc = 0
    yc = 0
    bc = 0
    for ρ in ρ_set
        xc += 1
        for povm in povms
            yc += 1
            for b in povm
                bc += 1
                index = Behavior.index((1,bc,xc,yc),α_expt,β_expt,"generalized")
                p_abxy = real(tr(ρ*b))
                gen_behavior[index] = p_abxy
            end
            bc = 0
        end
        yc = 0
    end

    behavior = gen_behavior
    if rep == "normalized"
        norm_proj = Behavior.gen_to_norm_proj(α_expt,β_expt)
        behavior = norm_proj*gen_behavior
    end

    behavior
end

end


text:
### Statistics

function countmap{T1,T2}(x::RLEVector{T1,T2})
  tally = Dict{T1,T2}()
  for i in 1:nrun(x)
    tally[x.runvalues[i]] = get(tally,x.runvalues[i],0) + x.runvalues[i]
  end
  return(tally)
end

function mode(x::RLEVector)
  tally = countmap(x)
  which_max = indmax(collect(values(tally)))
  max_val = collect(keys(tally))[ which_max ]
  return(max_val)
end


text:
function _make_spectogram_gui()

    popupmenu_spect = Menu()
    popupmenu_spect_freq = MenuItem("Frequency Range")
    popupmenu_spect_win = MenuItem("Window Size")
    popupmenu_spect_overlap = MenuItem("Window Overlap")
    push!(popupmenu_spect,popupmenu_spect_freq)
    push!(popupmenu_spect,popupmenu_spect_win)
    push!(popupmenu_spect,popupmenu_spect_overlap)

    popupmenu_spect_freq_select=Menu(popupmenu_spect_freq)
    if VERSION > v"0.7-"
        spect_f_handles=Array{MenuItemLeaf}(undef,0)
    else
        spect_f_handles=Array{MenuItemLeaf}(0)
    end
    spect_f_options=[300; 1000; 3000; 7500; 15000]

    push!(spect_f_handles,MenuItem(string(spect_f_options[1])))
    push!(popupmenu_spect_freq_select,spect_f_handles[1])

    for i=2:5
        push!(spect_f_handles,MenuItem(spect_f_handles[i-1],string(spect_f_options[i])))
        push!(popupmenu_spect_freq_select,spect_f_handles[i])
    end

    popupmenu_spect_win_select=Menu(popupmenu_spect_win)
    if VERSION > v"0.7-"
        spect_w_handles=Array{MenuItemLeaf}(undef,0)
    else
        spect_w_handles=Array{MenuItemLeaf}(0)
    end
    spect_w_options=[10; 50; 100]

    push!(spect_w_handles,MenuItem(string(spect_w_options[1])))
    push!(popupmenu_spect_win_select,spect_w_handles[1])

    for i=2:3
        push!(spect_w_handles,MenuItem(spect_w_handles[i-1],string(spect_w_options[i])))
        push!(popupmenu_spect_win_select,spect_w_handles[i])
    end

    Gtk.showall(popupmenu_spect)

    (spect_w_handles,spect_f_handles, popupmenu_spect)
end

function add_spect_callbacks(spect_f_handles,spect_w_handles,handles)

    for i=1:5
        signal_connect(spect_popup_freq_cb,spect_f_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    for i=1:3
        signal_connect(spect_popup_win_cb,spect_w_handles[i],"activate",Void,(),false,(handles,i-1))
    end

    nothing
end

function prepare_spectrogram(ctx,han)

    myheight=height(ctx)
    mywidth=width(ctx)

    set_source_rgb(ctx,1.0,1.0,1.0)

    move_to(ctx,50.0,myheight-35)
    show_text(ctx,"-1")

    move_to(ctx,50.0,myheight-35)
    show_text(ctx,"-1.0")

    move_to(ctx,(mywidth-50.0)/2+50,myheight-35)
    show_text(ctx,"-0.5")

    move_to(ctx,mywidth-20.0,myheight-35)
    show_text(ctx,"0.0")

    move_to(ctx,mywidth/2,myheight-20)
    show_text(ctx,"Time (s)")

    move_to(ctx,10.0,myheight-140)
    rotate(ctx,-pi/2)
    show_text(ctx,"Frequency")
    identity_matrix(ctx)

    move_to(ctx,35.0,myheight-50.0)
    rotate(ctx,-pi/2)
    show_text(ctx,"0.0")
    identity_matrix(ctx)

    move_to(ctx,35.0,myheight-50-125)
    rotate(ctx,-pi/2)
    show_text(ctx,string(round(Int64,han.spect.f_max*han.spect.f_div/2)))
    identity_matrix(ctx)

    nothing
end

function draw_spectrogram(rhd::RHD2000,han::Gui_Handles)

    ctx=Gtk.getgc(han.c)

    x=Gtk.cairo_surface(han.c)
    data = ccall((:cairo_image_surface_get_data,Cairo._jl_libcairo),Ptr{UInt32},(Ptr{Void},),x.ptr)

    c_h=round(Int64,height(ctx))
    c_w=round(Int64,width(ctx))-50

    for i=(size(rhd.v,1)+1):rhd.sr
        han.v_s[i-SAMPLES_PER_DATA_BLOCK] = han.v_s[i]
    end
    count=1

    if han.band_widgets.lfp_en[han.sc.spike]
        for i=(rhd.sr-SAMPLES_PER_DATA_BLOCK+1):rhd.sr
            han.v_s[i] = rhd.lfps[count,han.sc.spike]
            count+=1
        end
    else
        for i=(rhd.sr-SAMPLES_PER_DATA_BLOCK+1):rhd.sr
            han.v_s[i] = rhd.v[count,han.sc.spike]
            count+=1
        end
    end

    plot_spectrogram(han.v_s,rhd.sr,han.spect)

    in_w = han.spect.t_max
    in_h = han.spect.f_max

    scale_w = c_w / in_w
    scale_h = 250 / in_h

    mymin=minimum(han.spect.out)
    mymax=maximum(han.spect.out)
    myrange=mymax-mymin

    rgb_mat = zeros(UInt32,in_h,in_w)

    for h=1:in_h
        for w=1:in_w

            startround = (han.spect.out[h,w]-mymin)/myrange*255
            myinput::UInt8 = (startround >= 255) ? 255 : floor(UInt8,startround)+1

            myblue=jet_b[myinput]
            mygreen=jet_g[myinput]
            myred=jet_r[myinput]

            rgb_mat[h,w] = 0
            rgb_mat[h,w] |= myblue  # blue
            rgb_mat[h,w] |= (mygreen << 8) #green
            rgb_mat[h,w] |= (myred << 16) #red
        end
    end

    for h=1:250
        for w=1:c_w
            val = rgb_mat[ceil(Int,h/scale_h),ceil(Int,w/scale_w)]
            unsafe_store!(data,val,(c_h-h-50)*(c_w+50)+w+50)
        end
    end

    nothing
end

function plot_spectrogram(s,fs,spect)

    S = spectrogram(s,spect.win_width_s,spect.win_overlap_s; fs=fs, window=hanning)

    temp = power(S)

    for i=1:spect.f_max
        for j=1:spect.t_max
            spect.out[i,j]=log10(temp[i,j])
        end
    end

    nothing
end

function spect_popup_freq_cb(widget::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if event_id == 0
        f_max = 300
    elseif event_id == 1
        f_max = 1000
    elseif event_id == 2
        f_max = 3000
    elseif event_id == 3
        f_max = 7500
    else
        f_max = 15000
    end

    han.spect = Spectrogram(han.spect.fs; win_width_t = han.spect.win_width_t, win_overlap_t = han.spect.win_overlap_t, f_max = f_max)

    nothing
end

function spect_popup_win_cb(widget::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if event_id == 0
        mywin = .01
    elseif event_id == 1
        mywin = .05
    else
        mywin = .1
    end

    han.spect = Spectrogram(han.spect.fs; win_width_t = mywin, win_overlap_t = han.spect.win_overlap_t, f_max = han.spect.f_max*han.spect.f_div)

    nothing
end


text:
using DelimitedFiles
using LightGraphs

function contiguous_usa()
    es, _ = readdlm("networks/edgesUSA.csv", ',', Int, header=true)
    es .+= 1
    return  SimpleGraphFromIterator([Edge(es[i,1] => es[i,2]) for i in 1:size(es,1)])
end

function london_transport()
    es = readdlm("networks/LondonTransportEdgesSimplified.csv", ',', Int)
    es .+= 1
    return  SimpleGraphFromIterator([Edge(es[i,1] => es[i,2]) for i in 1:size(es,1)])
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
mutable struct LSDiagnostic{C}
    loc::UnitRange{Int}
    actions::Vector{Any}
    message::String
end

mutable struct Document
    _uri::String
    path::String
    _content::String
    _line_offsets::Union{Nothing,Vector{Int}}
    _open_in_editor::Bool
    _workspace_file::Bool
    cst::EXPR
    diagnostics::Vector{LSDiagnostic}
    _version::Int
    _runlinter::Bool
    server
    root::Union{Nothing,Document}
end

function Document(uri::AbstractString, text::AbstractString, workspace_file::Bool, server = nothing)
    path = uri2filepath(uri)
    cst = CSTParser.parse(text, true)
    doc = Document(uri, path, text, nothing, false, workspace_file, cst, [], 0, true, server, nothing)
    get_line_offsets(doc)
    cst.val = path
    doc.cst.ref = doc
    setroot(doc, doc)
    return doc
end

function get_text(doc::Document)
    return doc._content
end

function set_open_in_editor(doc::Document, value::Bool)
    doc._open_in_editor = value
end

function get_open_in_editor(doc::Document)
    return doc._open_in_editor
end

function is_workspace_file(doc::Document)
    return doc._workspace_file
end

"""
    get_offset(doc, line, char)

Returns the byte offset position corresponding to a line/character position.
This takes 0 based line/char inputs. Corresponding functions are available for
Position and Range arguments, the latter returning a UnitRange{Int}.
"""
function get_offset(doc::Document, line::Integer, character::Integer)
    line_offsets = get_line_offsets(doc)
    offset = line_offsets[line + 1]
    while character > 0
        offset = nextind(doc._content, offset)
        character -= 1
    end
    return offset
end
get_offset(doc, p::Position) = get_offset(doc, p.line, p.character)
get_offset(doc, r::Range) = get_offset(doc, r.start):get_offset(doc, r.stop)

"""
    get_line_offsets(doc::Document)

Updates the doc._line_offsets field, an n length Array each entry of which
gives the byte offset position of the start of each line. This always starts
with 0 for the first line (even if empty).
"""
function get_line_offsets(doc::Document, force = false)
    if force || doc._line_offsets === nothing
        doc._line_offsets = Int[0]
        text = doc._content
        ind = firstindex(text)
        while ind <= lastindex(text)
            c = text[ind]
            nl = c == '\n' || c == '\r'
            if c == '\r' && ind + 1 <= lastindex(text) && text[ind + 1] == '\n'
                ind += 1
            end
            nl && push!(doc._line_offsets, ind)
            ind = nextind(text, ind)
        end
    end
    return doc._line_offsets
end

function get_line_of(line_offsets::Vector{Int}, offset::Integer)
    nlines = length(line_offsets)
    if offset > last(line_offsets)
        line = nlines
    else
        line = 1
        while line < nlines
            if line_offsets[line] <= offset < line_offsets[line + 1]
                break
            end
            line += 1
        end
    end
    return line, line_offsets[line]
end

"""
    get_position_at(doc, offset)

Returns the 0-based line and character position within a document of a given
byte offset.
"""
function get_position_at(doc::Document, offset::Integer)
    offset > sizeof(doc._content) && error("offset[$offset] > sizeof(content)[$(sizeof(doc._content))]")
    line_offsets = get_line_offsets(doc)
    line, ind = get_line_of(doc._line_offsets, offset)
    char = 0
    while offset > ind
        ind = nextind(doc._content, ind)
        char += 1
    end
    return line - 1, char
end

"""
    Range(Doc, rng)
Converts a byte offset range to a LSP Range.
"""
function Range(doc::Document, rng::UnitRange)
    start_l, start_c = get_position_at(doc, first(rng))
    end_l, end_c = get_position_at(doc, last(rng))
    rng = Range(start_l, start_c, end_l, end_c)
end


text:
#-------------------------------------------------------------------------
function calculatekernelmatrix(X, Y, k, hyper)
#-------------------------------------------------------------------------

    K = zeros(length(X), length(Y))

    calculatekernelmatrix!(K, X, Y, k, hyper)

    return K

end

#-------------------------------------------------------------------------
function calculatekernelmatrix!(K, X, Y, k, hyper)
#-------------------------------------------------------------------------
    @assert(size(K, 1) == length(X))
    @assert(size(K, 2) == length(Y))

    # calculate kernel matrix
    for (i, xᵢ) in enumerate(X), (j, xⱼ) in enumerate(Y)
        @inbounds K[i,j] = k(xᵢ, xⱼ, hyper) + ((i==j) ? 1e-12 : 0.0)
    end
end

#-------------------------------------------------------------------------
function calculatekernelmatrix!(K, dist2, k, hyper)
#-------------------------------------------------------------------------
    @assert(size(K, 1) == size(dist2, 1))
    @assert(size(K, 2) == size(dist2, 2))

    copyto!(K, k(dist2, hyper))

    nothing
end

#-------------------------------------------------------------------------
function rbf(xᵢ, xⱼ, hyper)
#-------------------------------------------------------------------------
    makepos(hyper[1]) * exp(-0.5 * makepos(hyper[2]) * (xᵢ-xⱼ)^2)
end

#-------------------------------------------------------------------------
function rbf(dist2, hyper)
#-------------------------------------------------------------------------
    makepos(hyper[1]) .* exp.(-0.5 .* makepos(hyper[2]) .* dist2)
end

#-------------------------------------------------------------------------
function OU(xᵢ, xⱼ, hyper)
#-------------------------------------------------------------------------
    makepos(hyper[1]) * exp(-0.5 * makepos(hyper[2]) * abs(xᵢ-xⱼ))
end

#-------------------------------------------------------------------------
function OU(dist2, hyper)
#-------------------------------------------------------------------------
    makepos(hyper[1]) .* exp.(-0.5 .* makepos(hyper[2]) .* sqrt.(dist2))
end

#-------------------------------------------------------------------------
function makepos(x)
#-------------------------------------------------------------------------
    exp(x)
end


text:
using Accessors: @set
include("./util.jl")

struct Pos
    horiz::Int
    depth::Int
    aim::Int
end

struct Forward
    amt::Int
end

struct Down
    amt::Int
end

struct Up
    amt::Int
end

Command = Union{Forward,Down,Up}

function apply(pos::Pos, cmd::Forward)::Pos
    @set pos.horiz += cmd.amt
end

function apply(pos::Pos, cmd::Down)::Pos
    @set pos.depth += cmd.amt
end

function apply(pos::Pos, cmd::Up)::Pos
    @set pos.depth -= cmd.amt
end

function parse_cmd(s::AbstractString)::Command
    cmd, s_amt = split(s, ' ')
    amt = parse(Int, s_amt)

    if cmd == "forward"
        Forward(amt)
    elseif cmd == "up"
        Up(amt)
    elseif cmd == "down"
        Down(amt)
    else
        throw(DomainError("Unknown command $cmd"))
    end
end

inputs = util.input_lines(2)
commands = map(parse_cmd, inputs)

state = Pos(0, 0, 0)

for c in commands
    global state = apply(state, c)
end

println("Part 1: $(state.depth * state.horiz)")

function apply2(pos::Pos, cmd::Forward)::Pos
    Pos(pos.horiz + cmd.amt, pos.depth + pos.aim * cmd.amt, pos.aim)
end

function apply2(pos::Pos, cmd::Down)::Pos
    @set pos.aim += cmd.amt
end

function apply2(pos::Pos, cmd::Up)::Pos
    @set pos.aim -= cmd.amt
end

state2 = Pos(0, 0, 0)

for c in commands
    global state2 = apply2(state2, c)
end

println("Part 2: $(state2.depth * state2.horiz)")


text:
using Revise
using MathOptInterface
const MOI = MathOptInterface

export YasolVariable, YasolConstraint

# JUMP extensions

# variable extension
struct YasolVariable
    info::JuMP.VariableInfo
    quantifier::String
    block::Int64
end

function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    ::Type{YasolVariable};
    quantifier::String,
    block::Int64,
    kwargs...,
)
    return YasolVariable(
        info,
        quantifier,
        block,
    )
end

function JuMP.add_variable(
    model::JuMP.Model,
    yasolVar::YasolVariable,
    name::String,
)
    var = JuMP.add_variable(
            model,
            JuMP.ScalarVariable(yasolVar.info),
            name,
        )

    # add variable attributes to variable
    MOI.set(model, YasolSolver.VariableAttribute("quantifier"), var, yasolVar.quantifier)
    MOI.set(model, YasolSolver.VariableAttribute("block"), var, yasolVar.block)

    # print warning, if variable in first block is not existential
    if(yasolVar.block == 1 && yasolVar.quantifier != "exists")
        @error string("Variables in the first block need to be existential! Please add a dummy variable!")
        return
    end

    # check if quantifier is "exists" or "all"
    if((yasolVar.quantifier != "exists") && (yasolVar.quantifier != "all"))
        @error string("Variable quantifier has to be either 'exists' or 'all'!")
    end

    # check if block is an integer
    if(!isinteger(yasolVar.block))
        @error string("Variable blocks need to be of type integer!")
    end

    return var
end

# constraint extension
struct YasolConstraint
    f::AffExpr
    s::MOI.AbstractScalarSet
    quantifier::String
end

function JuMP.build_constraint(
    _error::Function,
    f::AffExpr,
    s::MOI.AbstractScalarSet,
    ::Type{YasolConstraint};
    quantifier::String,
)
    return YasolConstraint(f, s, quantifier)
end

function JuMP.add_constraint(
    model::Model,
    yasolCon::YasolConstraint,
    name::String,
)
    con = JuMP.add_constraint(
        model,
        ScalarConstraint(yasolCon.f, yasolCon.s),
        name,
    )

    # add constarint attributes to constraint
    MOI.set(model, YasolSolver.ConstraintAttribute("quantifier"), con, yasolCon.quantifier)

    # check if quantifier is "exists" or "all"
    if((yasolCon.quantifier != "exists") && (yasolCon.quantifier != "all"))
        @error string("Constraint quantifier has to be either 'exists' or 'all'!")
    end

    return con
end


text:
using Random
using GroupsCore

struct DirectProduct{Gt,Ht,GEl,HEl} <: GroupsCore.Group
    first::Gt
    last::Ht

    function DirectProduct(G::GroupsCore.Group, H::GroupsCore.Group)
        return new{typeof(G),typeof(H),eltype(G),eltype(H)}(G, H)
    end
end

struct DirectProductElement{GEl,HEl,Gt,Ht} <: GroupsCore.GroupElement
    elts::Tuple{GEl,HEl}
    parent::DirectProduct{Gt,Ht,GEl,HEl}
end

DirectProductElement(g, h, G::DirectProduct) = DirectProduct((g, h), G)

Base.one(G::DirectProduct) =
    DirectProductElement((one(G.first), one(G.last)), G)

Base.eltype(::Type{<:DirectProduct{Gt,Ht,GEl,HEl}}) where {Gt,Ht,GEl,HEl} =
    DirectProductElement{GEl,HEl,Gt,Ht}

function Base.iterate(G::DirectProduct)
    itr = Iterators.product(G.first, G.last)
    res = iterate(itr)
    @assert res !== nothing
    elt = DirectProductElement(first(res), G)
    return elt, (iterator = itr, state = last(res))
end

function Base.iterate(G::DirectProduct, state)
    itr, st = state.iterator, state.state
    res = iterate(itr, st)
    res === nothing && return nothing
    elt = DirectProductElement(first(res), G)
    return elt, (iterator = itr, state = last(res))
end

function Base.IteratorSize(::Type{<:DirectProduct{Gt,Ht}}) where {Gt,Ht}
    Gi = Base.IteratorSize(Gt)
    Hi = Base.IteratorSize(Ht)
    if Gi isa Base.IsInfinite || Hi isa Base.IsInfinite
        return Base.IsInfinite()
    elseif Gi isa Base.SizeUnknown || Hi isa Base.SizeUnknown
        return Base.SizeUnknown()
    else
        return Base.HasShape{2}()
    end
end

Base.size(G::DirectProduct) = (length(G.first), length(G.last))

GroupsCore.order(::Type{I}, G::DirectProduct) where {I<:Integer} =
    convert(I, order(I, G.first) * order(I, G.last))

GroupsCore.ngens(G::DirectProduct) = ngens(G.first) + ngens(G.last)

function GroupsCore.gens(G::DirectProduct)
    gens_first = [DirectProductElement((g, one(G.last)), G) for g in gens(G.first)]

    gens_last = [DirectProductElement((one(G.first), g), G) for g in gens(G.last)]

    return [gens_first; gens_last]
end

Base.isfinite(G::DirectProduct) = isfinite(G.first) && isfinite(G.last)

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{<:DirectProduct},
)
    G = rs[]
    return DirectProductElement((rand(rng, G.first), rand(rng, G.last)), G)
end

GroupsCore.parent(g::DirectProductElement) = g.parent
GroupsCore.parent_type(
    ::Type{DirectProductElement{GEl,HEl,Gt,Ht}},
) where {GEl,HEl,Gt,Ht} = DirectProduct{Gt,Ht,GEl,HEl}

Base.:(==)(g::DirectProductElement, h::DirectProductElement) =
    (parent(g) === parent(h) && g.elts == h.elts)

Base.hash(g::DirectProductElement, h::UInt) = hash(g.elts, hash(parent(g), h))

Base.deepcopy_internal(g::DirectProductElement, stackdict::IdDict) =
    DirectProductElement(Base.deepcopy_internal(g.elts, stackdict), parent(g))

Base.inv(g::DirectProductElement) =
    DirectProductElement(inv.(g.elts), parent(g))

function Base.:(*)(g::DirectProductElement, h::DirectProductElement)
    @assert parent(g) === parent(h)
    return DirectProductElement(g.elts .* h.elts, parent(g))
end

GroupsCore.order(::Type{I}, g::DirectProductElement) where {I<:Integer} =
    convert(I, lcm(order(I, first(g.elts)), order(I, last(g.elts))))

Base.isone(g::DirectProductElement) = all(isone, g.elts)

Base.show(io::IO, G::DirectProduct) =
    print(io, "Direct product of $(G.first) and $(G.last)")
Base.show(io::IO, g::DirectProductElement) =
    print(io, "( $(join(g.elts, ",")) )")


text:
using Test, PhysicalCommunications

testfiles = [ "prbs.jl", "eyediag.jl"]

for testfile in testfiles
	include(testfile)
end


text:
const USE_GPU = false
const do_visu = true
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using ImplicitGlobalGrid, Plots, Printf, LinearAlgebra
import MPI

# global reduction
norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))
# macros to avoid array allocation
macro qHx(ix,iy)  esc(:( -(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]))*(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]))*(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1])) * (H[$ix+1,$iy+1]-H[$ix,$iy+1])*_dx )) end
macro qHy(ix,iy)  esc(:( -(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]))*(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]))*(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1])) * (H[$ix+1,$iy+1]-H[$ix+1,$iy])*_dy )) end
macro dtau(ix,iy) esc(:(  (1.0/(min_dxy2 / (H[$ix+1,$iy+1]*H[$ix+1,$iy+1]*H[$ix+1,$iy+1]) / 4.1) + _dt)^-1  )) end

# @parallel_indices (ix,iy) function compute_update!(H2::Data.Array, dHdtau::Data.Array, H, Hold::Data.Array, _dt::Data.Number, damp::Data.Number, min_dxy2::Data.Number, _dx::Data.Number, _dy::Data.Number)
@parallel_indices (ix,iy) function compute_update!(H2, dHdtau, H, Hold, _dt, damp, min_dxy2, _dx, _dy)
    if (ix<=size(dHdtau,1) && iy<=size(dHdtau,2)) dHdtau[ix,iy] = -(H[ix+1, iy+1] - Hold[ix+1, iy+1])*_dt +
                                                                   (-(@qHx(ix+1,iy)-@qHx(ix,iy))*_dx -(@qHy(ix,iy+1)-@qHy(ix,iy))*_dy) +
                                                                   damp*dHdtau[ix,iy] end                       # damped rate of change
    if (ix<=size(dHdtau,1) && iy<=size(dHdtau,2)) H2[ix+1,iy+1] = H[ix+1,iy+1] + @dtau(ix,iy)*dHdtau[ix,iy] end # update rule, sets the BC as H[1]=H[end]=0
    return
end

# @parallel_indices (ix,iy) function compute_residual!(ResH::Data.Array, H::Data.Array, Hold::Data.Array, _dt::Data.Number, _dx::Data.Number, _dy::Data.Number)
@parallel_indices (ix,iy) function compute_residual!(ResH, H, Hold, _dt, _dx, _dy)
    if (ix<=size(ResH,1) && iy<=size(ResH,2)) ResH[ix,iy] = -(H[ix+1, iy+1] - Hold[ix+1, iy+1])*_dt +
                                                             (-(@qHx(ix+1,iy)-@qHx(ix,iy))*_dx -(@qHy(ix,iy+1)-@qHy(ix,iy))*_dy) end
    return
end

# @parallel_indices (ix,iy) function assign!(Hold::Data.Array, H::Data.Array)
@parallel_indices (ix,iy) function assign!(Hold, H)
    if (ix<=size(H,1) && iy<=size(H,2)) Hold[ix,iy] = H[ix,iy] end
    return
end

@views function diffusion_2D_damp_perf_multixpu()
    # Physics
    lx, ly = 10.0, 10.0   # domain size
    ttot   = 1.0          # total simulation time
    dt     = 0.2          # physical time step
    # Numerics
    # nx, ny = 16*32*16, 16*32*16 # numerical grid resolution
    nx, ny = 32, 32 # numerical grid resolution
    nout   = 100          # check error every nout
    tol    = 1e-6         # tolerance
    itMax  = 1e5          # max number of iterations
    # Derived numerics
    me, dims = init_global_grid(nx, ny, 1)  # MPI initialisation
    @static if USE_GPU select_device() end  # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy = lx/nx_g(), ly/ny_g()           # grid size
    damp   = 1-35/nx_g()                    # damping (this is a tuning parameter, dependent on e.g. grid resolution)
    # Array allocation
    ResH   = @zeros(nx-2, ny-2) # normal grid, without boundary points
    dHdtau = @zeros(nx-2, ny-2) # normal grid, without boundary points
    H      = @zeros(nx  , ny  )
    # Initial condition
    H     .= Data.Array([exp(-(x_g(ix,dx,H)+dx/2 -lx/2)^2 -(y_g(iy,dy,H)+dy/2 -ly/2)^2) for ix=1:size(H,1), iy=1:size(H,2)])
    Hold   = copy(H)
    H2     = copy(H)
    _dx, _dy, _dt = 1.0/dx, 1.0/dy, 1.0/dt
    min_dxy2 = min(dx,dy)^2
    len_ResH_g = ((nx-2-2)*dims[1]+2)*((ny-2-2)*dims[2]+2)
    # Preparation of visualisation
    if do_visu
        if me==0
            ENV["GKSwstype"]="nul"; if isdir("viz2D_xpu_out")==false mkdir("viz2D_xpu_out") end; loadpath = "./viz2D_xpu_out/"; anim = Animation(loadpath,String[])
            println("Animation directory: $(anim.dir)")
        end
        nx_v, ny_v = (nx-2)*dims[1], (ny-2)*dims[2]
        if (nx_v*ny_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        H_v   = zeros(nx_v, ny_v) # global array for visu
        H_inn = zeros(nx-2, ny-2) # no halo local array for visu
        Xi_g, Yi_g = LinRange(dx+dx/2, lx-dx-dx/2, nx_v), LinRange(dy+dy/2, ly-dy-dy/2, ny_v) # inner points only
    end
    t = 0.0; it = 0; ittot = 0; t_tic = 0.0; niter = 0
    # Physical time loop
    while t<ttot
        iter = 0; err = 2*tol
        # Picard-type iteration
        while err>tol && iter<itMax
            if (it==1 && iter==0) t_tic = Base.time(); niter = 0 end
            @hide_communication (8, 4) begin # communication/computation overlap
                @parallel compute_update!(H2, dHdtau, H, Hold, _dt, damp, min_dxy2, _dx, _dy)
                H, H2 = H2, H
                update_halo!(H)
            end
            if iter % nout == 0
                @parallel compute_residual!(ResH, H, Hold, _dt, _dx, _dy)
                err = norm_g(ResH)/len_ResH_g
            end
            iter += 1; niter += 1
        end
        ittot += iter; it += 1; t += dt
        @parallel assign!(Hold, H)
        # Visualize
        if do_visu
            H_inn .= H[2:end-1,2:end-1]; gather!(H_inn, H_v)
            if (me==0)
                fontsize = 12
                opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
                        ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
                        xlabel="Lx", ylabel="Ly", xlims=(Xi_g[1], Xi_g[end]), ylims=(Yi_g[1], Yi_g[end]), clims=(0.,1.))
                heatmap(Xi_g, Yi_g, Array(H_v)'; c=:davos, title="diffusion MPI (nt=$it)", opts...); frame(anim)
            end
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (2*2+1)/1e9*nx_g()*ny_g()*sizeof(Data.Number) # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                           # Execution time per iteration [s]
    T_eff = A_eff/t_it                            # Effective memory throughput [GB/s]
    if (me==0) @printf("Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=2), niter) end
    if (do_visu && me==0) gif(anim, "diffusion_2D_multixpu.gif", fps = 5)  end
    finalize_global_grid()
    return
end

diffusion_2D_damp_perf_multixpu()


text:
function compressLZW(decompressed::String)
    dictsize = 256
    dict     = Dict{String,Int}(string(Char(i)) => i for i in 0:dictsize)
    result   = Vector{Int}(undef, 0)
    w        = ""
    for c in decompressed
        wc = string(w, c)
        if haskey(dict, wc)
            w = wc
        else
            push!(result, dict[w])
            dict[wc]  = dictsize
            dictsize += 1
            w        = string(c)
        end
    end
    if !isempty(w) push!(result, dict[w]) end
    return result
end

function decompressLZW(compressed::Vector{Int})
    dictsize = 256
    dict     = Dict{Int,String}(i => string('\0' + i) for i in 0:dictsize)
    result   = IOBuffer()
    w        = string(Char(popfirst!(compressed)))
    write(result, w)
    for k in compressed
        if haskey(dict, k)
            entry = dict[k]
        elseif k == dictsize
            entry = string(w, w[1])
        else
            error("bad compressed k: $k")
        end
        write(result, entry)
        dict[dictsize] = string(w, entry[1])
        dictsize += 1
        w = entry
    end
    return String(take!(result))
end

original     = ["0123456789", "TOBEORNOTTOBEORTOBEORNOT", "dudidudidudida"]
compressed   = compressLZW.(original)
decompressed = decompressLZW.(compressed)

for (word, comp, decomp) in zip(original, compressed, decompressed)
    comprate = (length(word) - length(comp)) / length(word) * 100
    println("Original: $word\n-> Compressed: $comp (compr.rate: $(round(comprate, digits=2))%)\n-> Decompressed: $decomp")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import RestrictedBoltzmannMachines as RBMs
using Test: @test, @testset, @inferred
using Statistics: mean
using LinearAlgebra: dot
using Random: bitrand, rand!
using Zygote: gradient, jacobian
using RestrictedBoltzmannMachines: visible, hidden, weights, free_energy, energy, interaction_energy
using RestrictedBoltzmannMachines: RBM, BinaryRBM, grad2ave, subtract_gradients
using RestrictedBoltzmannMachines: Binary, Spin, Potts, Gaussian, ReLU, dReLU, pReLU, xReLU
using RestrictedBoltzmannMachines: sample_v_from_v, sample_h_from_h, ∂free_energy

struct CenteredRBM{V<:Binary,H<:Binary,W<:AbstractArray,Av<:AbstractArray,Ah<:AbstractArray}
    visible::V
    hidden::H
    w::W
    λv::Av
    λh::Ah
end

CenteredRBM(rbm::RBM, λv::Vector, λh::Vector) = CenteredRBM(
    visible(rbm), hidden(rbm), weights(rbm), λv, λh
)

CenteredRBM(a::Vector, b::Vector, w::Matrix, λv::Vector, λh::Vector) = CenteredRBM(
    BinaryRBM(a, b, w), λv, λh
)

RBM(rbm::CenteredRBM) = BinaryRBM(rbm.visible.θ, rbm.hidden.θ, rbm.w)

center(rbm::RBM{<:Binary,<:Binary}, λv::Vector, λh::Vector) = center(
    rbm.visible.θ, rbm.hidden.θ, rbm.w, λv, λh
)
center(a::Vector, b::Vector, w::Matrix, λv::Vector, λh::Vector) = CenteredRBM(
    a + w * λh, b + w' * λv, w, λv, λh
)
uncenter(rbm::CenteredRBM) = uncenter(rbm.visible.θ, rbm.hidden.θ, rbm.w, rbm.λv, rbm.λh)
uncenter(a::Vector, b::Vector, w::Matrix, λv::Vector, λh::Vector) = BinaryRBM(
    a - w * λh, b - w' * λv, w
)

function RBMs.energy(rbm::CenteredRBM, v::AbstractArray, h::AbstractArray)
    Ev = energy(rbm.visible, v)
    Eh = energy(rbm.hidden, h)
    Ew = interaction_energy(RBM(rbm), v .- rbm.λv, h .- rbm.λh)
    return Ev .+ Eh .+ Ew
end

function RBMs.free_energy(rbm::CenteredRBM, v::AbstractArray)
    inputs = RBMs.inputs_h_from_v(RBM(rbm), v .- rbm.λv)
    E = energy(rbm.visible, v) + energy(Binary(-rbm.λh), inputs)
    F = free_energy(rbm.hidden, inputs)
    return E + F
end

ΔE(rbm::CenteredRBM) = interaction_energy(RBM(rbm), rbm.λv, rbm.λh)

@testset "centered RBM" begin
    rbm = BinaryRBM(randn(3), randn(2), randn(3,2))
    rbm_c = center(rbm, randn(3), randn(2))
    v = bitrand(3, 100)
    h = bitrand(2, 100)
    @test energy(rbm_c, v, h) ≈ energy(rbm, v, h) .+ ΔE(rbm_c)
    @test free_energy(rbm_c, v) ≈ free_energy(rbm, v) .+ ΔE(rbm_c)
    hs = [[0,0], [0,1], [1,0], [1,1]]
    @test exp.(-free_energy(rbm_c, v)) ≈ sum(exp.(-energy(rbm_c, v, h)) for h in hs)

    rbm_u = uncenter(rbm_c)
    @test rbm_u.visible.θ ≈ rbm.visible.θ
    @test rbm_u.hidden.θ ≈ rbm.hidden.θ
    @test rbm_u.w ≈ rbm.w
end

@testset "center_gradient" begin
    rbm = BinaryRBM(randn(3), randn(2), randn(3,2))
    λv = randn(3)
    λh = randn(2)
    rbmc = center(rbm, λv, λh)
    v = bitrand(3, 1000)

    # subtraction eliminates constant energy displacement
    ∂, = gradient(rbm) do rbm
        mean(free_energy(rbm, v)) - free_energy(rbm, falses(3))
    end
    ∂c = RBMs.center_gradient(rbm, ∂, λv, λh)

    gs = gradient(rbmc) do rbmc
        mean(free_energy(rbmc, v)) - free_energy(rbmc, falses(3))
    end

    @test ∂c.visible.θ ≈ only(gs).visible.θ
    @test ∂c.hidden.θ ≈ only(gs).hidden.θ
    @test ∂c.w ≈ only(gs).w
end

@testset "centered gradient" begin
    rbm = BinaryRBM(randn(3), randn(2), randn(3,2))
    λv = randn(3)
    λh = randn(2)
    rbmc = center(rbm, λv, λh)
    v = bitrand(3, 1000)

    # subtraction eliminates constant energy displacement
    ∂ = subtract_gradients(∂free_energy(rbm, v), ∂free_energy(rbm, falses(3)))
    ∂c = @inferred RBMs.centered_gradient(rbm, ∂, rbmc.λv, rbmc.λh)

    ad, = gradient(rbmc) do rbmc
        mean(free_energy(rbmc, v)) - free_energy(rbmc, falses(3))
    end

    # since the centering transformation is linear, we can just uncenter the gradient
    ad_u = uncenter(ad.visible.θ, ad.hidden.θ, ad.w, λv, λh)

    @test ∂c.visible.θ ≈ ad_u.visible.θ
    @test ∂c.hidden.θ ≈ ad_u.hidden.θ
    @test ∂c.w ≈ ad_u.w
end

@testset "centered gradient via jacobian" begin
    rbm = BinaryRBM(randn(3), randn(2), randn(3,2))
    λv = randn(3)
    λh = randn(2)
    rbmc = center(rbm, λv, λh)
    v = bitrand(3, 1000)

    # subtraction eliminates constant energy displacement
    ∂ = subtract_gradients(∂free_energy(rbm, v), ∂free_energy(rbm, falses(3)))
    ∂c = @inferred RBMs.centered_gradient(rbm, ∂, rbmc.λv, rbmc.λh)

    ad, = gradient(rbmc) do rbmc
        mean(free_energy(rbmc, v)) - free_energy(rbmc, falses(3))
    end

    Ja, Jb, Jw = jacobian(rbmc.visible.θ, rbmc.hidden.θ, rbmc.w) do a, b, w
        uncenter(a, b, w, λv, λh).visible.θ
    end
    @test ∂c.visible.θ ≈ Ja * ad.visible.θ + Jb * ad.hidden.θ + Jw * vec(ad.w)

    Ja, Jb, Jw = jacobian(rbmc.visible.θ, rbmc.hidden.θ, rbmc.w) do a, b, w
        uncenter(a, b, w, λv, λh).hidden.θ
    end
    @test ∂c.hidden.θ ≈ Ja * ad.visible.θ + Jb * ad.hidden.θ + Jw * vec(ad.w)

    Ja, Jb, Jw = jacobian(rbmc.visible.θ, rbmc.hidden.θ, rbmc.w) do a, b, w
        uncenter(a, b, w, λv, λh).w
    end
    @test vec(∂c.w) ≈ Ja * ad.visible.θ + Jb * ad.hidden.θ + Jw * vec(ad.w)
end

"""
    struct2nt(obj)

Converts a struct to a NamedTuple with the same fields.
"""
struct2nt(s) = NamedTuple{propertynames(s)}(([getproperty(s, p) for p in propertynames(s)]...,))
@test struct2nt(1 + 2im) == (re = 1, im = 2)
@test struct2nt(Gaussian(ones(5), ones(5))) == (θ = ones(5), γ = ones(5))

_layers = (Binary, Spin, Potts, Gaussian, ReLU, dReLU, pReLU, xReLU)

@testset "shifting $Layer gradients" for Layer in _layers
    layer = Layer(5)
    ∂ = struct2nt(layer)
    λ = randn(size(layer))
    applicable(RBMs.uncenter_step, layer, ∂, λ) || continue
    for p in propertynames(layer)
        rand!(getproperty(layer, p))
    end
    ∂c = RBMs.uncenter_step(layer, ∂, λ)
    layerc = deepcopy(layer)
    for p in propertynames(layer)
        getproperty(layerc, p) .= getproperty(∂c, p)
    end
    x = randn(size(layer)..., 10)
    @test energy(layer, x) ≈ energy(layerc, x) - x' * λ
end

@testset "grad2ave $Layer" for Layer in _layers
    layer = Layer(5)
    rbm = RBM(layer, Binary(randn(3)), randn(5,3))
    v = sample_v_from_v(rbm, randn(5,100); steps=100)
    ∂ = ∂free_energy(rbm, v)
    @test grad2ave(visible(rbm), ∂.visible) ≈ dropdims(mean(v; dims=2); dims=2)
end


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

function code_coverage_test()
    if rand(1:2) == 3
        return "hello" # never reached
    else
        r = Int[]
        for i = 1:3
            push!(r, i)
        end
        nothing
        return r
    end
    not_reached
end

@eval short_form_func_coverage_test(x) = begin
    $(Expr(:line, 1234))
    y() = begin
        x
    end
    x * y()
end

success = code_coverage_test() == [1, 2, 3] &&
          short_form_func_coverage_test(2) == 4
exit(success ?  0 : 1)

# end of file

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using BinaryBuilder

name = "MKL"
version = v"2020.1.216"

sources = [
    ArchiveSource("https://anaconda.org/intel/mkl/2020.1/download/linux-64/mkl-2020.1-intel_217.tar.bz2",
                  "8814e952c0b4f28079361adac8bec1051b97d62dee621666798a3302e70d75e0"; unpack_target = "mkl-x86_64-linux-gnu"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.1/download/osx-64/mkl-2020.1-intel_216.tar.bz2",
                  "4dab8dd1aa12b02cd121228ba881a887448399c799ee5a2d53942716a03f880e"; unpack_target = "mkl-x86_64-apple-darwin14"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.1/download/win-32/mkl-2020.1-intel_216.tar.bz2",
                  "b43ca7a8f5aed51e197af1f9165b3a72e5dd0a7f036fd9364ccfe90a8645a235"; unpack_target = "mkl-i686-w64-mingw32"),
    ArchiveSource("https://anaconda.org/intel/mkl/2020.1/download/win-64/mkl-2020.1-intel_216.tar.bz2",
                  "5dd8eff29b390ddb1cf15f391aba7fa3bdc4034f79b98d226865a5e331060f76"; unpack_target = "mkl-x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/mkl-${target}
if [[ ${target} == *mingw* ]]; then
    cp -r Library/bin/* ${libdir}
else
    cp -r lib/* ${libdir}
fi
install_license info/*.txt
"""

platforms = [
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libmkl_core", "mkl_core"], :libmkl_core),
    LibraryProduct(["libmkl_rt", "mkl_rt"], :libmkl_rt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("IntelOpenMP_jll"),
]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)
include("../../fancy_toys.jl")
no_autofix_platforms = [Windows(:i686), Windows(:x86_64), MacOS(:x86_64)]
autofix_platforms = [Linux(:x86_64)]
if any(should_build_platform.(triplet.(no_autofix_platforms)))
    # Need to disable autofix: updating linkage of libmkl_intel_thread.dylib on
    # macOS causes runtime issues:
    # https://github.com/JuliaPackaging/Yggdrasil/issues/915.
    build_tarballs(non_reg_ARGS, name, version, sources, script, no_autofix_platforms, products, dependencies; lazy_artifacts = true, autofix = false)
end
if any(should_build_platform.(triplet.(autofix_platforms)))
    # ... but we need to run autofix on Linux, because here libmkl_rt doesn't
    # have a soname, so we can't ccall it without specifying the path:
    # https://github.com/JuliaSparse/Pardiso.jl/issues/69
    build_tarballs(ARGS, name, version, sources, script, autofix_platforms, products, dependencies; lazy_artifacts = true)
end


text:
module Lambda
export LEval

function LEval(str) #... it's unreadable, why did I wrote it this way???
    o = split(str, r"\r?\n") .|> String
    if length(o) == 0
        return
    elseif length(o) == 1
        s = o[1]
        if (r = match(r":=", s)) != nothing
            i = r.offset
            return Expr(:(=), Symbol(s[1:prevind(o[1], i)]), LEval(s[nextind(o[1], i, 2):end]))
        end
        if (r = match(r"\(", s)) != nothing
            if (r2 = match(r"λ.*\[.*\]", s)) != nothing
                if r2.offset < r.offset
                    @goto END_P
                end
            end
            i = r.offset
            u = mfc(s, ')', '(', i)
            z = Expr(:block, LEval(s[nextind(s, i):prevind(s, u)]))
            if u != lastindex(s)
                z = call_AAAAA(deepcopy(z), LEval(s[nextind(s, u):end]))
            end
            if i != firstindex(s)
                z = call_AAAAA(LEval(s[1:prevind(s, i)]), deepcopy(z))
            end
            return z
            @label END_P
        end
        if (r = match(r"λ.*\[.*\]", s)) != nothing
            i = r.offset
            u = mff(s, '[', i+1)
            e = mfc(s, ']', '[', u)
            ex = Expr(:->, Symbol(s[nextind(s, i):prevind(s, u)]), LEval(s[nextind(s, u):prevind(s, e)]))
            if e == lastindex(s)
                return ex
            else
                z = LEval(s[nextind(s, e):end])
                return call_AAAAA(ex, z) # too much mental pain from that function, thus the name
            end
        end
        if length(split(s, " ")) == 1
            return Meta.parse(idk(s))
        else
            throw("owo not implemedeeded") # coz idk
        end
    else
        println("not implmeented lolr")
        sleep(2)
        println("just joking :P")
        LEval.(o)
    end
end

function call_AAAAA(ex, z::Expr)
    l = 0
    if z.head == :call
        return Expr(:call, call_AAAAA(ex, z.args[1]), z.args[2])
    else
        return Expr(:call, ex, z)
    end
end

function call_AAAAA(ex, z)
    Expr(:call, ex, z)
end

macro LEval_str(str)
    LEval(str)
end

function mfc(str, c, n, i = 1)
    o = 0
    for d in eachindex(str)
        if d >= i
            if str[d] == n
                o -= 1
            elseif str[d] == c
                o += 1
                if o == 0
                    return d
                end
            end
        end
    end
    missing
end
function mff(str, c, i = 1)
    for d in eachindex(str)
        if d >= i
            str[d] == c && return d
        end
    end
    missing
end

function idk(s)
    o = ""
    for i in s
        if i ∉ "()"
            o = o * "($i)"
        else
            o = o * "$i"
        end
    end
    o
end

end


text: Teardown = Tuple{Base.Generator, UInt8}


text:
module UnitTests

using Test
using PBC

using PBC.Curve: BN, FP, EP, EP2, curve_gen, LIB
using PBC.Config: ORDER, G1, G2, @EP, @EP2
using PBC.Shamir: evalpoly, LagrangeCoeffGenerator, BarycentricWeightGenerator, lagrange_interpolate_c0

@testset "lagrange_interpolate_c0 - trivial" begin
    ids = 1:2
    for c in (BN(1), curve_gen(EP), curve_gen(EP2))
        poly = [c, c] # derive secret coeffs from c
        shares = Dict(x=>evalpoly(poly, x) for x in ids)
        @test lagrange_interpolate_c0(BarycentricWeightGenerator(keys(shares)), shares) == c
        @test lagrange_interpolate_c0(LagrangeCoeffGenerator(keys(shares)), shares) == c
        @test lagrange_interpolate_c0(shares) == c
    end
end

@testset "Shamir: lagrange_interpolate_c0 - complex" begin
    ids = 1:5
    n = length(ids)
    for c in (BN(123), 123 * curve_gen(EP), 123 * curve_gen(EP2))
        poly = [c, 2c, 3c] # derive secret coeffs from c
        t = length(poly)
        shares = Dict(x=>evalpoly(poly, x) for x in ids)
        for i in 1:n-t+1
            r = i:i+t-1
            @test lagrange_interpolate_c0(BarycentricWeightGenerator(keys(shares)), shares) == c
            @test lagrange_interpolate_c0(LagrangeCoeffGenerator(keys(shares)), shares) == c
            @test lagrange_interpolate_c0(shares) == c
        end
    end
end

@testset "Shamir: LagrangeCoeffGenerator" begin
    coeffs = LagrangeCoeffGenerator()
    # generate first coeff
    @test push!(coeffs, 1) == BN(1)

    # generate second coeff
    @test push!(coeffs, 2) == ORDER - 1
    @test push!(coeffs, 1) == BN(2)

    # generate third share
    @test push!(coeffs, 3) == BN(1)
    @test push!(coeffs, 2) == ORDER - 3
    @test push!(coeffs, 1) == BN(3)

    # delete third share
    delete!(coeffs, 3)
    @test length(coeffs) == 2
    @test push!(coeffs, 2) == ORDER - 1
    @test push!(coeffs, 1) == BN(2)
end

@testset "Shamir: BarycentricWeightGenerator" begin
    weights = BarycentricWeightGenerator()
    # generate first coeff
    @test push!(weights, 1) == ORDER - 1

    # generate second coeff
    @test push!(weights, 2) == BN(2)
    @test push!(weights, 1) == ORDER - 1

    # generate third share
    @test push!(weights, 3) == ORDER - 6
    @test push!(weights, 2) == BN(2)
    @test push!(weights, 1) == ORDER - 2

    # delete third share
    delete!(weights, 3)
    @test length(weights) == 2
    @test push!(weights, 2) == BN(2)
    @test push!(weights, 1) == ORDER - 1
end

@testset "PrivateKey" begin
    # generate random with default RNG
    @test rand(PrivateKey).sk < PBC.Config.PRIME

    # from bytes
    sk = PrivateKey(UInt8[1, 2, 3])
    @test sk.sk < PBC.Config.PRIME
    @test sk.sk.dp[1] == 1 << 16 | 2 << 8 | 3
    @test !iszero(sk) && isvalid(sk)

    # serialization
    @test length(Vector{UInt8}(sk)) == PBC.Config.PRIVATE_KEY_SIZE
    @test PrivateKey(Vector{UInt8}(sk)) == sk
end

@testset "PublicKey" begin
    # generate from secret
    sk = PrivateKey(UInt8[1, 2, 3])
    pk = PublicKey(sk)
    @test !isinf(pk) && isvalid(pk)

    # serialization
    @test length(Vector{UInt8}(pk)) == PBC.Config.PUBLIC_KEY_SIZE
    @test PublicKey(Vector{UInt8}(pk)) == pk
end

@testset "Hash" begin
    hash = Hash("foo")
    @test !isinf(hash) && isvalid(hash)
end

@testset "Identity" begin
    @test !signbit(Identity(-1).id)
    @test !iszero(Identity(PublicKey(PrivateKey(UInt8[1, 2, 3]))).id)
end

@testset "PrivateKeyPoly" begin
    poly = PrivateKeyPoly(2)
    @test length(poly) == 3
    # test iterate
    @test length(collect(poly)) == 3
    # test index
    @test isa(poly[1], BN)
end

@testset "PublicKeyPoly" begin
    skpoly = PrivateKeyPoly(2)
    pkpoly = PublicKeyPoly(skpoly)
    @test length(skpoly) == length(pkpoly)
    # test iterate
    @test length(collect(pkpoly)) == 3
    # test index
    @test isa(pkpoly[1], @EP2)
end

@testset "create_share & validate_share" begin
    # create a secret polynomial
    pk1, pk2 = PublicKey(PrivateKey(UInt8[1, 2, 3])), PublicKey(PrivateKey(UInt8[2, 3, 4]))
    skpoly = PrivateKeyPoly(10)
    share1, share2 = PBC.create_share(skpoly, pk1), PBC.create_share(skpoly, pk2)
    @test share1 != share2
    @test isa(share1, PBC.Model.AbstractPrivateKey)

    pkpoly = PublicKeyPoly(skpoly)
    @test PBC.verify_share(pkpoly, pk1, share1)
    @test PBC.verify_share(pkpoly, pk2, share2)
    # TODO: explore if this can be made to work?
    @test_broken PBC.verify_share(pkpoly, pk1+pk2, share1+share2)
end

@testset "batch validate_share" begin
    pk = PublicKey(PrivateKey(UInt8[1, 2, 3]))
    # create a secret polynomial
    skpoly1, skpoly2 = PrivateKeyPoly(10), PrivateKeyPoly(10)
    # create shares for the same id from two different secret polynomials
    share1, share2 = PBC.create_share(skpoly1, pk), PBC.create_share(skpoly2, pk)
    pkpoly1, pkpoly2 = PublicKeyPoly(skpoly1), PublicKeyPoly(skpoly2)
    # make sure that the aggregate secret share validates with the aggregate public key polynomial
    @test PBC.verify_share(pkpoly1+pkpoly2, pk, share1+share2)
end

@testset "share sign & recover" begin
    alice, bob, charlotte  = PublicKey(rand(PrivateKey)), PublicKey(rand(PrivateKey)), PublicKey(rand(PrivateKey))
    # create a secret polynomial of degree 2 => 3 shares
    skpoly = PrivateKeyPoly(2)
    pkpoly = PublicKeyPoly(skpoly)
    # create three private keys from the secret polynomial
    skalice, skbob, skcharlotte = PBC.create_share(skpoly, alice), PBC.create_share(skpoly, bob), PBC.create_share(skpoly, charlotte)
    # the signature shares are collected in a SignatureShares object
    sigshares = SignatureShares()
    # each user signs with each secret key share, thus creates separate signature shares
    sigshares[alice] = PBC.sign(skalice, "foo")
    sigshares[bob] = PBC.sign(skbob, "foo")
    # not enough signatures - the recovered signature is NOT valid
    @test !PBC.verify(Signature(sigshares), PublicKey(pkpoly), "foo")
    # enough signatures - the recovered signature is valid
    sigshares[charlotte] = PBC.sign(skcharlotte, "foo")
    @test PBC.verify(Signature(sigshares), PublicKey(pkpoly), "foo")
end

@testset "sign & verify" begin
    sk = PrivateKey(UInt8[1, 2, 3])
    sig = PBC.sign(sk, "foo")
    @test !isinf(sig) && isvalid(sig)
    @test PBC.verify(sig, PublicKey(sk), "foo")
    @test !PBC.verify(sig, PublicKey(sk), "bar")
end

@testset "simple aggregation" begin
    sk1, sk2, sk3 = PrivateKey(UInt8[1, 2, 3]), PrivateKey(UInt8[2, 3, 4]), PrivateKey(UInt8[3, 4, 5])
    pk1, pk2, pk3 = PublicKey(sk1), PublicKey(sk2), PublicKey(sk3)
    sig1, sig2, sig3 = PBC.sign(sk1, "foo"), PBC.sign(sk2, "foo"), PBC.sign(sk3, "foo")
    @test PBC.verify(sig1 + sig2, pk1 + pk2, "foo")
    @test !PBC.verify(sig2, pk1 + pk2, "foo")
    @test pk1 - pk1 + pk2 == pk1 + pk2 - pk1
    @test PBC.verify(sig2, pk1 - pk1 + pk2, "foo")
end

@testset "advanced aggregation" begin
    sk1, sk2, sk3 = PrivateKey(UInt8[1, 2, 3]), PrivateKey(UInt8[2, 3, 4]), PrivateKey(UInt8[3, 4, 5])
    pk1, pk2, pk3 = PublicKey(sk1), PublicKey(sk2), PublicKey(sk3)
    msg1, msg2, msg3 = "foo", "bar", "baz"
    sig1, sig2, sig3 = PBC.sign(sk1, msg1), PBC.sign(sk2, msg2), PBC.sign(sk3, msg3)
    @test PBC.verify(sig1 + sig2, (pk1 => Hash(msg1), pk2 => Hash(msg2)))
    @test !PBC.verify(sig1 + sig2, (pk1 => Hash(msg1),))
end

@testset "PrivateKey aggregation" begin
    sk1, sk2 = PrivateKey(UInt8[1, 2, 3]), PrivateKey(UInt8[2, 3, 4])
    pk1, pk2 = PublicKey(sk1), PublicKey(sk2)
    @test PBC.sign(sk1 + sk2, "foo") == PBC.sign(sk1, "foo") + PBC.sign(sk2, "foo")
    @test iszero((sk1 - sk1).sk)
end

@testset "invalid aggregation" begin
    sk = PrivateKey(UInt8[1, 2, 3])
    @test_throws ErrorException PBC.verify(PBC.sign(sk, "foo"), [])
end

@testset "Shamir: evalpoly" begin
    for c in (BN(123), curve_gen(EP), curve_gen(EP2))
        poly = [c, 2c, 3c]
        @test evalpoly(poly, 0) == c
        @test evalpoly(poly, 1) == c + 2c + 3c
        @test evalpoly(poly, 2) == c + 2c * 2 + 3c * 2^2
    end
end

end


text:
using MOCNeutronTransport
using HDF5
@testset "XDMF" begin
    @testset "c5g7 pin - triangles" begin
        vtk_to_xdmf_type = Dict(
            # triangle
            5  => 4,
            # triangle6
            22 => 36,
            # quadrilateral
            9 => 5,
            # quad8
            23 => 37
           )
        ref_points = [Point_2D(0.0, 0.0),
                      Point_2D(1.26, 0.0),
                      Point_2D(0.0, 1.26),
                      Point_2D(1.26, 1.26),
                      Point_2D(1.17, 0.63),
                      Point_2D(0.63, 0.0),
                      Point_2D(0.0, 0.63),
                      Point_2D(1.26, 0.63),
                      Point_2D(0.63, 1.26),
                      Point_2D(0.96668449300372, 1.0521890005327),
                      Point_2D(0.50983869566359, 1.1564610725782),
                      Point_2D(0.14347681133269, 0.86429721912348),
                      Point_2D(0.14347681133269, 0.39570278087652),
                      Point_2D(0.50983869566359, 0.10353892742181),
                      Point_2D(0.96668449300372, 0.20781099946726),
                      Point_2D(0.63, 0.63)]
        ref_faces = [
                     SVector{4, UInt64}(UInt64[5, 12, 16, 11]),
                     SVector{4, UInt64}(UInt64[5, 13, 16, 12]),
                     SVector{4, UInt64}(UInt64[5, 14, 16, 13]),
                     SVector{4, UInt64}(UInt64[5, 11, 16, 10]),
                     SVector{4, UInt64}(UInt64[5, 10, 16, 5]),
                     SVector{4, UInt64}(UInt64[5, 15, 16, 14]),
                     SVector{4, UInt64}(UInt64[5, 5, 16, 15]),
                     SVector{4, UInt64}(UInt64[5, 7, 1, 13]),
                     SVector{4, UInt64}(UInt64[5, 7, 12, 3]),
                     SVector{4, UInt64}(UInt64[5, 1, 6, 14]),
                     SVector{4, UInt64}(UInt64[5, 1, 14, 13]),
                     SVector{4, UInt64}(UInt64[5, 6, 2, 15]),
                     SVector{4, UInt64}(UInt64[5, 15, 2, 8]),
                     SVector{4, UInt64}(UInt64[5, 3, 11, 9]),
                     SVector{4, UInt64}(UInt64[5, 3, 12, 11]),
                     SVector{4, UInt64}(UInt64[5, 10, 8, 4]),
                     SVector{4, UInt64}(UInt64[5, 9, 10, 4]),
                     SVector{4, UInt64}(UInt64[5, 10, 5, 8]),
                     SVector{4, UInt64}(UInt64[5, 15, 8, 5]),
                     SVector{4, UInt64}(UInt64[5, 14, 6, 15]),
                     SVector{4, UInt64}(UInt64[5, 7, 13, 12]),
                     SVector{4, UInt64}(UInt64[5, 11, 10, 9])
                   ]
        ref_face_sets = Dict{String, Set{UInt64}}()
        ref_face_sets["MATERIAL_UO2-3.3"] = Set([1, 2, 3, 4, 5, 6, 7])
        ref_face_sets["GRID_L1_1_1"] = Set([1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                                            11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22])
        ref_face_sets["MATERIAL_WATER"] = Set([8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
                                               18, 19, 20, 21, 22])
        ref_face_sets["test_set"] = Set([1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
                                         13, 14, 15, 16, 17, 18, 19, 20, 21, 22])
        mesh = UnstructuredMesh_2D{Float64, UInt64}(
                                   points = ref_points,
                                   faces = ref_faces,
                                   face_sets = ref_face_sets,
                                   name = "c5g7_UO2_pin")
        # write_xdmf
        # redirect stdout to keep test results pretty
        redirect_stdout(devnull) do
            write_xdmf_2d("c5g7_UO2_pin.xdmf", mesh)
        end
        # check xdmf
        ref_file = open("./mesh/mesh_files/c5g7_UO2_pin.xdmf", "r")
        test_file = open("./c5g7_UO2_pin.xdmf", "r")
        while !eof(ref_file)
            ref_line = readline(ref_file)
            test_line = readline(test_file)
            @test ref_line == test_line
        end
        close(ref_file)
        close(test_file)
        # check h5
        h5_file = h5open("c5g7_UO2_pin.h5","r")

        # points
        point_array = zeros(Float64, length(mesh.points), 3)
        npoints = length(mesh.points)
        for i = 1:npoints
            point_array[i, 1] = mesh.points[i][1]
            point_array[i, 2] = mesh.points[i][2]
        end
        point_array = copy(transpose(point_array))
        test_points = read(h5_file["c5g7_UO2_pin"]["points"])
        @test test_points == point_array

        # cells
        topo_array = UInt64[]
        for face in mesh.faces
            # convert face to vector for mutability
            face_xdmf = [x for x in face]
            # adjust vtk to xdmf type
            face_xdmf[1] = vtk_to_xdmf_type[face_xdmf[1]]
            # adjust 1-based to 0-based indexing
            face_xdmf[2:length(face_xdmf)] = face_xdmf[2:length(face_xdmf)] .- 1
            for value in face_xdmf
                push!(topo_array, value)
            end
        end
        test_cells = read(h5_file["c5g7_UO2_pin"]["cells"])
        @test test_cells == topo_array

        # Cell sets
        test_grid_l1_1_1 = read(h5_file["c5g7_UO2_pin"]["GRID_L1_1_1"])
        ref_grid_l1_1_1 = [ x - 1 for x in ref_face_sets["GRID_L1_1_1"]]
        @test test_grid_l1_1_1 == ref_grid_l1_1_1
        test_test_set = read(h5_file["c5g7_UO2_pin"]["test_set"])
        ref_test_set = [ x - 1 for x in ref_face_sets["test_set"]]
        @test test_test_set == ref_test_set

        # Materials
        test_material_IDs = read(h5_file["c5g7_UO2_pin"]["material_id"])
        ref_material_IDs = [ i <= 7 ? 0 : 1 for i = 1:22 ]
        @test test_material_IDs == ref_material_IDs

        close(h5_file)
        rm("./c5g7_UO2_pin.xdmf")
        rm("./c5g7_UO2_pin.h5")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
This module contains classical rewriting functions and utilities
"""
module Classic

import ..Rule
import ..Theory
# import ..gettheory
import ..closure_generator
import ..@metatheory_init
using ..Util
using Base.Meta

@metatheory_init

include("matchcore_compiler.jl")
include("rewrite.jl")
include("match.jl")

export rewrite
export @rewrite
export @esc_rewrite
export @compile_theory
export @matcher
export @rewriter

end


text:
using MoltenSaltThermodynamics
using Test
using DataFrames, CSV

mydata1 = CSV.File("cocl2_naclkcl.csv") |> DataFrame
n = 2.0
T = 1073.0
edata = mydata1[!, Symbol("800C")].* -1
xdata = mydata1[!, Symbol("CoCl2_NaClKCl")]
model = infinitedilution_extrapolate(xdata, edata, n, T)

@testset "MoltenSaltThermodynamics.jl" begin
  @test isapprox(model.standardpotential, -1.1259; atol=1e-4)
  @test isapprox(model.slope, 0.04623; atol=1e-5)
end


text:
module PyPlotRecipes
using PyPlot
using RecipesBase

function replace_attr!(r::RecipeData, src::Symbol, dst::Symbol, trf=identity)
    if src ∈ keys(r.plotattributes)
        r.plotattributes[dst] = trf(r.plotattributes[src])
        delete!(r.plotattributes, src)
    end
end

function rplot(v, args...; kwargs...)
    map(RecipesBase.apply_recipe(Dict{Symbol,Any}(), v, args...)) do r
        replace_attr!(r, :markershape, :marker,
                      s -> Dict(:circle => "o")[s])
        plot(r.args...; r.plotattributes..., kwargs...)
    end
end

export rplot

end # module


text:
# identity function
id(x) = x

# remove elemnt from Array
#remove!(x,e) = filter!(n -> n ≢ e, x)

# returns a list of n times the vec concatenated
function cycle(vec, n)
    return collect(Base.Iterators.flatten((Base.Iterators.repeated(vec, n))))
end

# calculate the nernst potential
function vkf(ke,ki)
    8.314 * 293.15/96485.332 * log(ke/ki) * 1000
end

function vkf(ki)
    8.314 * 293.15/96485.332 * log(15/ki+1) * 1000
end

# is used to get the get the intermediate stepts for variables used
# to calculate the changes in the main variables
function setupstats()
    return Dict([
        ("t", []),
        ("vk", []),
        ("gkatpbar", []),
        ("gca", []),
        ("gkca", []),
    ])
end

function savestats(i)
    global stats
    push!(stats["t"]    , i.t)
    push!(stats["vk"]    , i.p["vk"])
    push!(stats["gkatpbar"]    , i.p["gkatpbar"])
end

# takes a dictionary and sorts it alphanumericaly by keys (returns a SortedDict)
function dsort(d::Dict)
    sorted_keys = map(x -> x[1], sort(pairs(d) |> collect, by=x->x[1]))
    sd = SortedDict()
    map(x -> sd[x] = d[x], sorted_keys)
    return sd
end

function get_keys_or_fields(v)
    try
        return keys(v)
    catch
        return fieldnames(typeof(v))
    end
end

# convert struct to named tuple
# https://discourse.julialang.org/t/get-fieldnames-and-values-of-struct-as-namedtuple/8991
gentup(struct_T) = NamedTuple{( fieldnames(struct_T)...,), Tuple{(fieldtype(struct_T,i) for i=1:fieldcount(struct_T))...}}
@generated function to_named_tuple(x)
    nt = Expr(:quote, gentup(x))
    tup = Expr(:tuple)
    for i=1:fieldcount(x)
        push!(tup.args, :(getfield(x, $i)) )
    end
    return :($nt($tup))
end

# calculates the difference between each element of a vector
function vecdiff(vec::Vector)
    map(v -> v[1] - v[2], zip(vec, [vec[2:end]; vec[end]]))
end

# enter the total time interval to get dy/dt value for the vector
function vecdiff(vec::Vector, interval::Number)
    vecdiff(vec) ./ (interval ./ length(vec))
end

# returns max/min of vector (there is also a builtin findmax/findmin function)
function max(vec::Vector)
    reduce(vec) do x,y
       x>y ? x : y end
end

function min(vec::Vector)
    reduce(vec) do x,y
       y>x ? x : y end
end

# returns a list of the values stored in a dict
function values(d::Dict; sorted=false)
    if sorted
        return sorted_values(d)
    end
    map(x -> x[2], collect(pairs(d)))
end

# returns a list of values in a dict sorted by key value
function sortedvalues(d::Dict)
    map(x -> x[2], sort(pairs(d) |> collect, by=x->x[1]))
end

function sortedvalues(d::Dict{Any, Dict})
    reduce(vcat, map(x -> sortedvalues(x), map(x -> x[2], sort(pairs(d) |> collect, by=x->x[1]))))
end

# By default the numerical integration returns an array for each timestep
# that contains the values for all variables at that time.
# For plotting and further calculations its easier to have a seperate
# time series for each variable.
function rowtocolumn(vec::AbstractVector{T}) where T <: AbstractVector
    dim1 = length(vec)
    dim2 = length(vec[1])
    vects = Vector{Vector}(undef, dim2)
    map(x -> vects[x] = Array{eltype(vec[1])}(undef, dim1), 1:dim2)
    @inbounds @fastmath for i in 1:dim1, j in 1:dim2
        vects[j][i] = vec[i][j]
    end
    return vects
end

# obsolete, but could be usefull for some calculations
# this converts a vector of vectors to a matrix.
function vector_to_matrix(vecvec::AbstractVector{T}) where T <: AbstractVector
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = Array{eltype(vecvec[1]), 2}(undef, dim1, dim2)
    @inbounds @fastmath for i in 1:dim1, j in 1:dim2
        my_array[i,j] = vecvec[i][j]
    end
    return my_array
end

function vec_to_matrix(vecvec::AbstractVector{T}) where T <: AbstractVector
    return vector_to_matrix(vecvec)
end

# scales the solutions to all be visible in the same plot
function scale_solution_columns(sol)
    #mat = vec_to_matrix(sol.u)
    mat = rowtocolumn(sol.u)
    mat[2] *= 100
    mat[3] *= 750
    mat[6] /= 7
    #if size(mat)[1] == 9
        #mat[9] *= 10
    #end
    return mat
end


text:
## --------------------------------------------------------------------------
_default_outfile(input_file) = replace(input_file, ".tex" => ".expanded.tex")

## --------------------------------------------------------------------------
_get_indent(input_line) = first(split(input_line, "\\input"))

## --------------------------------------------------------------------------
function _expand_inputs_parse_args(argsv::Vector=ARGS)

    # arg settings
    argset = ArgParse.ArgParseSettings()
    @ArgParse.add_arg_table! argset begin
    "--srcfile", "-f"
        help = "the path to the main file"
        arg_type = String
        default = load_setting("maintex", pwd())
    "--outfile", "-o"
        help = "the path to the main file"
        arg_type = String
        default = ""
    "--dry-run", "-d"
        help = "Only verbosity, no action"
        action = :store_true
    "--commented", "-c"
        help = "The extended file is all commented"
        action = :store_true
    "--verbose", "-v"
        help = "Print info while executing"
        action = :store_true
    end
    args = ArgParse.parse_args(argsv, argset)
    isempty(args["outfile"]) && (args["outfile"] = _default_outfile(args["srcfile"]))
    to_symbol_arg_dict(args)
end

## --------------------------------------------------------------------------
function expand_inputs(;
        srcfile::String="",
        outfile::String=_default_outfile(srcfile),
        dry_run::Bool=false,
        verbose=true,
        commented=true
    )

    # Set up src and dest
    srcdir = abspath(dirname(srcfile))
    outdir = abspath(dirname(outfile))
    dry_run || mkpath(outdir)
    !isfile(srcfile) && error(srcfile, " dont found!!!")

    verbose && @info("Args", srcfile, outfile, dry_run, verbose, commented)
    verbose && println()

    # expand infiles
    input_file_found = true
    lines = readlines(srcfile)
    while input_file_found
        input_file_found = false
        for li in eachindex(lines)
            line = lines[li]
            input_file = extract_input_file(line)
            isempty(input_file) && continue

            # if input file found, I'll inserted in place
            verbose && @info("Found input file", input_file)
            input_file_found = true

            input_file = joinpath(srcdir, input_file)
            !isfile(input_file) && error(input_file, " dont found. src line ", li, ": ", line)

            input_file_lines = readlines(input_file)
            for i in reverse(eachindex(input_file_lines))
                new_line = input_file_lines[i]

                isfirst = (i == firstindex(input_file_lines))
                isfirst ? lines[li] = new_line : insert!(lines, li + 1, new_line)
            end
            break
        end
    end
    verbose && println()

    # write to dest
    # dry_run || write(outfile, join(lines, "\n"))
    if !dry_run
        commented && (lines = map((line) -> string("% ", line), lines))
        open(outfile, "w") do io
            foreach((line) -> println(io, line), lines)
        end
    end
    verbose && @info("Output file", outfile)
    verbose && println()

    return lines

end


text:
using ClimateMachine.DGMethods: DGModel
using ClimateMachine.MPIStateArrays: MPIStateArray
using ClimateMachine.DGMethods.NumericalFluxes: NumericalFluxSecondOrder
using ClimateMachine.Mesh.Geometry: LocalGeometry
using ClimateMachine.Mesh.Grids: Direction

import ClimateMachine.BalanceLaws:
    flux_second_order!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

@inline function linearized_pressure(ρ, ρe, Φ)
    FT = eltype(ρ)
    γ = FT(gamma(param_set))
    if total_energy
        (γ - 1) * (ρe - ρ * Φ)
    else
        (γ - 1) * ρe
    end
end

abstract type DryAtmosLinearModel <: BalanceLaw end

function vars_state(lm::DryAtmosLinearModel, ::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
    end
end
vars_state(lm::DryAtmosLinearModel, st::Auxiliary, FT) =
    vars_state(lm.atmos, st, FT)

function update_auxiliary_state!(
    dg::DGModel,
    lm::DryAtmosLinearModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end
function flux_second_order!(
    lm::DryAtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
integral_load_auxiliary_state!(
    lm::DryAtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
integral_set_auxiliary_state!(lm::DryAtmosLinearModel, aux::Vars, integ::Vars) =
    nothing
reverse_integral_load_auxiliary_state!(
    lm::DryAtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
reverse_integral_set_auxiliary_state!(
    lm::DryAtmosLinearModel,
    aux::Vars,
    integ::Vars,
) = nothing
flux_second_order!(
    lm::DryAtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
) = nothing
function wavespeed(
    lm::DryAtmosLinearModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ref = aux.ref_state
    return soundspeed(ref.ρ, ref.p)
end

boundary_conditions(lm::DryAtmosLinearModel) = (1, 2)
function boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc,
    lm::DryAtmosLinearModel,
    args...,
)
    boundary_state!(nf, bc, lm.atmos, args...)
end
function boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc,
    lm::DryAtmosLinearModel,
    args...,
)
    nothing
end
init_state_auxiliary!(lm::DryAtmosLinearModel, aux::Vars, geom::LocalGeometry) =
    nothing
init_state_prognostic!(
    lm::DryAtmosLinearModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
) = nothing

struct DryAtmosAcousticGravityLinearModel{M} <: DryAtmosLinearModel
    atmos::M
    function DryAtmosAcousticGravityLinearModel(atmos::M) where {M}
        if atmos.ref_state === NoReferenceState()
            error("DryAtmosAcousticGravityLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end
function flux_first_order!(
    lm::DryAtmosAcousticGravityLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    ref = aux.ref_state

    flux.ρ = state.ρu
    pL = linearized_pressure(state.ρ, state.ρe, aux.Φ)
    flux.ρu += pL * I
    flux.ρe = ((ref.ρe + ref.p) / ref.ρ) * state.ρu
    nothing
end
function source!(
    lm::DryAtmosAcousticGravityLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    ::NTuple{1, Dir},
) where {Dir <: Direction}
    if Dir === VerticalDirection || Dir === EveryDirection
        ∇Φ = aux.∇Φ
        source.ρu -= state.ρ * ∇Φ
        if !total_energy
            source.ρe -= state.ρu' * ∇Φ
        end
    end
    nothing
end


text:
module Lif

export BaseNeuron, LIFNeuron, TSNeuron,
       show, reset!, step! #is_refractory, is_spiking

using ..LifConfig

import Base: show

abstract type BaseNeuron end
abstract type BaseParameters end
# ============================================================================ #
mutable struct LIFParameters <: BaseParameters
    # parameters in mV
    vrest::Float      #resting potential (mV)
    vreset::Float     #post-spike reset (mV)
    vspike::Float     #spike amplitude (mV)
    threshold::Float  #spike threshold (mV)

    # misc parameters
    rm::Float         #membrane resistance (Mohm)
    #cm::Float         #membrane capacitance (uF)
    tau::Float        #membrane time constant (ms)
    refractory_period::Float #refractory period (ms)
    xi::Float         #noise factor (std of Gaussian) (mV)

    # Ohm = kg * m^2 * s^-3 * A^-2
    # Farad = s^4 * A^2 * m^-2 * kg^-1

    function LIFParameters()
        return new(
            0.0,        #vrest
            -0.2,       #vreset
            1.22 * 6.0, #vspike
            1.22,       #threshold
            2.5,        #rm
            4.5,        #tau
            3.0,        #refractory_period
            0.2         #xi
        )
    end
end
# ============================================================================ #
mutable struct LIFNeuron <: BaseNeuron
    p::LIFParameters
    vm::Float
    refractory_start::Float
    function LIFNeuron()
        p = LIFParameters()
        return new(p, p.vrest, -Inf)
    end
end
# ---------------------------------------------------------------------------- #
function set!(x::LIFNeuron, f::Symbol, v::Float)
    if f in fieldnames(typeof(x))
        setfield!(x, f, v)
    elseif f in fieldnames(typeof(x.p))
        setfield!(x.p, f, v)
    else
        warn("Failed to set field '" * string(f) * "': field not found")
    end
end
# ---------------------------------------------------------------------------- #
function show(io::IO, x::LIFNeuron)
    print(io,x.vm)
end
# ---------------------------------------------------------------------------- #
function reset!(x::LIFNeuron)
    x.vm = x.p.vrest
    x.refractory_start = -Inf
end
# ---------------------------------------------------------------------------- #
function is_refractory(x::LIFNeuron, tnow::Number)
    return tnow <= (x.refractory_start + x.p.refractory_period)
end
# ---------------------------------------------------------------------------- #
function is_spiking(x::LIFNeuron)
    return x.vm > x.p.threshold
end
# ---------------------------------------------------------------------------- #
function step!(x::LIFNeuron, inp::Float, dt::Float, tnow::Float)
    spike = false
    if !is_refractory(x, tnow)

        # ######################
        # # NOISE ADDED HERE   #
        # ######################
        xi = randn()*x.p.xi / sqrt(dt)

        # NOTE: the correct equation is (Vrest - Vm + (Isyn * Rm)) / tau
        # see Vm equation derivation in quals2 proposal file vmeq_derivation.md
        # x.vm += (((x.p.vrest - x.vm + (inp * x.p.rm)) / x.p.tau) + xi)  * dt

        vinf = x.p.vrest + inp * x.p.rm
        x.vm = vinf + (x.vm - vinf) * exp(-dt / x.p.tau) + (xi*dt)

        if is_spiking(x)
            x.vm = x.p.vreset
            x.refractory_start = tnow
            spike = true
        end
    else
        x.vm = x.p.vrest #+ xi
    end

    return x.vm, spike
end
# ============================================================================ #
# abstract IndexGenerator
#
# # cycle through indicies 1:N
# type SerialIndex <: IndexGenerator
#     idx::UnitRange
#     ptr::Int
#     SerialIndex(n::Int) = new(1:n, 1)
# end
# # cycle through a permuted verison of 1:N, re-shuffling after reaching the end
# type ShuffledIndex <: IndexGenerator
#     idx::Vector{Int}
#     ptr::Int
#     ShuffledIndex(n::Int) = new(randperm(n), 1)
# end
# # pick a random index with replacement
# type RandomIndex <: IndexGenerator
#     idx::UnitRange
#     cur::Int
#     function RandomIndex(n::Int)
#         self = new(1:n, 0)
#         increment!(self)
#         return self
#     end
# end
# # ---------------------------------------------------------------------------- #
# reset!(x::ShuffledIndex) = shuffle!(x.idx)
# reset!(x::IndexGenerator) = nothing
# # ---------------------------------------------------------------------------- #
# function increment!(x::IndexGenerator)
#     if x.ptr == length(x.idx)
#         x.ptr = 1
#         reset!(x)
#     else
#         x.ptr += 1
#     end
# end
# function increment!(x::RandomIndex)
#     x.cur = rand(x.idx, 1)[1]
# end
# index(x::IndexGenerator) = x.idx[x.ptr]
# index(x::RandomIndex) = x.cur
# # ============================================================================ #
# type TSParameters <: BaseParameters
#     vspike::Float
#     TSParameters() = new(1.0)
# end
# type TSNeuron{T<:IndexGenerator} <: BaseNeuron
#     ts::Vector{Vector{Float}}
#     p::TSParameters
#     ptr::Int
#     idx::T
#     kt::Int
#     function TSNeuron{T<:IndexGenerator}(ts::Vector{Vector{Float}},
#         idx::T=SerialIndex(length(ts)))
#
#         #NOTE timestamps are in seconds but simulations run in ms
#         return new(ts * 1e3, TSParameters(), 1, idx, 0)
#     end
# end
# # ---------------------------------------------------------------------------- #
# function reset!(x::TSNeuron)
#     x.ptr = 1
#     x.kt = index(x.idx)
# end
# # ---------------------------------------------------------------------------- #
# function set!(x::TSNeuron, f::Symbol, v::Float)
#     #nothing to set
# end
# # ---------------------------------------------------------------------------- #
# function show(io::IO, x::TSNeuron)
#     print(io, @sprintf("TSNeuron(%d)", length(x.ts)))
# end
# # ---------------------------------------------------------------------------- #
# function step!(x::TSNeuron, inp, dt, tnow)
#     spike = false
#     if x.ptr > 0 && ((tnow + dt) > x.ts[x.kt][x.ptr])
#         spike = true
#         x.ptr = x.ptr < length(x.ts[x.kt]) ? x.ptr+1 : -1
#     end
#     return 0.0, spike
# end
# ============================================================================ #
end #END MODULE

# vinf = x.p.vrest + stim * x.p.rm
# x.vm = vinf + (x.vm - vinf) * exp(-dt / x.p.tau)

# x.vm += dt * ((x.p.vrest - x.vm) + (inp * x.p.rm) + xi) / x.p.tau


text:
export DynaAgent

"""
    DynaAgent(;kwargs...)

`DynaAgent` is first introduced in: *Sutton, Richard S. "Dyna, an integrated architecture for learning, planning, and reacting." ACM Sigart Bulletin 2.4 (1991): 160-163.*

# Keywords & Fields

- `π`::[`AbstractPolicy`](@ref): the policy to use
- `model`::[`AbstractEnvironmentModel`](@ref): describe the environment to interact with
- `buffer`::[`AbstractTurnBuffer`](@ref): used to store transitions between agent and environment
- `role=:DEFAULT`: used to distinguish different agents
- `plan_step::Int=10`: the count of planning steps

The main difference between [`DynaAgent`](@ref) and [`Agent`](@ref) is that an environment model is involved. It is best described in the book: *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*

![](../assets/img/RL_book_fig_8_1.png)

![](../assets/img/RL_book_fig_8_2.png)
"""
Base.@kwdef struct DynaAgent{P<:AbstractPolicy, B<:AbstractTurnBuffer, M<:AbstractEnvironmentModel, R} <: AbstractAgent
    π::P
    model::M
    buffer::B
    role::R = :DEFAULT
    plan_step::Int = 10
end

(agent::DynaAgent)(obs::Observation) = agent.π(obs)

function update!(agent::DynaAgent, experience::Pair)
    π, model, buffer = agent.π, agent.model, agent.buffer

    push!(buffer, experience)
    update!(model, buffer, π)  # model learning
    update!(π, buffer, model; plan_step = agent.plan_step)
end

update!(model::AbstractEnvironmentModel, buffer::AbstractTurnBuffer, π::AbstractPolicy) =
    update!(model, buffer)

function update!(model::AbstractEnvironmentModel, buffer::AbstractTurnBuffer)
    transitions = extract_transitions(buffer, model)
    if !isnothing(transitions)
        update!(model, transitions)
    end
end

function update!(
    model::PrioritizedSweepingSampleModel,
    buffer::AbstractTurnBuffer,
    π::AbstractPolicy,
)
    transition = extract_transitions(buffer, model)
    if !isnothing(transition)
        update!(model, transition, priority(transition, π))
    end
end

function update!(
    π::AbstractPolicy,
    buffer::AbstractTurnBuffer,
    model::AbstractEnvironmentModel;
    plan_step = 1,
)
    update!(π, buffer)  # direct RL
    update!(π, model; plan_step = plan_step)  # planning
end


text:
using LatinHypercubeSampling
using Random
using StableRNGs
using Test

@testset "AudzeEglais" begin

    LHC = [1 3; 3 1; 2 2]
    n = size(LHC,1)
    dist = zeros(Float64,Int(n*(n-1)*0.5))

    @test_logs (:warn,"AudzeEglaisObjective!(dist,LHC) is deprecated and does not differ from AudzeEglaisObjective(LHC)") AudzeEglaisObjective!(dist,LHC)
    @test AudzeEglaisObjective(LHC) ≈ 0.88888888888888888
end

@testset "mutateLHC" begin

    rng = StableRNGs.StableRNG(1)
    a = [0,0]
    LHC = randomLHC(rng,3,2)

    rng = StableRNGs.StableRNG(1)
    LatinHypercubeSampling.mutateLHC!(rng,LHC,a)
    @test LHC == [  2  1
                    3  2
                    1  3]
end

@testset "tournament!" begin

    pop = Array{Array{Float64}}(undef, 15)
    popfit = Array{Float64}(undef, 15)
    fitnessInds = Vector{Int64}(undef, 15)
    rng = StableRNGs.StableRNG(1)

    for i = 1:15
        LHC = randomLHC(rng,15,3)
        popfit[i] = AudzeEglaisObjective(LHC)
    end

    sortperm!(fitnessInds,popfit)

    tour_inds = Array{Int}(1:15)     #Storage of indices for tournament selection
    tour_fitinds = Array{Int}(1:15)  #Storage of fitness for tournament selection

    @test LatinHypercubeSampling.tournament!(rng,fitnessInds,15,tour_inds,1) == 1
    @test LatinHypercubeSampling.tournament!(rng,fitnessInds,15,tour_inds,0) == 6
end

@testset "fixedcross" begin

    parone = [1,2,3,4,5,6,7,8]
    partwo = [4,2,5,1,6,8,3,7]
    offsprone = similar(parone)
    offsprtwo = similar(parone)
    randlock = ReentrantLock()
    rng = StableRNGs.StableRNG(1)
    @test [1,2,3,4,5,6,8,7] == LatinHypercubeSampling._fixedcross!(rng,offsprone,parone,partwo,randlock)
    rng = StableRNGs.StableRNG(1)
    @test [4,2,5,1,3,6,7,8]== LatinHypercubeSampling._fixedcross!(rng,offsprone,partwo,parone,randlock)

    rng = StableRNGs.StableRNG(1)
    @test [1,2,3,4,5,6,8,7] == LatinHypercubeSampling.fixedcross!(rng,offsprone,offsprtwo,parone,partwo,randlock)[1]
    rng = StableRNGs.StableRNG(1)
    @test [4,2,5,1,3,6,7,8] == LatinHypercubeSampling.fixedcross!(rng,offsprone,offsprtwo,parone,partwo,randlock)[2]
end

@testset "inversion" begin

    individual = [1,2,3,4]
    rng = StableRNGs.StableRNG(1)
    LatinHypercubeSampling.inversion!(rng,individual)

    @test [1,2,4,3] == individual

end

@testset "is optimization result still an LHC" for threading in (true, false)
    numPoints = 64
    numDims = 10
    numGens = 100

    X = LHCoptim(numPoints,numDims,numGens;popsize=20,ntour=2,ptour=0.8,periodic_ae=false,threading=threading)[1]
    n, d = size(X)
    for i = 1:d
        @test length(unique(X[:,i])) == n
    end

    X = LHCoptim(numPoints,numDims,numGens;popsize=20,ntour=2,ptour=0.8,periodic_ae=true,threading=threading)[1]
    n, d = size(X)
    for i = 1:d
        @test length(unique(X[:,i])) == n
    end

    X = subLHCoptim(X,numPoints÷2,numGens;popsize=20,ntour=2,ptour=0.8,periodic_ae=false)[1]
    n, d = size(X)
    for i = 1:d
        @test length(unique(X[:,i])) == n
    end

    X = subLHCoptim(X,numPoints÷2,numGens;popsize=20,ntour=2,ptour=0.8,periodic_ae=true)[1]
    n, d = size(X)
    for i = 1:d
        @test length(unique(X[:,i])) == n
    end

end

@testset "is the fittest individual kept" for threading in (true,false)
    numPoints = 64
    numDims = 10
    numGens = 100

    X,fitnesses = LHCoptim(numPoints,numDims,numGens;popsize=20,ntour=2,ptour=0.8,threading=threading)
    @test sort(fitnesses) == fitnesses
    fitnesses = subLHCoptim(X,numPoints÷2,numGens;popsize=20,ntour=2,ptour=0.8)[2]
    @test sort(fitnesses) == fitnesses
end

@testset "categorical LHC" for threading in (true, false)
    numPoints = 64
    numCat = 4
    W = 1.0
    catW= 0.01
    dims = [Continuous(),Continuous(),Categorical(numCat,catW)]

    numGens = 100

    X = LHCoptim!(randomLHC(numPoints,dims),numGens;dims=dims,interSampleWeight=W,threading=threading)[1]

    n, d = size(X)
    for (i, dim) in enumerate(dims)
        if typeof(dim) == Continuous
            @test length(unique(X[:,i])) == n
        else
            @test length(unique(X[:,i])) == numCat
        end
    end

    @test Categorical(1) == Categorical(1, 0.0)

end

@testset "subLHCindex" begin
    X = randomLHC(5,2)
    Xsub = X[1:2:end,:]
    subInds = subLHCindex(X,Xsub)
    @test subInds == [1,3,5]
end

@testset "scaleLHC" begin
    plan = randomLHC(5,2)
    scaled_plan = scaleLHC(plan,[(-1,1),(10,100)])
    @test extrema(scaled_plan; dims=1) == [(-1.0, 1.0)  (10.0, 100.0)]
end

@testset "Global rng test" begin
    Random.seed!(1)
    before_lhc = rand(100)
    Random.seed!(1)

    LHCoptim(10,3,100,rng=Random.MersenneTwister(1))

    after_lhc = rand(100)

    @test before_lhc == after_lhc
end

max_stars_repo_name: MrUrq/LatinHypercubeSampling

max_issues_repo_name: MrUrq/LatinHypercubeSampling

max_forks_repo_name: MrUrq/LatinHypercubeSampling


text:
using Hiop
using LinearAlgebra

const INFINITY = 1e+20

function eval_f(x::Vector{Float64}, prob::HiopProblem)
    xrange = 1:prob.ns
    yrange = 2*prob.ns+1:2*prob.ns + prob.nd
    srange = prob.ns+1:2*prob.ns
    return obj = 0.5 * ( sum(x[xrange] .* (x[xrange] .- 1.))
                        + sum((prob.user_data.Q * x[yrange]) .* x[yrange])
    + sum(x[srange] .* x[srange])
   )
end

function eval_g(x::Vector{Float64}, g::Vector{Float64}, prob::HiopProblem)
    xrange = 1:prob.ns
    yrange = 2*prob.ns+1:2*prob.ns + prob.nd
    srange = prob.ns+1:2*prob.ns

    @assert 3+prob.ns == prob.m

    s = x[srange]
    y = x[yrange]
    for con_idx in 1:prob.m
        if con_idx - 1 < prob.ns
            g[con_idx] = x[con_idx] + s[con_idx]
        else
            @assert con_idx - 1 < prob.ns + 3
            if con_idx - 1 == prob.ns
                g[con_idx] = x[1]
                g[con_idx] += sum(x[srange])
                g[con_idx] += sum(x[yrange])
            elseif con_idx - 1 == prob.ns + 1
                g[con_idx] = x[2]
                g[con_idx] += sum(x[yrange])
            elseif con_idx - 1 == prob.ns + 2
                g[con_idx] = x[3]
                g[con_idx] += sum(x[yrange])
            else
                @assert false
            end
        end
    end
    g[1:prob.ns] += prob.user_data.Md * x[yrange]
    return Int32(1)
end

function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64}, prob::HiopProblem)
    xrange = 1:prob.ns
    yrange = 2*prob.ns+1:2*prob.ns + prob.nd
    srange = prob.ns+1:2*prob.ns

    grad_f[xrange] .= x[xrange] .- 0.5
    grad_f[yrange] .= prob.user_data.Q * x[yrange]
    grad_f[srange] = x[srange]
    return Int32(1)
end

function eval_jac_g(mode::Vector{Symbol}, x::Vector{Float64},
                    iJacS::Vector{Int32}, jJacS::Vector{Int32}, MJacS::Vector{Float64},
                    JacD::Vector{Float64}, prob::HiopProblem)
    @assert m == prob.ns + 3
    if :Structure in mode
        nnzit = 1
        for con_idx in 1:prob.ns
            # x
            iJacS[nnzit] = con_idx - 1
            jJacS[nnzit] = con_idx - 1
            nnzit += 1
            # y
            iJacS[nnzit] = con_idx - 1
            jJacS[nnzit] = con_idx + prob.ns - 1
            nnzit += 1
        end
        for con_idx in prob.ns+1:m
            if con_idx == prob.ns + 1
                # wrt x1
                iJacS[nnzit] = con_idx - 1
                jJacS[nnzit] = 0
                nnzit += 1
                # wrt s
                for i in 1:prob.ns
                    iJacS[nnzit] = con_idx - 1
                    jJacS[nnzit] = prob.ns + i - 1
                    nnzit += 1
                end
            else
                # wrt x2 or x3
                if con_idx-prob.ns - 1 == 1 || con_idx-prob.ns - 1 == 2
                    iJacS[nnzit] = con_idx - 1
                    jJacS[nnzit] = con_idx - prob.ns - 1
                    nnzit += 1
                else
                    @assert false
                end
            end
        end
        @assert nnzit == length(iJacS) + 1
    end
    # values for sparse Jacobian if requested by the solver
    if :Sparse in mode
        nnzit = 1
        for con_idx in 1:prob.ns
            # sparse Jacobian EQ w.r.t. x and s
            # x
            MJacS[nnzit] = 1.0
            nnzit += 1
            # s
            MJacS[nnzit] = 1.0
            nnzit += 1
        end
        for con_idx in prob.ns+1:m
            # parse Jacobian INEQ w.r.t x and s
            if con_idx - prob.ns - 1 == 0
                # wrt x1
                MJacS[nnzit] = 1.0
                nnzit += 1
                for i in 1:prob.ns
                    MJacS[nnzit] = 1.0
                    nnzit += 1
                end
            else
                # wrt x2 or x3
                if con_idx - prob.ns - 1 == 1 || con_idx - prob.ns - 1 == 2
                    MJacS[nnzit] = 1.0
                    nnzit += 1
                else
                    @assert false
                end
            end
        end
        @assert nnzit == length(iJacS) + 1
    end
    if :Dense in mode
        Mdvec = reshape(Md,(prob.ns * prob.nd,))
        JacD[1:prob.ns * prob.nd] .= Mdvec
        @assert prob.ns + 3 == prob.m
        JacDmat = reshape(JacD, (prob.nd, prob.m))
        for i in 1:3
            for j in 1:prob.nd
                JacDmat[j,prob.ns + i] = 1.0
            end
        end
    end
    return Int32(1)
end

function eval_h(mode::Vector{Symbol}, x::Vector{Float64}, obj_factor::Float64, lambda::Vector{Float64},
                iHSS::Vector{Int32}, jHSS::Vector{Int32}, MHSS::Vector{Float64}, HDD::Vector{Float64},
                iHSD::Vector{Int32}, jHSD::Vector{Int32}, MHSD::Vector{Float64}, prob)
    @assert length(MHSS) == 2 * prob.ns
    @assert length(MHSD) == 0
    if :Structure in mode
        for i in 1:2*prob.ns
            iHSS[i] = i - 1
            jHSS[i] = i - 1
        end
    end
    if :Sparse in mode
        for i in 1:2*prob.ns
            MHSS[i] = obj_factor
        end
    end
    if :Dense in mode
        Qvec = reshape(prob.user_data.Q, (prob.nd*prob.nd,))
        for i in 1:prob.nd^2
            HDD[i] = obj_factor * Qvec[i]
        end
    end
    return Int32(1)
end

struct User_data
    ns::Float64
    Q::Matrix{Float64}
    Md::Matrix{Float64}
    buf_y::Vector{Float64}
end

if length(ARGS) == 2
    ns = parse(Int, ARGS[1])
    nd = parse(Int, ARGS[2])
elseif @isdefined(ns) && @isdefined(nd)
    println("ns and nd predefined.")
elseif length(ARGS) == 0
    ns = 100
    nd = 100
else
    error("Wrong number of arguments.")
end

println("ns: $ns, nd: $nd")

Q = Matrix{Float64}(undef, nd, nd)
Q .= 1e-8
for i in 1:nd
    Q[i,i] += 2.0
end
for i in 2:nd-1
    Q[i,i+1] += 1.0
    Q[i+1,i] += 1.0
end

Md = Matrix{Float64}(undef, ns, nd)
Md .= -1.0
buf_y = Vector{Float64}(undef, nd)

user_data = User_data(ns, Q, Md, buf_y)

n = 2*ns + nd
m = ns+3
nx_sparse = 2*ns;
nx_dense = nd;
nnz_sparse_Jaceq = 2*ns;
nnz_sparse_Jacineq = 3+ns;
nnz_sparse_Hess_Lagr_SS = 2*ns;
nnz_sparse_Hess_Lagr_SD = 0;

x_L = Vector{Float64}(undef, n)
x_U = Vector{Float64}(undef, n)

g_L = Vector{Float64}(undef, m)
g_U = Vector{Float64}(undef, m)

# x
# for(int i=0; i<ns; ++i) xlow_[i] = -INFINITY;
for i in 1:ns
    x_L[i] = -INFINITY
end
# s
# for(int i=ns; i<2*ns; ++i) xlow_[i] = 0.;
for i in ns+1:2*ns
    x_L[i] = 0.0
end
# y
# xupp_[2*ns] = -4.;
# for(int i=2*ns+1; i<n; ++i) xlow_[i] = -INFINITY;
x_L[2*ns+1] = -4.0
for i in 2*ns+2:n
    x_L[i] = -INFINITY
end

# x
# for(int i=0; i<ns; ++i) xupp_[i] = 3.;
for i in 1:ns
    x_U[i] = 3.0
end
# s
# for(int i=ns; i<2*ns; ++i) xupp_[i] = +INFINITY;
for i in ns+1:2*ns
    x_U[i] = INFINITY
end
# y
# xupp_[2*ns] = 4.;
# for(int i=2*ns+1; i<n; ++i) xupp_[i] = +INFINITY;
x_U[2*ns+1] = 4.0
for i in 2*ns+2:n
    x_U[i] = INFINITY
end

g_L .= 0.0
g_U .= 0.0
g_L[end-2] = -2.0;    g_U[end-2] = 2.0
g_L[end-1] = -INFINITY; g_U[end-1] = 2.0
g_L[end]   = -2.0;    g_U[end]   = INFINITY

nlp = HiopProblem(ns, nd,
                  Int32(nx_sparse), Int32(nx_dense),
                  Int32(nnz_sparse_Jaceq), Int32(nnz_sparse_Jacineq),
                  Int32(nnz_sparse_Hess_Lagr_SS), Int32(nnz_sparse_Hess_Lagr_SD),
                  n, x_L, x_U,
                  m, g_L, g_U,
                  eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, user_data)
nlp.x0 .= 1.0

solveProblem(nlp)