text:
g = [
  [08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08]
  [49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00]
  [81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65]
  [52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91]
  [22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80]
  [24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50]
  [32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70]
  [67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21]
  [24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72]
  [21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95]
  [78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92]
  [16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57]
  [86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58]
  [19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40]
  [04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66]
  [88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69]
  [04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36]
  [20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16]
  [20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54]
  [01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48]
]

l = 4
m, n = size(g)
maxsofar = 0
for i = 1:m-(l-1)
  global maxsofar
  for j = 1:n-(l-1)
    p_col = prod(g[i:i+(l-1), j])
    p_row = prod(g[i,j:j+(l-1)])
    p_dia1 = g[i,j] * g[i+1,j+1] * g[i+2,j+2] * g[i+3,j+3]
    p_dia2 = g[i,j+3] * g[i+1,j+2] * g[i+2,j+1] * g[i+3,j]
    maxsofar = max(maxsofar, p_col)
    maxsofar = max(maxsofar, p_row)
    maxsofar = max(maxsofar, p_dia1)
    maxsofar = max(maxsofar, p_dia2)
  end
end
println(maxsofar)

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

# %%
tr = Data.lorenz96_nonstationary(;N=1e4, D = 4)
x,y,z = columns(tr)
figure()
plot3D(x, y, z; lw = 1.0, color = "C0", alpha = 0.2)
scatter3D(x, y, z; lw = 1.0, color = "C0", alpha = 0.2)
# compare this with actual attractor
tr = Data.lorenz96_chaotic(;N=1e4, D = 4)
x,y,z = columns(tr)
figure()
plot3D(x, y, z; lw = 1.0, color = "C0", alpha = 0.2)
scatter3D(x, y, z; color = "C0", alpha = 0.2)

# %%
datas = Symbol[]
labels = String[]

N = 1e5

push!(datas, :lorenz96_nonstationary)
push!(labels, "non-stationary Lorenz96")

push!(datas, :lorenz96_chaotic)
push!(labels, "standard Lorenz96")

qH = 1
qC = 2
# Calculate values for H, C
es, Hs, Cs = Vector{Float64}[], Vector{Float64}[], Vector{Float64}[]
for i in 1:length(datas)
    data = datas[i]
    params = @strdict N qH qC data
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = false,
        ignores = ["data"],
    )
    @unpack e, H, C = output
    push!(es, e); push!(Hs, H); push!(Cs, C)
end

# Do the actual plot
legendtitle = "inappropriate data"
fig, axs = mainplot(Hs, Cs, es, labels, legendtitle; qH, qC, tol = 0.25)

ext: jl

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using Documenter, QuDiffEq

const PAGES = [
    "Home" => "index.md",
    "Tutorial" => ["tutorial/lin.md","tutorial/nlin.md",] ,
    "Manual" => ["man/algs.md", "man/taylor.md", ]
]

makedocs(sitename="QuDiffEq.jl",
       modules = [QuDiffEq],
       format = Documenter.HTML(
             prettyurls = ("deploy" in ARGS),
             canonical = ("deploy" in ARGS) ? "https://quantumbfs.github.io/QuDiffEq.jl/latest/" : nothing,
             assets = ["assets/favicon.ico"],
             ),
         doctest = ("doctest=true" in ARGS),
         clean = false,
         linkcheck = !("skiplinks" in ARGS),
         pages = PAGES)

deploydocs(
    repo = "github.com/QuantumBFS/QuDiffEq.jl.git",
    target = "build",
    )


text:
# DESCRIPTION
# Simulates closed-loop dynamics using a dynamic inverse controller.
#
# INPUT
# - acfun: name of system dynamics function
# - finp: name of control input function
# - Tsim: total simulation time
# - dt: time step
# - x0: initial state vector
# - xdot0: initial state derivative vector
# - u0: initial control vector
#
# OUTPUT
# - state: state vector history
# - time: time vector history
#
# ------------------------------------------------------------------------------

function simulate_DI(acfun,finp,Tsim,dt,x0,xdot0,u0)
    # set up simulation
    # load DI parameters
    global k_lat = load("DI.jld","k_lat")
    global k_lon = load("DI.jld","k_lon")
    global k_ped = load("DI.jld","k_ped")
    global kp_roll = load("DI.jld","kp_roll")
    global ki_roll = load("DI.jld","ki_roll")
    global kp_pitch = load("DI.jld","kp_pitch")
    global ki_pitch = load("DI.jld","ki_pitch")
    global kp_yaw = load("DI.jld","kp_yaw")
    global ki_yaw = load("DI.jld","ki_yaw")
    global τ_p = load("DI.jld","tau_p")
    global τ_q = load("DI.jld","tau_q")
    global τ_r = load("DI.jld","tau_r")
    global τ_u = load("DI.jld","tau_u")
    global CATAB = load("DI.jld","CATAB")
    global CBinvTAB = load("DI.jld","CBinvTAB")
    # vector of times
    time=collect(0:dt:Tsim)
    # initialize state time history
    state=zeros(NSTATES+NDISTATES,length(time))
    dstate=zeros(NSTATES+NDISTATES,length(time))
    # set initial condition
    state[:,1]=x0_DI
    dstate[:,1]=xdot0_DI
    # simulate closed-loop system
    for i=1:length(time)-1
        # integrate
        sol=rk4(acfun,finp,time[i],dt,dstate[:,i],state[:,i],u0)
        state[:,i+1]=sol[1]
        dstate[:,i+1]=sol[2]
    end
    #
    return state, time
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
# To be used instead of _PMACDC.variable_converter_ne() - supports deduplication of variables
function variable_ne_converter_indicator(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            Z_dc_conv_ne = _PM.var(pm, nw)[:conv_ne] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne",
                binary = true,
                start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "branchdc_tnep_start", 1.0)
            )
        else
            Z_dc_conv_ne = _PM.var(pm, nw)[:conv_ne] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne",
                lower_bound = 0,
                upper_bound = 1,
                start = _PM.comp_start_value(_PM.ref(pm, nw, :convdc_ne, i), "branchdc_tnep_start", 1.0)
            )
        end
    else
        Z_dc_conv_ne = _PM.var(pm, nw)[:conv_ne] = _PM.var(pm, first_n)[:conv_ne]
    end
    report && _PM.sol_component_value(pm, nw, :convdc_ne, :isbuilt, _PM.ids(pm, nw, :convdc_ne), Z_dc_conv_ne)
end

function variable_ne_converter_investment(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            investment = _PM.var(pm, nw)[:conv_ne_investment] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne_investment",
                binary = true,
                start = 0
            )
        else
            investment = _PM.var(pm, nw)[:conv_ne_investment] = JuMP.@variable(pm.model,
                [i in _PM.ids(pm, nw, :convdc_ne)], base_name="$(nw)_conv_ne_investment",
                lower_bound = 0,
                upper_bound = 1,
                start = 0
            )
        end
    else
        investment = _PM.var(pm, nw)[:conv_ne_investment] = _PM.var(pm, first_n)[:conv_ne_investment]
    end
    report && _PM.sol_component_value(pm, nw, :convdc_ne, :investment, _PM.ids(pm, nw, :convdc_ne), investment)
end

# To be used instead of _PMACDC.variable_dc_converter_ne() - supports deduplication of variables
function variable_dc_converter_ne(pm::_PM.AbstractPowerModel; kwargs...)
    _PMACDC.variable_conv_tranformer_flow_ne(pm; kwargs...)
    _PMACDC.variable_conv_reactor_flow_ne(pm; kwargs...)
    variable_ne_converter_indicator(pm; kwargs..., relax=true) # FlexPlan version: replaces _PMACDC.variable_converter_ne().
    variable_ne_converter_investment(pm; kwargs...)

    _PMACDC.variable_converter_active_power_ne(pm; kwargs...)
    _PMACDC.variable_converter_reactive_power_ne(pm; kwargs...)
    _PMACDC.variable_acside_current_ne(pm; kwargs...)
    _PMACDC.variable_dcside_power_ne(pm; kwargs...)
    # _PMACDC.variable_converter_firing_angle_ne(pm; kwargs...)

    _PMACDC.variable_converter_filter_voltage_ne(pm; kwargs...)
    _PMACDC.variable_converter_internal_voltage_ne(pm; kwargs...)
    #
    _PMACDC.variable_converter_to_grid_active_power_ne(pm; kwargs...)
    _PMACDC.variable_converter_to_grid_reactive_power_ne(pm; kwargs...)
end


text:
export field_add, field_sub, field_neg, field_mul, field_sqr, field_inv, field_exp,
    curve_add, curve_dbl, curve_gen, curve_map, curve_mul, curve_neg, curve_order, curve_sub,
    md_sha256

Base.:(==)(a::BN, b::BN) = iszero(ccall((:bn_cmp, LIB), Cint, (Ref{BN}, Ref{BN}), a, b))
Base.:(>)(a::BN, b::BN) = isone(ccall((:bn_cmp, LIB), Cint, (Ref{BN}, Ref{BN}), a, b))
Base.:(<)(a::BN, b::BN) = b > a
Base.:(<=)(a::BN, b::BN) = !(a > b)
Base.:(>=)(a::BN, b::BN) = !(a < b)

Base.:(*)(a::BN, b) = bn_mul_comba!(BN(), a, BN(b))
Base.:(*)(a::BN, b::BN) = bn_mul_comba!(BN(), b, a)
Base.:(*)(a::Integer, b::BN) = b * a
Base.:(*)(a::BN, b::ULimbMax) = bn_mul_dig!(BN(), a, convert(Limb, b))
function Base.:(*)(a::BN, b::SLimbMax)
    negb = -b
    c = BN()
    if signbit(b)
        bn_mul_dig!(c, a, convert(Limb, negb))
        c.sign = xor(c.sign, 1)
    else # side-channel resistance by balancing the if else branches
        c.sign = xor(c.sign, 1) # side-channel resistance (no logical effect)
        bn_mul_dig!(c, a, convert(Limb, b))
    end
    return c
end

Base.:(+)(a::BN, b) = bn_add!(BN(), a, BN(b))
Base.:(+)(a::BN, b::BN) = bn_add!(BN(), a, b)
Base.:(+)(a::Integer, b::BN) = b + a
Base.:(+)(a::BN, b::ULimbMax) = bn_add_dig!(BN(), a, convert(Limb, b))
function Base.:(+)(a::BN, b::SLimbMax)
    negb = -b
    return signbit(b) ? bn_sub_dig!(BN(), a, convert(Limb, negb)) : bn_add_dig!(BN(), a, convert(Limb, b))
end

Base.:(-)(a::BN) = bn_neg!(BN(), a)
Base.:(-)(a::BN, b) = bn_sub!(BN(), a, BN(b))
Base.:(-)(a::BN, b::BN) = bn_sub!(BN(), a, b)
Base.:(-)(a::Integer, b::BN) = -b + a
Base.:(-)(a::BN, b::ULimbMax) = bn_sub_dig!(BN(), a, convert(Limb, b))
function Base.:(-)(a::BN, b::SLimbMax)
    negb = -b
    return signbit(b) ? bn_add_dig!(BN(), a, convert(Limb, negb)) : bn_sub_dig!(BN(), a, convert(Limb, b))
end

# function Base.invmod(a::BN, m::BN)
#     c, d = BN(), BN()
#     ccall((:bn_gcd_ext_lehme, LIB), Cvoid, (Ref{BN}, Ref{BN}, Ptr{BN}, Ref{BN}, Ref{BN}), c, d, C_NULL, a, m)
#     if !isone(c)
#         val = (BigInt(a), BigInt(m))
#         throw(DomainError(val, "Greatest common divisor is $(val[1])."))
#     end
#     return bn_mod_basic!(c, isone(a.sign) ? bn_sub!(d, m, d) : bn_add!(d, m, d), m)
# end

# TODO: this is slower than the BigInt variant, but faster then bn_gcd_ext_lehme above
Base.invmod(a::BN, m) = BN(invmod(BigInt(a), m))
Base.invmod(a, m::BN) = BN(invmod(a, BigInt(m)))
Base.invmod(a::BN, m::BN) = BN(invmod(BigInt(a), BigInt(m)))

Base.:(==)(a::FP, b::FP) = iszero(ccall((:fp_cmp, LIB), Cint, (Ref{FP}, Ref{FP}), a, b))
Base.:(==)(a::FP2, b::FP2) = iszero(ccall((:fp2_cmp, LIB), Cint, (Ref{FP2}, Ref{FP2}), a, b))
Base.:(==)(a::FP6, b::FP6) = iszero(ccall((:fp6_cmp, LIB), Cint, (Ref{FP6}, Ref{FP6}), a, b))
Base.:(==)(a::FP12, b::FP12) = iszero(ccall((:fp12_cmp, LIB), Cint, (Ref{FP12}, Ref{FP12}), a, b))
Base.:(==)(p::EP, q::EP) = iszero(ccall((:ep_cmp, LIB), Cint, (Ref{EP}, Ref{EP}), p, q))
Base.:(==)(p::EP2, q::EP2) = iszero(ccall((:ep2_cmp, LIB), Cint, (Ref{EP2}, Ref{EP2}), p, q))

Base.:(+)(a::FPX, b::FPX) = field_add(a, b)
Base.:(-)(a::FPX, b::FPX) = field_sub(a, b)
Base.:(-)(a::FPX) = field_neg(a)
Base.:(*)(a::FPX, b::FPX) = field_mul(a, b)
Base.:(*)(a::FPX, b::Integer) = field_mul(a, typeof(a)(b))
Base.:(*)(a::Integer, b::FPX) = field_mul(b, typeof(b)(a))
Base.:(÷)(a::FPX, b::FPX) = field_mul(a, field_inv(b))
Base.:(//)(a::FPX, b::FPX) = field_mul(a, field_inv(b))
Base.:(^)(a::FPX, b::BN) = field_exp(a, b)
Base.:(^)(a::FPX, b::Integer) = field_exp(a, BN(b))
Base.inv(a::FPX) = field_inv(a)

Base.:(+)(a::EPX, b::EPX) = curve_add(a, b)
Base.:(-)(a::EPX, b::EPX) = curve_sub(a, b)
Base.:(-)(a::EPX) = curve_neg(a)
Base.:(*)(a::EPX, b::BN) = curve_mul(a, b)
Base.:(*)(a::BN, b::EPX) = curve_mul(b, a)
Base.:(*)(a::EPX, b::Integer) = curve_mul(a, BN(b))
Base.:(*)(a::Integer, b::EPX) = curve_mul(b, BN(a))

Base.rand(::Type{BN}; bits=MD_LEN_BITS) = bn_rand!(BN(), bits)
Base.rand(::Type{FP}) = fp_rand!(FP(undef))
Base.rand(::Type{FP2}) = fp2_rand!(FP2(undef))
Base.rand(::Type{FP6}) = fp6_rand!(FP6(undef))
Base.rand(::Type{FP12}) = fp12_rand!(FP12(undef))
Base.rand(::Type{EP}) = ep_rand!(EP(undef))
Base.rand(::Type{EP2}) = ep2_rand!(EP2(undef))

# optimised iszero methods
Base.iszero(a::BN) = isone(ccall((:bn_is_zero, LIB), Cint, (Ref{BN},), a))
Base.iszero(a::FP) = isone(ccall((:fp_is_zero, LIB), Cint, (Ref{FP},), a))
Base.iszero(a::FP2) = isone(ccall((:fp2_is_zero, LIB), Cint, (Ref{FP2},), a))
Base.iszero(a::FP6) = isone(ccall((:fp6_is_zero, LIB), Cint, (Ref{FP6},), a))
Base.iszero(a::FP12) = isone(ccall((:fp12_is_zero, LIB), Cint, (Ref{FP12},), a))

Base.zero(::Type{BN}) = BN()
Base.zero(::Type{FP}) = FP()
Base.zero(::Type{FP2}) = FP2()
Base.zero(::Type{FP6}) = FP6()
Base.zero(::Type{FP12}) = FP12()
Base.zero(::Type{EP}) = EP()
Base.zero(::Type{EP2}) = EP2()
Base.zero(::BN) = BN()
Base.zero(::FP) = FP()
Base.zero(::FP2) = FP2()
Base.zero(::FP6) = FP6()
Base.zero(::FP12) = FP12()
Base.zero(::EP) = EP()
Base.zero(::EP2) = EP2()

Base.one(::Type{BN}) = BN(1)
Base.one(::Type{FP}) = FP(1)
Base.one(::Type{FP2}) = FP2(1)
Base.one(::Type{FP6}) = FP6(1)
Base.one(::Type{FP12}) = FP12(1)
Base.one(::BN) = BN(1)
Base.one(::FP) = FP(1)
Base.one(::FP2) = FP2(1)
Base.one(::FP6) = FP6(1)
Base.one(::FP12) = FP12(1)

Base.isvalid(a::EP) = isone(ccall((:ep_is_valid, LIB), Cint, (Ref{EP},), a))
Base.isvalid(a::EP2) = isone(ccall((:ep2_is_valid, LIB), Cint, (Ref{EP2},), a))

Base.isinf(a::EP) = isone(ccall((:ep_is_infty, LIB), Cint, (Ref{EP},), a))
Base.isinf(a::EP2) = isone(ccall((:ep_is_infty, LIB), Cint, (Ref{EP2},), a))

Base.string(bn::BN) = string(BigInt(bn))

Base.Vector{UInt8}(a::BN) = bn_write_bin!(Vector{UInt8}(undef, bn_size_bin(a)), a)
Base.Vector{UInt8}(a::EP; pack=true) = ep_write_bin!(Vector{UInt8}(undef, ep_size_bin(a, pack)), a, pack)
Base.Vector{UInt8}(a::EP2; pack=true) = ep2_write_bin!(Vector{UInt8}(undef, ep2_size_bin(a, pack)), a, pack)

function Base.BigInt(bn::BN)
    bigint = Base.GMP.MPZ.realloc2((bn.used * sizeof(Limb)) << 3)
    @assert bigint.alloc >= bn.used
    if isone(bn.sign)
        bigint.size = -bn.used
    else
        bigint.size = iszero(bn) ? 0 : bn.used
    end
    unsafe_copyto!(bigint.d, unsafe_bnptr(bn), bn.used)
    return bigint
end

# NOTE: can be hand optimised but its probably a poor RoI
#Base.BigInt(fp::FP) = BigInt(BN(fp))

# automatic promotion rules
#Base.promote_rule(::Type{BN}, ::Type{Int}) = BN
#Base.promote_rule(::Type{BN}, ::Type{BigInt}) = BigInt

Base.mod(a::BN, m::BN) = bn_mod_basic!(BN(), a, m)
Base.mod(a, m::BN) = bn_mod_basic!(BN(), BN(a), m)
Base.mod(a::BN, m) = bn_mod_basic!(BN(), a, BN(m))

field_add(a::FP, b::FP) = fp_add_basic!(FP(undef), a, b)
field_sub(a::FP, b::FP) = fp_sub_basic!(FP(undef), a, b)
field_neg(a::FP) = fp_neg_basic!(FP(undef), a)
field_mul(a::FP, b::FP) = fp_mul_comba!(FP(undef), a, b)
field_sqr(a::FP) = fp_sqr_comba!(FP(undef), a)
field_inv(a::FP) = fp_inv_lower!(FP(undef), a)
field_exp(a::FP, b::BN) = fp_exp_slide!(FP(undef), a, b)

field_add(a::FP2, b::FP2) = fp2_add_integ!(FP2(undef), a, b)
field_sub(a::FP2, b::FP2) = fp2_sub_integ!(FP2(undef), a, b)
field_neg(a::FP2) = fp2_neg!(FP2(undef), a)
field_mul(a::FP2, b::FP2) = fp2_mul_integ!(FP2(undef), a, b)
field_sqr(a::FP2) = fp2_sqr_integ!(FP2(undef), a)
field_inv(a::FP2) = fp2_inv!(FP2(undef), a)
field_exp(a::FP2, b::BN) = fp2_exp!(FP2(undef), a, b)

field_add(a::FP6, b::FP6) = fp6_add!(FP6(undef), a, b)
field_sub(a::FP6, b::FP6) = fp6_sub!(FP6(undef), a, b)
field_neg(a::FP6) = fp6_neg!(FP6(undef), a)
field_mul(a::FP6, b::FP6) = fp6_mul_lazyr!(FP6(undef), a, b)
field_sqr(a::FP6) = fp6_sqr_lazyr!(FP6(undef), a)
field_inv(a::FP6) = fp6_inv!(FP6(undef), a)
field_exp(a::FP6, b::BN) = fp6_exp!(FP6(undef), a, b)

field_add(a::FP12, b::FP12) = fp12_add!(FP12(undef), a, b)
field_sub(a::FP12, b::FP12) = fp12_sub!(FP12(undef), a, b)
field_neg(a::FP12) = fp12_neg!(FP12(undef), a)
field_mul(a::FP12, b::FP12) = fp12_mul_lazyr!(FP12(undef), a, b)
field_sqr(a::FP12) = fp12_sqr_lazyr!(FP12(undef), a)
field_inv(a::FP12) = fp12_inv!(FP12(undef), a)
field_exp(a::FP12, b::BN) = fp12_exp!(FP12(undef), a, b)

curve_add(a::EP, b::EP) = ep_add_basic!(EP(undef), a, b)
curve_dbl(a::EP) = ep_dbl_basic!(EP(undef), a)
curve_gen(::Type{EP}) = ep_curve_get_gen!(EP(undef))
curve_map(::Type{EP}, msg::Vector{UInt8}) = ep_map!(EP(undef), msg)
curve_mul(a::EP, b::BN) = ep_mul_lwnaf!(EP(undef), a, b)
curve_neg(a::EP) = ep_neg_basic!(EP(undef), a)
curve_order(::Type{EP}) = ep_curve_get_ord!(BN())
curve_sub(a::EP, b::EP) = ep_sub_basic!(EP(undef), a, b)

curve_add(a::EP2, b::EP2) = ep2_add_basic!(EP2(undef), a, b)
curve_dbl(a::EP2) = ep2_dbl_basic!(EP2(undef), a)
curve_gen(::Type{EP2}) = ep2_curve_get_gen!(EP2(undef))
curve_map(::Type{EP2}, msg::Vector{UInt8}) = ep2_map!(EP2(undef), msg)
curve_mul(a::EP2, b::BN) = ep2_mul_lwnaf!(EP2(undef), a, b)
curve_neg(a::EP2) = ep2_neg_basic!(EP2(undef), a)
curve_order(::Type{EP2}) = ep2_curve_get_ord!(BN())
curve_sub(a::EP2, b::EP2) = ep2_sub_basic!(EP2(undef), a, b)

md_sha256(msg::Vector{UInt8}) = md_map_sh256(msg)


text:
# YAS style !!!
# This is a WIP

struct YASStyle <: AbstractStyle end
@inline getstyle(s::YASStyle) = s

yasformat(s::AbstractString; kwargs...) = format_text(
    s;
    kwargs...,
    always_for_in = true,
    whitespace_ops_in_indices = true,
    whitespace_typedefs = false,
    style = YASStyle(),
)

nestable(::YASStyle, cst::CSTParser.EXPR) = false

function p_kw(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for a in cst
        add_node!(t, pretty(ys, a, s), s, join_lines = true)
    end
    t
end

function p_curly(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
        elseif CSTParser.is_comma(a) && i == length(cst) - 1
            # remove trailing comma
        elseif is_closer(n) || is_opener(n)
            add_node!(t, n, s, join_lines = true)
        else
            add_node!(t, n, s, join_lines = t.endline == n.startline)
        end
    end
    t
end
@inline p_braces(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_curly(ys, cst, s)

function p_call(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))

    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
        elseif CSTParser.is_comma(a) && i == length(cst) - 1
            # remove trailing comma
        elseif is_closer(n) || is_opener(n)
            add_node!(t, n, s, join_lines = true)
        elseif n.typ === CSTParser.Parameters
            join_lines = t.endline == n.startline
            add_node!(t, n, s, join_lines = join_lines)
        else
            join_lines = t.endline == n.startline
            if join_lines && !is_opener(t[end])
                add_node!(t, Whitespace(1), s)
            end
            add_node!(t, n, s, join_lines = join_lines)
        end
    end
    t
end
@inline p_ref(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)
@inline p_vect(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)
@inline p_tupleh(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)
@inline p_parameters(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)
@inline p_comprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)

function p_whereopcall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    add_node!(t, pretty(ys, cst[1], s), s)
    add_node!(t, Whitespace(1), s)
    add_node!(t, pretty(ys, cst[2], s), s, join_lines = true)
    add_node!(t, Whitespace(1), s)

    add_braces =
        !CSTParser.is_lbrace(cst[3]) && cst.parent.typ !== CSTParser.Curly &&
        cst[3].typ !== CSTParser.Curly && cst[3].typ !== CSTParser.BracesCat

    brace = FST(CSTParser.PUNCTUATION, t.endline, t.endline, "{")
    add_braces && add_node!(t, brace, s, join_lines = true)

    for (i, a) in enumerate(cst.args[3:end])
        n = a.typ === CSTParser.BinaryOpCall ? pretty(ys, a, s, nospace = true) :
            pretty(ys, a, s)
        if is_closer(n) || is_opener(n)
            add_node!(t, n, s, join_lines = true)
        elseif CSTParser.is_comma(a) && !is_punc(cst[i+3])
            add_node!(t, n, s, join_lines = true)
        elseif a.typ === CSTParser.BinaryOpCall
            add_node!(t, n, s, join_lines = t.endline == n.startline)
        else
            add_node!(t, n, s, join_lines = t.endline == n.startline)
        end
    end
    brace = FST(CSTParser.PUNCTUATION, t.endline, t.endline, "}")
    add_braces && add_node!(t, brace, s, join_lines = true)
    t
end

function p_generator(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = a.typ === CSTParser.BinaryOpCall ? pretty(ys, a, s, nonest = true) :
            pretty(ys, a, s)
        if a.typ === CSTParser.KEYWORD
            add_node!(t, Whitespace(1), s)
            add_node!(t, n, s, join_lines = true)
            # add_node!(t, Whitespace(1), s)
            if a.kind === Tokens.FOR
                for j = i+1:length(cst)
                    eq_to_in_normalization!(cst[j], s.opts.always_for_in)
                end
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
        else
            join_lines = t.endline == n.startline
            join_lines && add_node!(t, Whitespace(1), s)
            add_node!(t, n, s, join_lines = join_lines)
        end
    end
    t
end
@inline p_filter(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_generator(ys, cst, s)

#
# Nesting
#

function nest_if_over_margin!(style, fst::FST, s::State, i::Int)
    margin = s.line_offset
    idx = findnext(n -> n.typ === PLACEHOLDER, fst.nodes, i + 1)
    if idx === nothing
        margin += sum(length.(fst[i+1:end])) + fst.extra_margin
    else
        margin += sum(length.(fst[i+1:idx]))
    end
    if margin > s.margin || is_comment(fst[i+1])
        fst[i] = Newline(length = fst[i].len)
        s.line_offset = fst.indent
    else
        nest!(style, fst[i], s)
    end
end

function n_call!(ys::YASStyle, fst::FST, s::State)
    line_offset = s.line_offset
    fst.indent = line_offset + sum(length.(fst[1:2]))

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === TRAILINGSEMICOLON
            n.val = ""
            n.len = 0
            nest!(ys, n, s)
        elseif n.typ === CSTParser.Parameters || n.typ === CSTParser.Generator
            n.indent = fst.indent
            n.extra_margin = 1
            nest!(ys, n, s)
        else
            n.extra_margin = 1
            nest!(ys, n, s)
        end
    end
end
@inline n_curly!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)
@inline n_ref!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)

function n_tupleh!(ys::YASStyle, fst::FST, s::State)
    line_offset = s.line_offset
    fst.typ !== CSTParser.Parameters && (fst.indent = line_offset)
    length(fst.nodes) > 0 && is_opener(fst[1]) && (fst.indent += 1)

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === TRAILINGSEMICOLON
            n.val = ""
            n.len = 0
            nest!(ys, n, s)
        elseif n.typ === CSTParser.Generator
            n.indent = fst.indent
            n.extra_margin = 1
            nest!(ys, n, s)
        else
            n.extra_margin = 1
            nest!(ys, n, s)
        end
    end
end
@inline n_braces!(ys::YASStyle, fst::FST, s::State) = n_tupleh!(ys, fst, s)
@inline n_vect!(ys::YASStyle, fst::FST, s::State) = n_tupleh!(ys, fst, s)
@inline n_parameters!(ys::YASStyle, fst::FST, s::State) = n_tupleh!(ys, fst, s)
@inline n_invisbrackets!(ys::YASStyle, fst::FST, s::State) = n_tupleh!(ys, fst, s)
@inline n_comprehension!(ys::YASStyle, fst::FST, s::State) = n_tupleh!(ys, fst, s)

function n_import!(ys::YASStyle, fst::FST, s::State)
    fst.indent = s.line_offset + sum(length.(fst[1:2]))
    for (i, n) in enumerate(fst.nodes)
        if n.typ === PLACEHOLDER
            nest_if_over_margin!(ys, fst, s, i)
        elseif n.typ === NEWLINE
            s.line_offset = fst.indent
        else
            nest!(ys, n, s)
        end
    end
end
@inline n_export!(ys::YASStyle, fst::FST, s::State) = n_import!(ys, fst, s)
@inline n_using!(ys::YASStyle, fst::FST, s::State) = n_import!(ys, fst, s)

n_whereopcall!(ys::YASStyle, fst::FST, s::State) =
    nest!(ys, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
n_chainopcall!(ys::YASStyle, fst::FST, s::State) =
    n_block!(DefaultStyle(ys), fst, s, custom_indent = s.line_offset)


text:
module HTML5 # Bukdu

include("HTML5/Form.jl")

end # module Bukdu.HTML5


text:
defs(mod::LLVM.Module)  = filter(f -> !isdeclaration(f), collect(functions(mod)))
decls(mod::LLVM.Module) = filter(f ->  isdeclaration(f) && !LLVM.isintrinsic(f),
                                 collect(functions(mod)))

## lazy module loading

using UUIDs

struct LazyModule
    pkg::Base.PkgId
    LazyModule(name, uuid) = new(Base.PkgId(uuid, name))
end

function Base.getproperty(lazy_mod::LazyModule, sym::Symbol)
    pkg = getfield(lazy_mod, :pkg)
    mod = get(Base.loaded_modules, pkg, nothing)
    if mod === nothing
        error("This functionality requires the $(pkg.name) package, which should be installed and loaded first.")
    end
    getfield(mod, sym)
end


text:
module LinAlg

using Arpack: eigs

using SimpleTraits
using SparseArrays: SparseMatrixCSC
import SparseArrays: blockdiag, sparse
using LinearAlgebra: I, Symmetric, diagm, dot, eigen, eigvals, norm, rmul!, tril, triu
import LinearAlgebra: Diagonal, diag, issymmetric, mul!

using ..LightGraphs

import Base: convert, size, eltype, ndims, ==, *, length

export convert,
    SparseMatrix,
    GraphMatrix,
    Adjacency,
    adjacency,
    Laplacian,
    CombinatorialAdjacency,
    CombinatorialLaplacian,
    NormalizedAdjacency,
    NormalizedLaplacian,
    StochasticAdjacency,
    StochasticLaplacian,
    AveragingAdjacency,
    AveragingLaplacian,
    PunchedAdjacency,
    Noop,
    diag,
    degrees,
    symmetrize,
    prescalefactor,
    postscalefactor,
    perron,

    non_backtracking_matrix,
    Nonbacktracking,
    contract!,
    contract,

    adjacency_matrix,
    laplacian_matrix,
    incidence_matrix,
    adjacency_spectrum,
    laplacian_spectrum,
    coo_sparse,
    spectral_distance

include("./graphmatrices.jl")
include("./spectral.jl")
include("./nonbacktracking.jl")

end


text:
# Determine the local system's time zone
# Based upon Python's tzlocal https://pypi.python.org/pypi/tzlocal
using Mocking: Mocking, @mock

if Sys.iswindows()
    import TimeZones.WindowsTimeZoneIDs: WINDOWS_TRANSLATION
end

"""
    localzone() -> TimeZone

Returns a `TimeZone` object that is equivalent to the system's current time zone.
"""
function localzone()
    # Only allow creating a TimeZone using standard and legacy IANA time zone database
    # names. We allow the use of legacy names here as most operating systems still use the
    # legacy names.
    mask = Class(:STANDARD) | Class(:LEGACY)

    @static if Sys.isapple()
        name = @mock read(`systemsetup -gettimezone`, String)  # Appears to only work as root
        if startswith(name, "Time Zone: ")
            name = strip(replace(name, "Time Zone: " => ""))
        else
            # link will be something like /usr/share/zoneinfo/Europe/Warsaw
            name = @mock readlink("/etc/localtime")
            name = match(r"(?<=zoneinfo/).*$", name).match
        end
        return TimeZone(name, mask)
    elseif Sys.isunix()
        name = ""

        # Try getting the time zone from the "TZ" environment variable
        # http://linux.die.net/man/3/tzset
        if haskey(ENV, "TZ")
            name = ENV["TZ"]

            # When the TZ format starts with a colon this indicates that the time zone information
            # should be read from a file.
            if startswith(name, ':')
                name = name[2:end]
            else
                return parse_tz_format(name)
            end

            if startswith(name, '/')
                return @mock open(name) do f
                    read_tzfile(f, "local")
                end
            else
                # Relative name matches pre-compiled time zone name
                istimezone(name, mask) && return TimeZone(name, mask)

                # The system time zone directory used depends on the (g)libc version
                tzdirs = ["/usr/lib/zoneinfo", "/usr/share/zoneinfo"]
                haskey(ENV, "TZDIR") && pushfirst!(tzdirs, ENV["TZDIR"])

                for dir in tzdirs
                    filepath = joinpath(dir, name)
                    (@mock isfile(filepath)) || continue
                    return @mock open(filepath) do f
                        read_tzfile(f, name)
                    end
                end

                error("unable to locate tzfile: $name")
            end
        end

        # Look for distribution specific configuration files that contain the time zone name.

        filename = "/etc/timezone"
        if @mock isfile(filename)
            @mock open(filename) do file
                name = read(file, String)

                # Get rid of host definitions and comments:
                name = strip(replace(name, r"#.*" => ""))
                name = replace(name, ' ' => '_')
            end

            istimezone(name, mask) && return TimeZone(name, mask)
        end

        # CentOS has a ZONE setting in /etc/sysconfig/clock,
        # OpenSUSE has a TIMEZONE setting in /etc/sysconfig/clock and
        # Gentoo has a TIMEZONE setting in /etc/conf.d/clock

        zone_re = r"(TIME)?ZONE\s*=\s*\"(?<name>.*?)\""
        for filepath in ("/etc/sysconfig/clock", "/etc/conf.d/clock")
            (@mock isfile(filepath)) || continue
            @mock open(filepath) do file
                for line in readlines(file)
                    matched = match(zone_re, line)
                    if matched != nothing
                        name = matched["name"]
                        name = replace(name, ' ' => '_')
                        break
                    end
                end
            end

            istimezone(name, mask) && return TimeZone(name, mask)
        end

        # systemd distributions use symlinks that include the zone name,
        # see manpage of localtime(5) and timedatectl(1)
        link = "/etc/localtime"
        if @mock islink(link)
            filepath = @mock readlink(link)
            pos = findnext(isequal('/'), filepath, 1)

            while pos !== nothing
                name = SubString(filepath, pos + 1)
                istimezone(name, mask) && return TimeZone(name, mask)
                pos = findnext(isequal('/'), filepath, pos + 1)
            end
        end

        # No explicit setting existed. Use localtime
        for filepath in ("/etc/localtime", "/usr/local/etc/localtime")
            (@mock isfile(filepath)) || continue
            return @mock open(filepath) do f
                read_tzfile(f, "local")
            end
        end
    elseif Sys.iswindows()
        # Windows powershell should be available on Windows 7 and above
        win_name = strip(@mock read(`powershell -Command "[TimeZoneInfo]::Local.Id"`, String))

        if haskey(WINDOWS_TRANSLATION, win_name)
            return TimeZone(WINDOWS_TRANSLATION[win_name], mask)
        else
            error("unable to translate to POSIX time zone name from: \"$win_name\"")
        end
    end

    error("Failed to find local time zone")
end

# Using conditional expressions `(?(condition)yes-regexp)` as control flow to indicate that
# that a captured group is dependent on a previous group begin matched. This could also be
# accomplished using nested groups.
const TZ_REGEX = r"""
    ^
    (?<name>[a-zA-Z]{3,})?
    (?(name)(?<sign>[+-]))?
    (?(name)(?<hour>\d+))?
    (?(hour)\:(?<minute>\d+))?
    (?(minute)\:(?<second>\d+))?
    $
    """x

"""
    parse_tz_format(str) -> TimeZone

Parse the time zone format typically provided via the "TZ" environment variable. Details on
the format can be found under the man page for
[tzset](http://man7.org/linux/man-pages/man3/tzset.3.html).

Currently this function handles only the first format which is a fixed time zone without
daylight saving time.
"""
function parse_tz_format(str::AbstractString)
    m = match(TZ_REGEX, str)

    # Currently the only supported TZ format is reading time zone information from
    # a file.
    if m === nothing
        throw(ArgumentError("Unhandled TZ environment variable format: \"$str\""))
    end

    parse_digits(s) = s === nothing ? 0 : Base.parse(Int, s)

    name = m[:name] === nothing ? "UTC" : m[:name]

    # Note: positive indidates the local time zone is west of the Prime Meridian and
    # negative if it is east. This is the opposite of what FixedTimeZone expects.
    sign_val = m[:sign] == "-" ? 1 : -1

    # The tzset specification indicates that hours must be between 0 and 24 and minutes and
    # seconds 0 and 59. If values exceed these bounds they are clamped rather than treating
    # the entire format as invalid.
    hour = clamp(parse_digits(m[:hour]), 0, 24)
    minute = clamp(parse_digits(m[:minute]), 0, 59)
    second = clamp(parse_digits(m[:second]), 0, 59)

    offset = sign_val * (hour * 3600 + minute * 60 + second)

    FixedTimeZone(name, offset)
end


text:
module HartreeFock

export extract_tij_Uijlk, solve_scf

using LinearAlgebra
import PyCall: pyimport

"""
Extract tij and Uijlk from a FermionOperator object, representing a Hamiltonian.
The Hamiltonian must be number conserving, and is allowed to contain upto two-body operators.
"""
function extract_tij_Uijlk(ham)
    ofermion = pyimport("openfermion")
    FermionOperator = ofermion.ops.operators.FermionOperator
    count_qubits = ofermion.utils.operator_utils.count_qubits
    normal_ordered = ofermion.transforms.opconversions.normal_ordered

    ham = normal_ordered(ham)
    if ham.many_body_order() > 4
        throw(DomainError("Allowed only upto two-body operators!"))
    end
    if !ham.is_two_body_number_conserving()
        throw(DomainError("ham must be number conserving!"))
    end
    nflavors = count_qubits(ham)
    T = ham.zero().constant isa Real ? Float64 : Complex{Float64}
    tij = zeros(T, nflavors, nflavors)
    Uijlk = zeros(T, nflavors, nflavors, nflavors, nflavors)
    for (ops, constant) in ham.terms
        flavors = collect((o[1]+1 for o in ops))
        if length(ops) == 2
            tij[flavors[1], flavors[2]] += constant
        elseif length(ops) == 4
            Uijlk[flavors[1], flavors[2], flavors[3], flavors[4]] += 2*constant
        else
            throw(DomainError("Invalid many_body_order!"))
        end
    end
    tij, Uijlk
end

"""
Compute Tij
"""
function compute_Tij(tij, barUijlk, rhoij)
    Tij = copy(tij)
    N = size(tij)[1]
    for i in 1:N, j in 1:N
        for k in 1:N, l in 1:N
            Tij[i,j] += barUijlk[i,k,j,l] * rhoij[k,l]
        end
    end
    Tij
end

function compute_rhoij!(rhoij, evecs, nelec)
    rhoij .= 0.0
    N = size(rhoij)[1]
    for e in 1:nelec
        for i in 1:N, j in 1:N
           rhoij[i, j] += conj(evecs[i,e]) * evecs[j,e]
        end
    end
end

"""
Solve the Hartree-Fock equation self-consistently:
   hat{H} = sum_{ij=1}^N t_{ij} hat{c}^dagger_i hat{c}_j
       + frac{1}{2} sum_{ijkl=1}^N U_{ijlk}hat{c}^dagger_i hat{c}^dagger_j hat{c}_k hat{c}_l
"""
function solve_scf(tij::Array{T,2}, Uijlk::Array{T,4}, rhoij0::Array{T,2}, nelec::Integer, niter::Integer, mixing::Float64) where T
    # Antisymmetric Coulomb tensor
    barUijlk = Uijlk - permutedims(Uijlk, [1, 2, 4, 3])
    rhoij = copy(rhoij0)
    rhoij_out = copy(rhoij)
    for iter in 1:niter
        Tij = compute_Tij(tij, Uijlk, rhoij)
        e = eigen(Hermitian(Tij))
        compute_rhoij!(rhoij_out, e.vectors, nelec)
        @. rhoij = (1-mixing) * rhoij + mixing * rhoij_out
    end
    rhoij
end

end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
module C10k

include("server.jl")

end


text:
@testset "preprocessors" begin

    @testset "ResizeImage" begin
        obs = (state = ones(4, 4),)
        p = ResizeImage(2, 2)
        @test get_state(p(obs)) == ones(2, 2)
    end

    @testset "StackFrames" begin
        A = ones(2, 2)
        p = StackFrames(2, 2, 3)

        for i in 1:3
            obs = (state = A * i,)
            p(obs)
        end

        obs = (state = A * 4,)
        @test get_state(p(obs)) == reshape(repeat([2, 3, 4]; inner = 4), 2, 2, :)

    end

end


text:
module Tools

export alias, aliasRegion!, setupIndex!, getMarkup, mergeMcosts
export checkConvg, fastFill!, whsaleprice, collectΠ!
using TypeDef
using NLsolve, LoopVectorization

function fastFill!(
  A    ::Matrix{Float64},
  nRow ::Int64,
  nCol ::Int64)

  @avx for j=1:nCol
    for i=1:nRow
      A[i,j] = 0
    end
  end
end

function alias(
  dim::dimType,
  sim::simType,
  optD2::Bool=false)

  s      = sim.s
  num    = sim.num
  D      = sim.D
  ms     = sim.ms
  ∂s∂p   = sim.∂s∂p
  ∂S∂p   = sim.∂S∂p
  ∂D∂p   = sim.∂D∂p

  fill!(ms,0)

  if optD2==true
    ∂²S∂p² = sim.∂²S∂p²
    ∂²D∂p² = sim.∂²D∂p²
    @inbounds for j=1:dim.product
      D[j] = 0
      for i=1:dim.product
        ∂D∂p[i,j] = 0
        for k=1:dim.product
          ∂²D∂p²[k,i,j] = 0
        end
      end
    end
    return s,num,D,ms,∂s∂p,∂S∂p,∂²S∂p²,∂D∂p,∂²D∂p²
  else
    @inbounds for j=1:dim.product
      D[j] = 0
      for i=1:dim.product
        ∂D∂p[i,j] = 0
      end
    end
    return s,num,D,ms,∂s∂p,∂S∂p,∂D∂p
  end
end

function aliasRegion!(
  num   ::Vector{Float64},
  indx  ::indxType,
  r     ::Int64,
  dim   ::dimType,
  mat   ::matType,
  par   ::parType,
  sim   ::simType)

  fill!(num,0)
  # r   = g - minimum(indx.market) + 1 # region index
  g   = indx.market[r]
  αᵢ  = view(par.αᵢ,:,r)
  μ̃   = view(sim.μ̃,:,:,r)
  pop = mat.pop[r]
  # active = view(indx.active,:,r)
  active = findall(view(indx.active,:,r))

  return αᵢ,μ̃,pop,active
end

function setupIndex!(
  sim ::simType,
  indx::indxType,
  pdline::BitArray{1},
  dim ::dimType,
  mat ::matType,
  mask::maskType,
  t   ::Int64)

  indx.market .= dim.region*(t-1) .+ collect(1:dim.region)
  sim.μ̃       .= view(mat.μ̃,:,:,indx.market)
  sim.p       .= view(mat.p,:,t)

  anyActive = findmax(view(mat.active,:,indx.market), dims=2)[1]
  for j=1:dim.product
    indx.anyActive[j] = anyActive[j] & pdline[j]
    indx.mno[j]       = mask.mno[j]  & indx.anyActive[j]
    indx.mvno[j]      = mask.mvno[j] & indx.anyActive[j]
    for r=1:dim.region
      g = indx.market[r]
      indx.active[j,r] = mat.active[j,g] & pdline[j]
    end
  end
  for j=1:3
    for i=1:dim.mvno
      indx.mvnp[i,j] = mat.mvnp[i,j] & indx.anyActive[dim.mno+i]
    end
  end

end

function getMarkup(
  p     ::Vector{Float64},
  wprice::Vector{Float64},
  mcost ::Vector{Vector{Float64}})

  markupR = p - mcost[2]
  markupW = wprice - mcost[1]
  markupV = p - wprice - mcost[2]

  return markupR, markupW, markupV
end

function mergeMcosts(
  mcostinc,
  mcostfree,
  dim::dimType,
  indx::indxType,
  mask::maskType)

  # mcost = zeros(dim.product)
  # mcost[indx.anyActive .& .!mask.free] = mcostinc
  # mcost[indx.anyActive .& mask.free] .= mcostfree

  mcost = [zeros(dim.product) for t=1:2]
  idx = findall(indx.anyActive .& .!mask.free)
  for j in eachindex(idx)
    i = idx[j]
    if indx.mno[i]
      mcost[2][i] = mcostinc[j]
    elseif indx.mvno[i]
      mcost[2][i] = 0
      mcost[1][i] = mcostinc[j]
    end
  end
  mcost[2][mask.free] = mcostfree
  return mcost
end

function checkConvg(
  iter::Int64,
  Op  ::NLsolve.SolverResults,
  fail::Bool=false)

  exitflag = (converged(Op)==true) & (maximum(Op.zero)<100) & (fail==false)

  return exitflag
end

function whsaleprice(
  p,
  dim::dimType,
  mat::matType,
  sim::simType,
  indx::indxType)

  any    = indx.anyActive
  D1own  = sim.∂D∂p' .* mat.own
  markup = -D1own[any,any] \ sim.D[any]
  wprice = zeros(dim.product)
  wprice[any] = p[any] - markup

  return wprice
end

function collectΠ!( # Collect profits both retail & wholesale
  Π,
  eq  ::equilType,
  vInt::Bool,
  t   ::Int64,
  dim ::dimType,
  mat ::matType,
  indx::indxType,
  mask::maskType)

  # Incumbents' profit
  prodIdx = findall(indx.anyActive)
  for j in prodIdx
    if vInt==false
      i_retail = mat.firm[j,2]
      Π[i_retail] += eq.profit[j,t,2]
      if mask.mvno[j]
        i_whsale = mat.firm[j,3]
        Π[i_whsale] += eq.profit[j,t,1]
      end
    else
      i_whsale = mat.firm[j,3]
      Π[i_whsale] += eq.profit[j,t,2]
    end
  end # end of j loop

end

end

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
function draw_dendrite_panel(gui_data::GuiData, app_data::Model.AppData, sock::Comm.SocClient)
    if CImGui.TreeNode("Dendrite")
        CImGui.PushItemWidth(80)

        # Row 1 *****************************************************
        gui_data.buffer = Model.prep_field(Model.tao_eff(app_data.model), 16)
        returned = CImGui.InputText("taoEff", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_tao_eff!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(250)

        gui_data.buffer = Model.prep_field(Model.dendrite_length(app_data.model), 16)
        returned = CImGui.InputText("Length", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_dendrite_length!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(450)

        gui_data.buffer = Model.prep_field(Model.dendrite_min_psp(app_data.model), 16)
        returned = CImGui.InputText("MinPSP", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_dendrite_min_psp!(app_data.model, gui_data.buffer)
        end

        CImGui.PopItemWidth()

        CImGui.TreePop()
    end
end


text:
"""
    KSDist(n)

Distribution of the (two-sided) Kolmogorov-Smirnoff statistic

```math
D_n = \\sup_x | \\hat{F}_n(x) -F(x)| \\sqrt(n)
```

``D_n`` converges a.s. to the Kolmogorov distribution.
"""
struct KSDist <: ContinuousUnivariateDistribution
    n::Int
end

@distr_support KSDist 1 / (2 * d.n) 1.0

#### Evaluation

# TODO: implement Simard and L'Ecuyer (2011) meta-algorithm
# requires Pomeranz and Pelz-Good algorithms
function cdf(d::KSDist,x::Float64)
    n = d.n
    b = x*n
    # known exact values
    if b <= 1/2
        return 0.0
    elseif b <= 1
        # accuracy could be improved
        return exp(lfact(n)+n*(log(2*b-1)-log(n)))
    elseif x >= 1
        return 1.0
    elseif b >= n-1
        return 1 - 2*(1-x)^n
    end

    a = b*x
    if a >= 18
        return 1.0
    elseif n <= 10_000
        if a <= 4
            return cdf_durbin(d,x)
        else
            return 1 - ccdf_miller(d,x)
        end
    else
        return cdf(Kolmogorov(),sqrt(a))
    end
end

function ccdf(d::KSDist,x::Float64)
    n = d.n
    b = x*n
    # Ruben and Gambino (1982) known exact values
    if b <= 0.5
        return 1.0
    elseif b <= 1
        return 1-exp(lfact(n)+n*(log(2*b-1)-log(n)))
    elseif x >= 1
        return 0.0
    elseif b >= n-1
        return 2*(1-x)^n
    end

    a = b*x
    if a >= 370
        return 0.0
    elseif a >= 4 || (n > 140 && a >= 2.2)
        return ccdf_miller(d,x)
    else
        return 1-cdf(d,x)
    end
end

# Durbin matrix CDF method, based on Marsaglia, Tsang and Wang (2003)
# modified to avoid need for exponent tracking
function cdf_durbin(d::KSDist,x::Float64)
    n = d.n
    k, ch, h = ceil_rems_mult(n,x)

    m = 2*k-1
    H = Matrix{Float64}(m, m)
    for i = 1:m, j = 1:m
        H[i,j] = i-j+1 >= 0 ? 1 : 0
    end
    r = 1.0
    for i = 1:m
        # (1-h^i) = (1-h)(1+h+...+h^(i-1))
        H[i,1] = H[m,m-i+1] = ch*r
        r += h^i
    end
    H[m,1] += h <= 0.5 ? -h^m : -h^m+(h-ch)
    for i = 1:m, j = 1:m
        for g = 1:max(i-j+1,0)
            H[i,j] /= g
        end
        # we can avoid keeping track of the exponent by dividing by e
        # (from Stirling's approximation)
        H[i,j] /= ℯ
    end
    Q = H^n
    s = Q[k,k]
    s*stirling(n)
end

# Miller (1956) approximation
function ccdf_miller(d::KSDist, x::Real)
    2*ccdf(KSOneSided(d.n),x)
end

## these functions are used in durbin and pomeranz algorithms
# calculate exact remainders the easy way
function floor_rems_mult(n,x)
    t = big(x)*big(n)
    fl = floor(t)
    lrem = t - fl
    urem = (fl+one(fl)) - t
    return convert(typeof(n),fl), convert(typeof(x),lrem), convert(typeof(x),urem)
end
function ceil_rems_mult(n,x)
    t = big(x)*big(n)
    cl = ceil(t)
    lrem = t - (cl - one(cl))
    urem = cl - t
    return convert(typeof(n),cl), convert(typeof(x),lrem), convert(typeof(x),urem)
end

# n!*(e/n)^n
function stirling(n)
    if n < 500
        s = 1.0
        for i = 1:n
            s *= i/n*ℯ
        end
        return s
    else
        # 3rd-order Stirling's approximation more accurate for large n
        twn = 12*n
        return sqrt(2*pi*n)*(1 + twn\(1 + (2*twn)\(1 - (15*twn)\139)))
    end
end


text:
function randa(T, K)
    data = valloc(T, 64, K)
    data .= rand(T, K)
    data
end


text:
function CellRef(n::AbstractString)
    @assert is_valid_cellname(n) "$n is not a valid CellRef."

    m_row = match(r"[0-9]+$", n)
    row_number = parse(Int, m_row.match)
    m_column_name = match(r"^[A-X]?[A-Z]?[A-Z]", n)

    return CellRef(n, row_number, decode_column_number(m_column_name.match))
end

CellRef(row::Int, col::Int) = CellRef(encode_column_number(col) * string(row))

"""
    decode_column_number(column_name::AbstractString) :: Int

Converts column name to a column number.

```julia
julia> XLSX.decode_column_number("D")
4
```

See also: `encode_column_number`.
"""
function decode_column_number(column_name::AbstractString) :: Int
    local result::Int = 0

    num_characters = length(column_name)

    iteration = 1
    for i in num_characters:-1:1
        column_char_as_int = Int(column_name[i])
        result += (26^(iteration-1)) * (column_char_as_int - 64) # From A to Z we have 26 values. 'A' Char is ASCII 65.
        iteration += 1
    end

    return result
end

"""
    encode_column_number(column_number::Int) :: String

Converts column number to a column name.

```julia
julia> XLSX.encode_column_number(4)
"D"
```

See also: `decode_column_number`.
"""
function encode_column_number(column_number::Int) :: String
    @assert column_number > 0 && column_number <= 16384 "Column number should be in the range from 1 to 16384."

    third_letter_sequence = div(column_number, 26^2)
    column_number = column_number - third_letter_sequence*(26^2)

    second_letter_sequence = div(column_number, 26) # 26^1
    column_number = column_number - second_letter_sequence*(26)

    first_letter_sequence = column_number # 26^0

    if third_letter_sequence > 0
        # result will have 3 letters
        return String([ Char(third_letter_sequence+64), Char(second_letter_sequence+64), Char(first_letter_sequence+64) ])

    elseif second_letter_sequence > 0
        # result will have 2 letters
        return String([ Char(second_letter_sequence+64), Char(first_letter_sequence+64) ])

    else
        # result will have 1 letter
        return String([ Char(first_letter_sequence+64) ])
    end
end

Base.string(c::CellRef) = c.name
Base.show(io::IO, c::CellRef) = print(io, string(c))

Base.:(==)(c1::CellRef, c2::CellRef) = c1.name == c2.name
Base.hash(c::CellRef) = hash(c.name)

const RGX_COLUMN_NAME = r"^[A-Z]?[A-Z]?[A-Z]$"
const RGX_CELLNAME = r"^[A-Z]+[0-9]+$"
const RGX_CELLRANGE = r"^[A-Z]+[0-9]+:[A-Z]+[0-9]+$"

function is_valid_column_name(n::AbstractString) :: Bool
    if !ismatch(RGX_COLUMN_NAME, n)
        return false
    end

    column_number = decode_column_number(n)
    if column_number < 1 || column_number > 16384
        return false
    end

    return true
end

# Cellname is bounded by A1 : XFD1048576
function is_valid_cellname(n::AbstractString) :: Bool

    if !ismatch(RGX_CELLNAME, n)
        return false
    end

    m_row = match(r"[0-9]+$", n)
    row = parse(Int, m_row.match)

    if row < 1 || row > 1048576
        return false
    end

    m_column = match(r"^[A-Z]+", n)
    if !is_valid_column_name(m_column.match)
        return false
    end

    return true
end

function is_valid_cellrange(n::AbstractString) :: Bool
    if !ismatch(RGX_CELLRANGE, n)
        return false
    end

    m_start = match(r"^[A-Z]+[0-9]+", n)
    start_name = m_start.match
    if !is_valid_cellname(start_name)
        return false
    end

    m_stop = match(r"[A-Z]+[0-9]+$", n)
    stop_name = m_stop.match
    if !is_valid_cellname(stop_name)
        return false
    end

    return true
end

macro ref_str(ref)
    CellRef(ref)
end

function CellRange(r::AbstractString)
    @assert ismatch(RGX_CELLRANGE, r) "Invalid cell range: $r."

    m_start = match(r"^[A-Z]+[0-9]+", r)
    start_name = CellRef(m_start.match)

    m_stop = match(r"[A-Z]+[0-9]+$", r)
    stop_name = CellRef(m_stop.match)

    return CellRange(start_name, stop_name)
end

Base.string(cr::CellRange) = "$(string(cr.start)):$(string(cr.stop))"
Base.show(io::IO, cr::CellRange) = print(io, string(cr))

Base.:(==)(cr1::CellRange, cr2::CellRange) = cr1.start == cr2.start && cr2.stop == cr2.stop
Base.hash(cr::CellRange) = hash(cr.start) + hash(cr.stop)

macro range_str(cellrange)
    CellRange(cellrange)
end

"""
    Base.in(ref::CellRef, rng::CellRange) :: Bool

Checks wether `ref` is a cell reference inside a range given by `rng`.
"""
function Base.in(ref::CellRef, rng::CellRange) :: Bool
    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    r = row_number(ref)

    if top <= r && r <= bottom
        left = column_number(rng.start)
        right = column_number(rng.stop)
        c = column_number(ref)

        if left <= c && c <= right
            return true
        end
    end

    return false
end

"""
    Base.issubset(subrng::CellRange, rng::CellRange)

Checks wether `subrng` is a cell range contained in `rng`.
"""
Base.issubset(subrng::CellRange, rng::CellRange) :: Bool = in(subrng.start, rng) && in(subrng.stop, rng)

function Base.size(rng::CellRange)
    top = row_number(rng.start)
    bottom = row_number(rng.stop)
    left = column_number(rng.start)
    right = column_number(rng.stop)

    return ( bottom - top + 1, right - left + 1 )
end

"""
    row_number(c::CellRef) :: Int

Returns the row number of a given cell reference.
"""
row_number(c::CellRef) :: Int = c.row_number

"""
    column_number(c::CellRef) :: Int

Returns the column number of a given cell reference.
"""
column_number(c::CellRef) :: Int = c.column_number

column_name(c::CellRef) :: String = encode_column_number(column_number(c))

"""
Returns (row, column) representing a `ref` position relative to `rng`.

For example, for a range "B2:D4", we have:

* "C3" relative position is (2, 2)

* "B2" relative position is (1, 1)

* "C4" relative position is (3, 2)

* "D4" relative position is (3, 3)

"""
function relative_cell_position(ref::CellRef, rng::CellRange)
    @assert ref ∈ rng "$ref is outside range $rng."

    top = row_number(rng.start)
    left = column_number(rng.start)

    r, c = row_number(ref), column_number(ref)

    return ( r - top + 1 , c - left + 1 )
end

#
# ColumnRange
#

const RGX_COLUMN_RANGE = r"^[A-Z]?[A-Z]?[A-Z]:[A-Z]?[A-Z]?[A-Z]$"

function is_valid_column_range(r::AbstractString) :: Bool
    if !ismatch(RGX_COLUMN_RANGE, r)
        return false
    end

    start_name = match(r"^[A-Z]+", r).match
    stop_name = match(r"[A-Z]+$", r).match

    if !is_valid_column_name(start_name)
        return false
    end

    if !is_valid_column_name(stop_name)
        return false
    end

    return true
end

function ColumnRange(r::AbstractString)
    @assert is_valid_column_range(r) "Invalid column range: $r."

    start_name = match(r"^[A-Z]+", r).match
    stop_name = match(r"[A-Z]+$", r).match

    return ColumnRange(decode_column_number(start_name), decode_column_number(stop_name))
end

convert(::Type{ColumnRange}, str::AbstractString) = ColumnRange(str)
convert(::Type{ColumnRange}, column_range::ColumnRange) = column_range

column_bounds(r::ColumnRange) = (r.start, r.stop)
Base.length(r::ColumnRange) = r.stop - r.start + 1

# ColumnRange iterator
Base.start(itr::ColumnRange) = itr.start
Base.done(itr::ColumnRange, column_index::Int) = column_index > itr.stop
Base.next(itr::ColumnRange, column_index::Int) = (encode_column_number(column_index), column_index + 1)

# CellRange iterator
struct CellRefIteratorState
    row::Int
    col::Int
end

Base.start(rng::CellRange) = CellRefIteratorState(row_number(rng.start), column_number(rng.start))
Base.done(rng::CellRange, state::CellRefIteratorState) = state.row > row_number(rng.stop)

function Base.length(rng::CellRange)
    (r, c) = size(rng)
    return r * c
end

# (i, state) = next(I, state)
function Base.next(rng::CellRange, state::CellRefIteratorState)
    local next_state::CellRefIteratorState
    if state.col == column_number(rng.stop)
        # reached last column. Go to the next row.
        next_state = CellRefIteratorState(state.row + 1, column_number(rng.start))
    else
        # go to the next column
        next_state = CellRefIteratorState(state.row, state.col + 1)
    end

    return CellRef(state.row, state.col), next_state
end


text:
using NonUniformRandomVariateGeneration
using Distributions
using BenchmarkTools

function benchNative(μ::Float64, N::Int64)
  m::Float64 = 0.0
  for i in 1:N
    m += samplePoisson(μ)
  end
  return m / N
end

function benchDistributions(μ::Float64, N::Int64)
  pois = Distributions.Poisson(μ)
  m::Float64 = 0.0
  for i in 1:N
    m += rand(pois)
  end
  return m / N
end

@btime benchNative(5.0, 1024*1024*8)
@btime benchDistributions(5.0, 1024*1024*8)

@btime benchNative(10.0, 1024*1024*8)
@btime benchDistributions(10.0, 1024*1024*8)

@btime benchNative(20.0, 1024*1024*8)
@btime benchDistributions(20.0, 1024*1024*8)


text:
using PowerModels
using Memento

# Suppress warnings during testing.
setlevel!(getlogger(PowerModels), "error")

using Ipopt
using Pajarito
using GLPKMathProgInterface
using SCS

# needed for Non-convex OTS tests
if (Pkg.installed("AmplNLWriter") != nothing && Pkg.installed("CoinOptServices") != nothing)
    using AmplNLWriter
    using CoinOptServices
end

using Base.Test

# default setup for solvers
ipopt_solver = IpoptSolver(tol=1e-6, print_level=0)
pajarito_solver = PajaritoSolver(mip_solver=GLPKSolverMIP(), cont_solver=ipopt_solver, log_level=0)
scs_solver = SCSSolver(max_iters=1000000, verbose=0)

include("common.jl")

include("pti.jl")

include("matpower.jl")

include("output.jl")

include("modify.jl")

include("data.jl")

include("pf.jl")

include("opf.jl")

include("ots.jl")

include("tnep.jl")

include("multinetwork.jl")

include("docs.jl")


text:
using Documenter
using Yao
using Yao: YaoBlocks, YaoArrayRegister, YaoBase, YaoSym
using YaoBase: BitBasis
using YaoBlocks: AD
using YaoBlocks: Optimise
using Documenter.Writers.HTMLWriter
using Documenter.Utilities.DOM
using Documenter.Utilities.DOM: Tag, @tags
# Evil Prirate

const base_url = raw"https://yaoquantum.org"

const top_nav = """

    <button class="navbar-toggler" type="button" data-toggle="collapse" data-target="#navbarSupportedContent" aria-controls="navbarSupportedContent" aria-expanded="false" aria-label="Toggle navigation">

    </button>

          Tutorial

            Documentation(current)

            Benchmark

          SoC

          Research

            GitHub

            Contribute

"""

function HTMLWriter.render_html(
    ctx,
    navnode,
    head,
    sidebar,
    navbar,
    article,
    footer,
    scripts::Vector{DOM.Node} = DOM.Node[],
)
    @tags html body div
    DOM.HTMLDocument(html[:lang=>"en"](
        head,
        body(
            Tag(Symbol("#RAW#"))(top_nav),
            div[".documenter-wrapper#documenter"](
                sidebar,
                div[".docs-main"](navbar, article, footer),
                HTMLWriter.render_settings(ctx),
            ),
        ),
        scripts...,
    ))
end

const PAGES = [
    "Home" => "index.md",
    "Manual" => Any[
        "man/array_registers.md",
        "man/symbolic.md",
        "man/blocks.md",
        "man/automatic_differentiation.md",
        "man/simplification.md",
        "man/base.md",
        "man/registers.md",
        "man/bitbasis.md",
        "man/extending_blocks.md",
    ],
    "Benchmark" => "benchmarks.md",
    "Developer Notes" => "dev/index.md",
]

makedocs(
    modules = [Yao, YaoBase, YaoArrayRegister, YaoBlocks, BitBasis, YaoSym, AD, Optimise],
    format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = ("deploy" in ARGS) ? "https://docs.yaoquantum.org/" : nothing,
        assets = [
            "assets/main.css",
            asset("https://yaoquantum.org/assets/main.css"),
            asset("http://yaoquantum.org/favicon.ico"),
            asset("https://stackpath.bootstrapcdn.com/bootstrap/4.4.1/css/bootstrap.min.css"),
        ],
    ),
    doctest = ("doctest=true" in ARGS),
    clean = false,
    sitename = "Documentation | Yao",
    linkcheck = !("skiplinks" in ARGS),
    pages = PAGES,
)

deploydocs(repo = "github.com/QuantumBFS/Yao.jl.git", target = "build")


text:
include("dynamic_symtridiagonal.jl")
include("dynamic_lanczos.jl")
using LinearAlgebra, Plots

n = 32
m = 10
A = Diagonal(LinRange(1e-3, 1, n).^3)
il = DIL(A, ones(n), n, true)
ChangeVelocityTol!(il.dynamic_symtridiagonal, 0.0)
fig1 = plot(legend=false)
for II in m:-1:1
    il()
    Ritzs = GetRitzValues(il)
    ConvergedRitzs = GetConverged(il.dynamic_symtridiagonal)
    scatter!(fig1, Ritzs, II*(Ritzs|>length|>ones), marker=:x)
    scatter!(
        fig1,
        ConvergedRitzs,
        II*(ConvergedRitzs|>length|>ones),
        marker=:square
    )
end
display(fig1)


text:
# ---
# title: 227. Basic Calculator II
# id: problem227
# author: Qling
# date: 2021-03-11
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/basic-calculator-ii/description/>
# hidden: true
# ---
#
# Given a string `s` which represents an expression, _evaluate this expression
# and return its value_.
#
# The integer division should truncate toward zero.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "3+2*2"
#     Output: 7
#
#
# **Example 2:**
#
#
#
#     Input: s = " 3/2 "
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: s = " 3+5 / 2 "
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 3 * 105`
#   * `s` consists of integers and operators `('+', '-', '*', '/')` separated by some number of spaces.
#   * `s` represents **a valid expression**.
#   * All the integers in the expression are non-negative integers in the range `[0, 231 - 1]`.
#   * The answer is **guaranteed** to fit in a **32-bit integer**.
#
#
## @lc code=start
using LeetCode

function calculate_ii(s::String)::Int
    num, stk, sign = 0, Int64[], '+'

    for (i, c) in enumerate(s)
        isdigit(c) && (num = 10 * num + parse(Int, c))
        if c in ['-', '+', '/', '*'] || i == length(s)
            if sign == '+'
                push!(stk, num)
            elseif sign == '-'
                push!(stk, -num)
            elseif sign == '*'
                push!(stk, pop!(stk) * num)
            elseif sign == '/'
                push!(stk, pop!(stk) ÷ num)
            end

            num = 0
            sign = c
        end
    end

    return sum(stk)
end

## @lc code=end


text:
import Random: shuffle

@generated function myshuffle(xs)
    canShuffleEmpty = try
        shuffle([])
        true
    catch
        false
    end
    if canShuffleEmpty
        :(shuffle(xs))
    else
        :(length(xs)>0 ? shuffle(xs) : xs)
    end
end

@generated function myshuffle(rng, xs)
    canShuffleEmpty = try
        shuffle([])
        true
    catch
        false
    end
    if canShuffleEmpty
        :(shuffle(rng, xs))
    else
        :(length(xs)>0 ? shuffle(rng, xs) : xs)
    end
end

@doc """
    loop_update!(model, param::Parameter)
    loop_update!(model, T::Real,
                 Jz::AbstractArray,
                 Jxy::AbstractArray,
                 Gamma:AbstractArray)

Updates spin configuration by loop algorithm
under the temperature `T = param["T"]` and coupling constants `Jz, Jxy` and transverse field `Gamma`
"""
@inline function loop_update!(model::QuantumXXZ, param::Parameter)
    p = convert_parameter(model, param)
    return loop_update!(model, p...)
end

function loop_update!(model::QuantumXXZ, T::Real,
                      Jzs::AbstractArray, Jxys::AbstractArray, Gs::AbstractArray)
    rng = model.rng
    lo_types = [LET_FMLink, LET_AFLink, LET_Vertex, LET_Cross]
    nsites = numsites(model)
    S2 = model.S2
    nspins = nsites*S2
    nbonds = numbonds(model)
    nst = numsitetypes(model)
    nbt = numbondtypes(model)
    weights = zeros(4*nbt+nst)
    for i in 1:nbt
        nb = numbonds(model,i)*S2*S2
        z = nb*Jzs[i]
        x = nb*abs(Jxys[i])
        if z > x
            ## AntiFerroIsing like
            weights[(4i-3):(4i)] .= 0.5*[0.0, z-x, x, 0.0]
        elseif z < -x
            ## FerroIsing like
            weights[(4i-3):(4i)] .= 0.5*[-z-x, 0.0, 0.0, x]
        else
            ## XY like
            weights[(4i-3):(4i)] .= 0.25*[0.0, 0.0, x+z, x-z]
        end
    end
    for i in 1:nst
        ns = numsites(model,i)*S2
        weights[4nbt+i] = 0.5*ns*Gs[i]
    end
    accumulated_weights = cumsum(weights)
    op_dt = T/accumulated_weights[end]

    spins = model.spins[:]
    currents = collect(1:nspins)
    uf = UnionFind(nspins)

    ops = LocalLoopOperator[]
    iops = 1 # index of the next operator in the original string
    t = randexp(rng)*op_dt
    while t <= 1.0 || iops <= length(model.ops)
        if iops > length(model.ops) || t < model.ops[iops].time
            ## INSERT
            ot = searchsortedfirst(accumulated_weights, rand(rng)*accumulated_weights[end])
            if ot <= 4nbt
                ## Bond
                bt = ceil(Int, ot/4)
                ot = mod1(ot,4)
                lo_type = lo_types[ot]
                b = bonds(model,bt)[rand(rng, 1:numbonds(model,bt))]
                s1,s2 = source(b), target(b)
                ss1, ss2 = rand(rng, 1:S2, 2)
                if ifelse(spins[site2subspin(s1,ss1,S2)] == spins[site2subspin(s2,ss2,S2)],
                          lo_type == LET_FMLink || lo_type == LET_Cross,
                          lo_type == LET_AFLink || lo_type == LET_Vertex)
                    push!(ops, LocalLoopOperator(lo_type, t, nsites+b.id, (ss1, ss2)))
                    t += randexp(rng)*op_dt
                else
                    t += randexp(rng)*op_dt
                    continue
                end
            else
                ## Site
                st = ot-4nbt
                s = sites(model,st)[rand(rng, 1:numsites(model,st))]
                ss = rand(rng, 1:S2)
                push!(ops, LocalLoopOperator(LET_Cut, t, s, (ss,ss)))
                t += randexp(rng)*op_dt
            end
        else
            op = model.ops[iops]
            iops += 1
            if op.isdiagonal
                ## REMOVE
                continue
            elseif op.let_type == LET_Cut
                push!(ops, op)
            else
                push!(ops, op)
                op = ops[end]
                b = op.space-nsites
                ss1,ss2 = op.subspace
                bt = bondtype(model,b)
                otype = ifelse(rand(rng)*(weights[4bt-1]+weights[4bt])<weights[4bt-1], LET_Vertex, LET_Cross)
                op.let_type = otype
            end
        end

        op = ops[end]
        if op.let_type == LET_Cut
            s = op.space
            ss = op.subspace[1]
            subspin = site2subspin(s,ss,S2)
            op.bottom_id = currents[subspin]
            c = addnode!(uf)
            currents[subspin] = c
            op.top_id = c
            spins[subspin] *= ifelse(op.isdiagonal, 1, -1)
        else
            b = op.space - nsites
            ss1,ss2 = op.subspace
            s1 = source(model, b)
            s2 = target(model, b)
            subspin1 = site2subspin(s1,ss1,S2)
            subspin2 = site2subspin(s2,ss2,S2)
            if op.let_type == LET_FMLink || op.let_type == LET_AFLink
                c = unify!(uf, currents[subspin1], currents[subspin2])
                currents[subspin1] = currents[subspin2] = c
                op.bottom_id = op.top_id = c
            elseif op.let_type == LET_Cross
                op.bottom_id = currents[subspin1]
                op.top_id = currents[subspin2]
                spins[subspin1], spins[subspin2] = spins[subspin2], spins[subspin1]
                currents[subspin1], currents[subspin2] = currents[subspin2], currents[subspin1]
            else # if op.let_type == LET_Vertex
                unify!(uf, currents[subspin1], currents[subspin2])
                op.bottom_id = currents[subspin1]
                c = addnode!(uf)
                op.top_id = c
                currents[subspin1] = currents[subspin2] = c
                spins[subspin1] *= ifelse(op.isdiagonal, 1, -1)
                spins[subspin2] *= ifelse(op.isdiagonal, 1, -1)
            end
        end
    end # of while loop

    ## PBC for imaginary time axis
    subspin = 0
    for s in 1:nsites
        ups0 = zeros(Int,0)
        ups1 = zeros(Int,0)
        downs0 = zeros(Int,0)
        downs1 = zeros(Int,0)
        for ss in 1:S2
            subspin += 1
            push!(ifelse(spins[subspin]==1, ups0, downs0), subspin)
            push!(ifelse(model.spins[subspin]==1, ups1, downs1), subspin)
        end
        @assert length(ups0) == length(ups1)
        @assert length(downs0) == length(downs1)
        for (u, u2) in zip(ups0, myshuffle(rng,ups1))
            unify!(uf, u, currents[u2])
        end
        for (d, d2) in zip(downs0, myshuffle(rng,downs1))
            unify!(uf, d, currents[d2])
        end
    end

    model.ops = ops

    nc = clusterize!(uf)
    flips = rand(rng, [1,-1], nc)
    for s in 1:nspins
        model.spins[s] *= flips[clusterid(uf, s)]
    end
    for op in model.ops
        if op.let_type == LET_Cross || op.let_type == LET_Vertex || op.let_type == LET_Cut
            bid = clusterid(uf, op.bottom_id)
            tid = clusterid(uf, op.top_id)
            op.isdiagonal ⊻= (flips[bid] != flips[tid])
        end
    end

    return uf
end

@gen_convert_parameter(QuantumXXZ, ("T", 1, 1.0),
                                   ("Jz", numbondtypes, 1.0),
                                   ("Jxy", numbondtypes, 1.0),
                                   ("Gamma", numsitetypes, 0.0),
                                  )


text:
polynomialize_vars(x, u...) = Dict(v.val => PolyVar{true}(string(v.val.name)) for v in vcat(x,u...))
polynomialize_vars(x) = Dict(v.val => PolyVar{true}(string(v.val.name)) for v in x)

# need to find a better check if expression has been expanded. Here I assume that every expression will be supplied as Num.
# If expression is supplied as NON-EXPANDED SymbolicUtils object, this whole transformation may fail
# For now good enough.
polynomialize_expr(ex::Num, vars::Dict) = polynomialize_expr(expand(ex).val, vars)
polynomialize_expr(ex::Pow, vars::Dict) = vars[ex.base]^ex.exp
polynomialize_expr(ex::Mul, vars::Dict) = prod(polynomialize_expr(sub_ex, vars) for sub_ex in arguments(ex))
polynomialize_expr(ex::Add, vars::Dict) = sum(polynomialize_expr(sub_ex, vars) for sub_ex in arguments(ex))
polynomialize_expr(ex::Sym, vars::Dict) = vars[ex]
polynomialize_expr(ex::Number, vars::Dict) = ex

function polynomialize_expr(exs::Array, vars::Dict)
    poly_exs = Array{Polynomial{true, Float64}}(undef, size(exs))
    for (i, ex) in enumerate(exs)
        poly_exs[i] = polynomialize_expr(ex, vars)
    end
    return poly_exs
end

function polynomialize_set(S::Vector, vars::Dict)
    ineq = Polynomial{true,Float64}[]
    eq = Polynomial{true,Float64}[]
    for entry in S
        polys = polynomialize_expr(entry.val.arguments, vars)
        if entry.val.f == >=
            push!(ineq, polys[1] - polys[2])
        elseif entry.val.f == <=
            push!(ineq, polys[2] - polys[1])
        else
            push!(eq, polys[1] - polys[2])
        end
    end
    return BasicSemialgebraicSet(algebraicset(eq), ineq)
end

polynomialize_set(S::Num, vars::Dict) = polynomialize_set([S], vars)

## Process Constructors
# Jump Processes
function JumpProcess(x::Vector{Num}, a::Vector{Num}, h::Vector{Vector{Num}}, X = []; iv = [], controls = [])
    poly_vars = polynomialize_vars(x, controls, iv)
    a = polynomialize_expr(a, poly_vars)
    h = [polynomialize_expr(hi, poly_vars) for hi in h]
    X_poly = isempty(X) ? FullSpace() : polynomialize_set(X, poly_vars)
    poly_x = [poly_vars[state] for state in x]
    poly_iv = isempty(iv) ? PV{true}("t") : poly_vars[iv]
    poly_controls = isempty(controls) ? PV{true}[] : [poly_vars[u] for u in controls]
    return JumpProcess(poly_x, a, h, X_poly; poly_vars = poly_vars, iv = poly_iv, controls = poly_controls)
end

JumpProcess(x::Num, a::Num, h::Num, X = []; iv = [], controls = []) =
            JumpProcess([x],[a],[[h]],X; iv = iv, controls = controls)
JumpProcess(x::Num, a::Vector{Num}, h::Vector{Num}, X = []; iv = [], controls = []) =
            JumpProcess([x],a,[[hi] for hi in h],X; iv = iv, controls = controls)
JumpProcess(x::Num, a::Vector{Num}, h::Vector{Vector{Num}}, X = []; iv = [], controls = []) =
            JumpProcess([x],a,h,X; iv = iv, controls = controls)

# Diffusion Processes
function DiffusionProcess(x::Vector{Num}, f::Vector{Num}, σ::Matrix{Num}, X = []; iv = [], controls = [])
    poly_vars = polynomialize_vars(x, controls, iv)
    f = polynomialize_expr(f, poly_vars)
    σ = polynomialize_expr(σ, poly_vars)
    X_poly = isempty(X) ? FullSpace() : polynomialize_set(X, poly_vars)
    poly_x = [poly_vars[state] for state in x]
    poly_iv = isempty(iv) ? PV{true}("t") : poly_vars[iv]
    poly_controls = isempty(controls) ? PV{true}[] : [poly_vars[u] for u in controls]
    return DiffusionProcess(poly_x, f, σ, X_poly, poly_vars = poly_vars, iv = poly_iv, controls = poly_controls)
end

DiffusionProcess(x::Num, f::Num, σ::Num, X = []; iv = [], controls = []) =
                 DiffusionProcess([x], [f], reshape([σ],1,1), X; iv = iv, controls = controls)

# Jump-Diffusion process
function JumpDiffusionProcess(x::Vector{Num}, a::Vector{Num}, h::Vector{Vector{Num}}, f::Vector{Num}, σ::Matrix{Num}, X = []; iv = [], controls = [])
    poly_vars = polynomialize_vars(x, iv, controls)
    a = polynomialize_expr(a, poly_vars)
    h = [polynomialize_expr(hi, poly_vars) for hi in h]
    f = polynomialize_expr(f, poly_vars)
    σ = polynomialize_expr(σ, poly_vars)
    X_poly = isempty(X) ? FullSpace() : polynomialize_set(X, poly_vars)
    poly_x = [poly_vars[state] for state in x]
    poly_iv = isempty(iv) ? PV{true}("t") : poly_vars[iv]
    poly_controls = isempty(controls) ? PV{true}[] : [poly_vars[u] for u in controls]
    return JumpDiffusionProcess(poly_x, a, h, f, σ, X_poly, poly_vars = poly_vars, iv = poly_iv, controls = poly_controls)
end

JumpDiffusionProcess(x::Num, a::Num, h::Num, f::Num, σ::Num, X = []; iv = [], controls = []) =
                     JumpDiffusionProcess([x], [a], [[h]], [f], reshape([σ],1,1), X; iv = iv, controls = controls)
JumpDiffusionProcess(x::Num, a::Vector{Num}, h::Vector{Num}, f::Num, σ::Num, X = []; iv = [], controls = []) =
                     JumpDiffusionProcess([x], a, [[hi] for hi in h], [f], reshape([σ],1,1), X; iv = iv, controls = controls)


text:
#############################################################################
# prox_convex.jl
# compute the prox of convex functions
#############################################################################

using Mosek

export prox

function prox(f::AbstractExpr, w::AbstractArray, abs_tol::Float64 = 1e-3)
  x = get_cv(f)[1]
  problem = Problem(:minimize, f + 0.5 * norm(x - w)^2)
  solve!(problem, MosekSolver(MSK_IPAR_LOG = 0, MSK_DPAR_INTPNT_CO_TOL_REL_GAP = abs_tol))
  return(x.value)
end


text:
module Tutorial

# Import the Mapper package.
using Mapper2

# Much of the exact behavior of the Mapper is controlled by parameterizing the
# TopLevel and Map types with a subtype of Architecture.
#
# This declarations allows us to extend various methods of placement and routing,
# should we want.
#
# If we don't extend those, default behavior will be used.
struct TutorialRuleset <: RuleSet end

# A tutorial on how to model a simple architecture.
include("Architecture.jl")

# Building a taskgraph for this simple architecture.
include("Taskgraph.jl")

# With the Architecture and Taskgraph defined, creating a Map is easy!.
function make_map(; width = 4, height = 4, ntasks = 10, nedges = 20)
    arch = build_arch(width, height)
    taskgraph = make_taskgraph(ntasks, nedges)

    return Map(TutorialRuleset(), arch, taskgraph)
end

end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("car", "SLID"), x="Wages", color="Language", Geom.histogram)


text:
## First Brillouin Zone

#function fbz_filter(n)
#  ks_set = Set()
#  for g=1:length(g_neighbors)
#    if dot(n,collect(g_neighbors[g]))/2 < norm(g_neighbors[g])/4
#    then
#     push!(ks_set,n)
#    end
#  end
#  return collect(ks_set)
#end


text:
function update_nu!(state::State, data::Data, prior::Prior, tuners::Tuners)
  # NOTE: This could be optimized.
  for k in eachindex(state.nu)
    update_nu!(k, state, data, prior, tuners)
  end
end

function update_nu!(k::Int, state::State, data::Data, prior::Prior,
                    tuners::Tuners)
  log_prob = function(nu_k)
    logprior = logpdf(prior.nu, nu_k)

    loglike = zero(state.p)
    for i in samplenames
      lami = ref_lambda(state, i)
      vi = ref_v(state, i)
      for n in eachindex(vi)
        (lami[n] == k) && (loglike += gammalogpdf(nu_k/2, 2/nu_k, vi[n]))
      end
    end

    return logprior + loglike
  end

  state.nu[k] = MCMC.metLogAdaptive(state.nu[k], log_prob, tuners.nu[k])
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
JSON.json(jws::JSONWorksheet) = JSON.json(data(jws))
JSON.json(jws::JSONWorksheet, indent) = JSON.json(data(jws), indent)

function dropnull(s)
    replace(s, r"(\"[\w]*\":null,)|(,?\"[\w]*\":null)" => "")
end

function write(file::Union{String, IO}, jws::JSONWorksheet; indent = 2, drop_null = false)
    open(file, "w") do io
        data = JSON.json(jws, indent)
        if drop_null
            data = dropnull(data)
        end

        Base.write(io, data)
    end
end

function write(path::String, jwb::JSONWorkbook; kwargs...)
    f = splitext(basename(xlsxpath(jwb)))[1]
    for s in sheetnames(jwb)
        write(joinpath(path, "$(f)_$(s).json"), jwb[s]; kwargs...)
    end
end


text:
module CommunalHelperCoreModeMusicController

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/CoreModeMusicController" Controller(
    x::Int,
    y::Int,
    params::String="",
    hot::Number=1.0,
    cold::Number=0.0,
    none::Number=0.0,
    disable::Bool=false
)

const placements = Ahorn.PlacementDict(
    "Core Mode Music Controller (Communal Helper)" => Ahorn.EntityPlacement(
        Controller,
    ),
    "Core Mode Music Controller (Communal Helper) (Disable)" => Ahorn.EntityPlacement(
        Controller,
        "point",
        Dict{String, Any}(
            "disable" => true,
        ),
    ),
)

const spriteEnable = "objects/CommunalHelper/coreModeMusicController/iconEnable"
const spriteDisable = "objects/CommunalHelper/coreModeMusicController/iconDisable"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Controller) = Ahorn.drawImage(ctx, Bool(get(entity.data, "disable", false)) ? spriteDisable : spriteEnable, -8, -11)

function Ahorn.selection(entity::Controller)
    x, y = Ahorn.position(entity)
    return Ahorn.Rectangle(x - 8, y - 11, 22, 22)
end

end


text:
get_parameter(sdf::SDFFile, p::Symbol) = getindex(sdf.param[], p)
get_parameter(sdf::SDFFile, p::Symbol, c::Symbol) = getindex(get_parameter(sdf, p), c)

get_time(sdf::SDFFile) = sdf.header.time * u"s"
get_npart(sdf::SDFFile, species) = sdf.blocks[Symbol("grid/"*species)].np

function Base.ndims(sdf::SDFFile)
    ks = string.(keys(sdf))
    i = findfirst(k->startswith(k, "grid"), ks)
    ndims(sdf.blocks[Symbol(ks[i])])
end

function get_data_description(file)
    h, blocks = open(file_summary, file)
    names = string.(keys(blocks))
    idx = findfirst(b->occursin("grid/", b), names)
    if isnothing(idx)
        @debug "No particle data found"
        N = 1
        T = Float64
    else
        N = ndims(blocks[Symbol(names[idx])])
        T = eltype(blocks[Symbol(names[idx])])
    end

    return N, T
end


text:
struct Reverse{D} end

"""
    Reverse{d}(A)
Lazy version of `reverse(A, dims=d)`.
"""
function Reverse{D}(A::AbstractArray{T,N}) where {D,T,N}
    1 ≤ D ≤ N || throw(ArgumentError("dimension $D is not 1 ≤ $D ≤ $N"))
    tup = ntuple(i -> i==D ? reverse(axes(A,D)) : Colon(), N)
    view(A, tup...)
end

Reverse(A::AbstractArray; dims=ndims(A)) = Reverse{dims}(A)

# The reverse() function in Base is twice as quick as copying this view,
# somewhere I had a real type in order to exploit this?

using Random

struct Shuffle{D} end

"""
    Shuffle{d}(A)
For a vector this is a lazy version of `shuffle(A)`,
for a tensor this shuffles along one axis only.
"""
function Shuffle{D}(A::AbstractArray{T,N}) where {D,T,N}
    1 ≤ D ≤ N || throw(ArgumentError("dimension $D is not 1 ≤ $D ≤ $N"))
    tup = ntuple(i -> i==D ? shuffle(axes(A,D)) : Colon(), N)
    view(A, tup...)
end

Shuffle(A::AbstractArray; dims=ndims(A)) = Shuffle{dims}(A)


text:
# __precompile__()

module OnlineAI

using Reexport
using Distributions
@reexport using QuickStructs
@reexport using OnlineStats
#@reexport using Qwt
#@reexport using CTechCommon
# using StatsBase
@reexport using StatsBase
# using Plots
using Requires
using ArrayViews
using Plots
@reexport using LearnBase
import LearnBase: value

import OnlineStats: nrows, ncols,
                    VecF, MatF, AVecF, AMatF,
                    standardize

export row, col, row!, col!, rows, cols, nrows, ncols,
       VecF, MatF, AVecF, AMatF,
       standardize

# represents a node in an arbitrary graph... typically representing a neuron within a neural net
abstract NetStat <: OnlineStat
nobs(o::NetStat) = 0
abstract AbstractNeuron

abstract NeuralNetLayer

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

include("utils.jl")

# export
#   Mapping,
#   IdentityMapping,
#   SigmoidMapping,
#   TanhMapping,
#   SoftsignMapping,
#   ReLUMapping,
#   LReLUMapping,
#   SoftmaxMapping
# include("nnet/activations.jl")

export
  DataPoint,
  Transformation,
  IdentityTransform,
  AbsTransform,
  LogPlus1Transform,
  SquareTransform,
  CubeTransform,
  SignSquareTransform,
  Transformer,
  IdentityTransformer,
  VectorTransformer,
  transform,
  transform!,
  DataPoints,
  splitDataPoints,
  DataSampler,
  SimpleSampler,
  SubsetSampler,
  splitDataSamplers,
  StratifiedSampler,
  crossValidationSets
include("nnet/data.jl")

export
  cost,
  totalCost
include("nnet/costs.jl")

export
  current_updater,
  current_mloss,
  current_ploss,
  current_updater!,
  current_mloss!,
  current_ploss!,
  LearningRateModel,
  FixedLearningRate,
  AdaptiveLearningRate
include("nnet/gradient.jl")

export
  DropoutStrategy,
  Dropout,
  NoDropout,
  NetParams
include("nnet/params.jl")

export
  SolverParams,
  SolverStats,
  solve!
include("nnet/solver.jl")

export
  NeuralNetLayer,
  Layer,
  # forward,
  # backward,
  forward!,
  backward!
include("nnet/layer.jl")

export
  NormalizedLayer
include("nnet/normlayer.jl")

export
  NeuralNet
include("nnet/net.jl")

export
  pretrain
include("nnet/pretrain.jl")

# include("nnet/lstm.jl")

export
  buildNet,
  buildClassificationNet,
  buildTanhClassificationNet,
  buildRegressionNet
include("nnet/build.jl")

export
  Constant,
  HiddenLayerSampler,
  VectorSampler,
  ParameterSampler,
  generateTransformer,
  generateModels,
  Ensemble
include("nnet/ensembles.jl")

# export
#   visualize,
#   track_progress

# function __init__()

# include("nnet/visualize.jl")

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

# NOTE: EXPERIMENTAL

export
  GaussianReceptiveField,
  # value,
  Synapse,
  DelaySynapse,
  fire!,    # checks for threshold crossing, then fires
  SpikingNeuron,
  DiscreteLeakyIntegrateAndFireNeuron,
  LiquidParams,
  Liquid,

  ImmediateSynapse,
  GRFNeuron,
  GRFInput,
  LiquidInput,
  LiquidInputs,

  Readout,
  FireReadout,
  StateReadout,
  FireWindowReadout,

  LiquidStateMachine,
  liquidState,

  LiquidVisualization,
  LiquidVisualizationNode

abstract Synapse
abstract SpikingNeuron <: AbstractNeuron
abstract LiquidInput

include("spiking/readout.jl")
include("spiking/liquid.jl")
include("spiking/srm.jl")
include("spiking/input.jl")
# include("spiking/visualize.jl")

include("spiking/skan.jl")

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

end


text:
push!(LOAD_PATH,"../src/")

using Documenter, CellMLToolkit

makedocs(
    sitename= "CellMLToolkit",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md"
    ]
)


text:
# Autogenerated wrapper script for LHAPDF_jll for i686-w64-mingw32-cxx11
export libLHAPDF

JLLWrappers.@generate_wrapper_header("LHAPDF")
JLLWrappers.@declare_library_product(libLHAPDF, "libLHAPDF-0.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libLHAPDF,
        "bin\\libLHAPDF-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Adapted from:
# SOSDEMO5 --- Upper bound for the structured singular value mu
# Section 3.5 of SOSTOOLS User's Manual

@testset "SOSDEMO5 with $(factory.constructor)" for factory in sdp_factories
    isscs(factory) && continue
    iscsdp(factory) && continue # See https://github.com/JuliaOpt/SumOfSquares.jl/issues/52
    @polyvar x[1:8]

    # The matrix under consideration
    alpha = 3 + sqrt(3);
    beta = sqrt(3) - 1;
    a = sqrt(2/alpha);
    b = 1/sqrt(alpha);
    c = b;
    d = -sqrt(beta/alpha);
    f = (1 + im)*sqrt(1/(alpha*beta));
    U = [a 0; b b; c im*c; d f];
    V = [0 a; b -b; c -im*c; -im*f -d];
    M = U*V';

    @testset "with γ=$γ it should be $(feasible ? "feasible" : "infeasible")" for (γ, feasible) in ((0.8722, false), (0.8724, true))
        Z = monomials(x, 1)

        function build_A(i)
            H = M[i,:]*M[i,:]' - (γ^2)*sparse([i],[i],[1],4,4)
            H = [real(H) -imag(H); imag(H) real(H)]
            dot(Z, H*Z)
        end
        A = build_A.(1:4)

        m = SOSModel(factory)

        # -- Q(x)'s -- : sums of squares
        # Monomial vector: [x1; ... x8]
        Q = Vector{GramMatrix{JuMP.VariableRef, monomialtype(x[1]),
                                 monovectype(x[1])}}(undef, 4)
        @variable m Q[1:4] SOSPoly(Z)

        # -- r's -- : constant sum of squares
        Z = monomials(x, 0)
        #r = Matrix{GramMatrix{JuMP.VariableRef}}(4,4) # FIXME doesn't work with 1x1 SDP matrix :(
        @variable m r[i=1:4,j=(i+1):4] >= 0

        # Constraint : -sum(Qi(x)*Ai(x)) - sum(rij*Ai(x)*Aj(x)) + I(x) >= 0
        expr = 0
        # Adding term
        for i = 1:4
            expr -= A[i]*Q[i]
        end
        for i = 1:4
            for j = (i+1):4
                expr -= A[i]*A[j]*r[i,j]
            end
        end
        # Constant term: I(x) = -(x1^4 + ... + x8^4)
        I = -sum(x.^4)
        expr = expr + I

        @constraint m expr >= 0

        JuMP.optimize!(m)

        # Program is feasible, thus 0.8724 is an upper bound for mu.
        if feasible
            @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        else
            @test JuMP.dual_status(m) == MOI.INFEASIBILITY_CERTIFICATE
        end
    end
end


text:
import Criterion
using FactCheck

include("stats.jl")
include("analysis.jl")


text:
module LigParGen2LAMMPS

include("read_data.jl")
include("print_mol.jl")

end # module


text:
using AugmentedGaussianProcesses;
using Distributions, LinearAlgebra
using Random: seed!;
seed!(123);
using Plots;
gr();

N = 100;
noise = 1e-1;
N_test = 1000
X = sort(rand(N))
X_test = collect(range(0; length=N_test, stop=1))
k = SqExponentialKernel() ∘ ScaleTransform(10.0)
K = kernelmatrix(k, X) + noise * I
f = rand(MvNormal(K))
λ = 10
y = rand.(Poisson.(λ * AGP.logistic.(f)))

poissonmodel = VGP(X, y, k, PoissonLikelihood(Float64(λ)), AnalyticVI(); optimiser=false);
train!(poissonmodel, 300)
poissonpred = predict_y(poissonmodel, X_test)

negbinmodel = VGP(X, y, k, NegBinomialLikelihood(λ), AnalyticVI(); optimiser=false);
train!(negbinmodel, 50)
negbinpred = predict_y(negbinmodel, X_test)
##
default(; legendfontsize=14.0, xtickfontsize=10.0, ytickfontsize=10.0)
p = plot(X, y; t=:scatter, lab="Training Points")
plot!(X_test, poissonpred; lab="Poisson Expectation", lw=7.0)
plot!(X_test, negbinpred; lab="Negative Binomial Expectation", lw=7.0, legend=:topright)

display(p)
savefig(p, joinpath(@__DIR__, "Events.png"))
# a,b = logitmodel.fstar(X_test)
# using Expectations, QuadGK, BenchmarkTools

# @btime begin; for i in 1:N_test; d = Normal(a[i],sqrt(max(b[i],1e-8))); expectation(logit,d); end; end;
# @btime begin; for i in 1:N_test; d = Normal(a[i],sqrt(max(b[i],1e-8))); quadgk(x->pdf(d,x)*logit(x),-Inf,Inf); end; end;


text:
using PolyChaos, Test, LinearAlgebra
import QuadGK: quadgk

function integration_test(op::AbstractOrthoPoly, dim::Int,name::String;tol::Float64=1e-6,
                          dom::Tuple{<:Real,<:Real}=op.measure.dom)
    N, w, ind = op.deg, op.measure.w, zeros(Int64,dim)
    @testset "$name" begin
        for ind_ in Iterators.product([collect(1:N) for i in 1:dim]...)
            ind[:] .= ind_[:]
            s1 = computeSP(ind, op)
            f(t) = prod(evaluate(ind[i], t, op) for i in 1:dim)*w(t)
            s2 = quadgk(f, dom[1], dom[2])[1]
            if abs(s1)<=1e-3
                @test isapprox(s1, s2; atol=tol)
            else
                @test isapprox(abs(s1/s2), 1; atol=tol)
            end
        end
    end
end

N, Nrec, dim, tol = 3, 1000, 3, 1e-4

op = GaussOrthoPoly(N;Nrec=Nrec)
set_gaussian = integration_test(op,dim,"gaussian";tol=tol,dom=(-10.,10.))

α, β = 1.32, 3.23
op = Beta01OrthoPoly(N, α, β; Nrec=Nrec)
set_beta = integration_test(op, dim, "beta01"; tol=tol)

op = GammaOrthoPoly(N, 1, 1; Nrec=Nrec)
set_gamma = integration_test(op,dim,"gamma";tol=tol)


text:
# Autogenerated wrapper script for MKL_jll for i686-w64-mingw32
export libmkl_core, libmkl_rt

using IntelOpenMP_jll
JLLWrappers.@generate_wrapper_header("MKL")
JLLWrappers.@declare_library_product(libmkl_core, "mkl_core.1.dll")
JLLWrappers.@declare_library_product(libmkl_rt, "mkl_rt.1.dll")
function __init__()
    JLLWrappers.@generate_init_header(IntelOpenMP_jll)
    JLLWrappers.@init_library_product(
        libmkl_core,
        "bin\\mkl_core.1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmkl_rt,
        "bin\\mkl_rt.1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()