text:
using BITEModel
using Base.Test
const BM = BITEModel

## Forward model
################

plotyes = false
verbose = true

global gid, hs1d, ivs1d, hs2d, ivs2d, glaciertype, # this is a bit ugly but
       gl, gb, pp, pm, pn                          # needed because of
                                                   # include("../scripts/forward.jl")

@testset "Forward model with update_cache=$update_cache" for update_cache in [false,true]
    # Synthetic (run by including scripts/forward.jl)
    ###########
    glaciertype = :synth
    plotyes = false
    include("../scripts/forward.jl")
    Base.Test.@inferred BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)
    Base.Test.@inferred BM.make_bands(gl, pp, pm, pn)
    fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

    hs1d_ = [0.0, 167.025, 182.466, 189.223, 188.814, 183.506, 175.65,
             166.279, 151.378, 141.485, 131.018, 119.681, 109.229, 98.3605,
             86.7481, 76.5545, 67.8539, 63.5029, 57.3474, 50.484, 1.0]
    ivs1d_ = [0.0, 3.3608, 5.31713, 6.91192, 8.26612, 9.03612, 9.67009,
              10.0828, 10.5786, 11.0246, 10.9542, 10.8031, 10.8764, 10.6826,
              9.68466, 8.8367, 7.85266, 6.36658, 4.69378, 2.97981, 0.0]
    @test norm(hs1d_-fwdsol.hs1d)<2e-3
    @test norm(ivs1d_-fwdsol.ivs1d)<5e-4

    # test 2D at some random indices
    @test length(fwdsol.hs2d)==15921
    inds = [8179, 15516, 3523, 5846, 3236, 13314, 4959, 11, 12505, 6528,
            11823, 13553, 1043, 1420, 4210, 1991, 34, 4056, 5848, 15447,
            13452, 390, 8071, 9023, 12655, 3717, 7865, 8238, 14551, 4782]
    hs2d_ = [176.314, 0.0, 129.79, 172.317, 108.296, 0.0, 0.0, 0.0,
             57.2049, 0.0, 104.003, 0.0, 0.0, 0.0, 55.4563, 0.0, 0.0, 161.232,
             174.862, 0.0, 109.54, 0.0, 95.0336, 222.902, 119.131, 82.9983,
             98.8966, 230.184, 0.0, 132.357]
    ivs2d_ = [8.22742, 0.0, 6.17785, 7.93149, 6.84787, 0.0, 0.0, 0.0, 0.0,
              0.0, 7.88523, 0.0, 0.0, 0.0, 7.28846, 0.0, 0.0, 5.7638,
              7.93149, 0.0, 4.89699, 0.0, 0.0, 6.58108, 6.01924,
              7.96476, 6.96876, 6.70248, 0.0, 8.22056]
    @test norm(hs2d_-fwdsol.hs2d[inds])<1e-3
    @test norm(ivs2d_-fwdsol.ivs2d[inds])<5e-4

    # # ITMIX (run by hand)
    # #####################
    # glacier=:Urumqi # pretty fast this one
    # pl_kws = Dict{Symbol,Any}(:use_glogem => false)
    # gid = ITMIXGlacier(glacier)
    # gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose;
    #                                     update_cache=update_cache, pl_kws...)
    # hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb =
    #     BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)
    # Base.Test.@inferred BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)
    # Base.Test.@inferred BM.make_bands(gl, pp, pm, pn)

    # hs1d_ = [0.0, 25.4951, 29.2112, 32.2141, 32.3665, 30.4775, 29.6186,
    #          29.2802, 27.4351, 28.3908, 29.75, 37.2109, 39.2642,
    #          40.3066, 42.8661, 44.964 , 48.2977, 53.8874, 56.7994,
    #          57.8052, 51.7617, 40.3284, 31.8477, 5.63301]
    # ivs1d_ = [0.0, 1.47265, 2.2546, 3.01979, 3.92654, 4.90726,
    #           4.87954, 4.17233, 3.82581, 3.51174, 3.2171, 2.61505, 2.57666,
    #           2.68062, 2.84157, 3.30573, 3.93735, 4.08862, 4.08278, 4.08057,
    #           4.2603, 4.31557, 2.45511, 0.00525341]
    # @test norm(hs1d_-hs1d)<5e-4
    # @test norm(ivs1d-ivs1d_)<5e-4

    # # test 2D at some random indices
    # @test length(hs2d)==133331
    # inds = [85082, 118470, 7808, 31087, 78719, 96341, 104214, 9715,
    #         91721, 131435, 49463, 127799, 36643, 131306, 71573, 88897,
    #         110465, 58186, 93094, 29907, 105007, 65361, 118730, 67826,
    #         118921, 66973, 114737, 84487, 19245, 126314]
    # hs2d_ = [0.0, 54.1206, 0.0, 0.0, 39.7876, 0.0, 0.0, 39.0787,
    #          23.5233, 0.0, 0.0, 19.2427, 0.0, 0.0, 22.237, 50.1744,
    #          0.0, 71.186, 0.0, 35.3194, 27.0517, 77.9176, 0.0, 21.425,
    #          0.0, 0.0, 62.8376, 82.3543, 43.1886, 28.9674]
    # ivs2d_ = [0.0, 1.73392, 0.0, 0.0, 3.50925, 0.0, 0.0, 1.3902,
    #           5.06693, 0.0, 0.0, 3.04547, 0.0, 0.0, 2.57997, 3.94191,
    #           0.0, 2.76635, 0.0, 1.6151, 6.13513, 3.14324, 0.0,
    #           2.31482, 0.0, 0.0, 1.87241, 4.19182, 1.40529, 1.68439]
    # @test norm(hs2d_-hs2d[inds])<5e-4
    # @test norm(ivs2d_-ivs2d[inds])<5e-4
end

# make sure this runs
include("../scripts/forward-advanced.jl")

## Inverse model
################
@testset "Inverse model" begin
    gid = SyntheticGlacier(:bench)
    gl,gb,pp,pm,pn,pl = BM.init_forward(gid)
    (theta0, logposterior, logposterior1d, logprior, logprior_debug,
     loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
     pmcmc_defaults, fit_target) =
         BM.init_inverse(gb, pp, pm, pn)
    Base.Test.@inferred fwdm_fn(theta0.th0)
    Base.Test.@inferred logprior(theta0.th0)
    Base.Test.@inferred loglikelihood(theta0.th0)
    Base.Test.@inferred logposterior(theta0.th0)

    pmcmc = BM.MCMCNum(;pmcmc_defaults...,
                       niter=5*10^4,
                       nthin = 50)

    # sample posterior
    srand(0)
    thetas, theta_expect, theta_mode, accept_ratio, blobs, lpost, mc_hs1d, mc_ivs1d,
      mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d = mcmc(logposterior, theta0, pmcmc; verbose=true)
    @test round(accept_ratio,5) == 0.11801

    @test norm(theta_expect- [1.81352, 1.40765, -0.29361, -0.170199,
                              -0.0292185, -0.15876, -0.3106, -0.165316, -0.0345219, 0.172756,
                              0.406927, 0.162265, -0.144996, 0.315479, 0.16641, -0.00986009,
                              0.0872772, -0.109316, 0.711055, -0.239801, -0.121702, -0.0558777,
                              -0.219639]) < 1e-5
    # test 2D at some random indices
    inds = [8179, 15516, 3523, 5846, 3236, 13314, 4959, 11, 12505, 6528,
            11823, 13553, 1043, 1420, 4210, 1991, 34, 4056, 5848, 15447,
            13452, 390, 8071, 9023, 12655, 3717, 7865, 8238, 14551, 4782]

    @test norm(mh2d[inds]-[195.185, 0.0, 112.593, 175.143, 110.134, 0.0, 0.0, 0.0, 48.5617,
                           0.0, 115.334, 0.0, 0.0, 0.0, 48.79, 0.0, 0.0, 127.909, 173.234,
                           0.0, 86.9006, 0.0, 80.6748, 169.704, 104.964, 85.447, 86.3577,
                           180.154, 0.0, 147.823])<2e-3

    @test norm(sh2d[inds]-[14.5854, 0.0, 21.0239, 18.2294, 11.4231, 0.0, 0.0, 0.0, 9.05472,
                           0.0, 10.1193, 0.0, 0.0, 0.0, 19.4267, 0.0, 0.0, 24.128, 21.1017,
                           0.0, 16.3925, 0.0, 15.0425, 27.2756, 19.6065, 9.78661, 33.134,
                           32.304, 0.0, 11.914])<2e-4

    @test norm(miv2d[inds]-[11.8817, 0.0, 8.42797, 10.3403, 10.106,
                            0.0, 0.0, 0.0, 0.0, 0.0, 12.5615, 0.0, 0.0, 0.0, 16.4126, 0.0,
                            0.0, 7.73834, 10.3403, 0.0, 7.61675, 0.0, 0.0, 7.85299, 8.40403,
                            14.5605, 16.3398, 8.09658, 0.0, 12.651])<1e-4
end


text:
# This file is a part of H3DDetectorSystems.jl, licensed under the MIT License (MIT).

struct PolarisEventHeader
    timestamp::UInt64     # Event timestamp
end

function Base.read(src::IO, ::Type{PolarisEventHeader})
    timestamp = ntoh(read(src, UInt64))

    PolarisEventHeader(timestamp)
end

struct PolarisHit
    detno::Int32 # Detector number
    x::Int32  # Position in µm,
    y::Int32  # Position in µm,
    z::Int32  # Position in µm,
    edep::Int32 # Energy deposition in eV
end

function Base.read(src::IO, ::Type{PolarisHit})
    detno = ntoh(read(src, UInt8))
    x = ntoh(read(src, Int32)) # µm,
    y = ntoh(read(src, Int32)) # µm,
    z = ntoh(read(src, Int32)) # nm, zero at anode
    edep = ntoh(read(src, Int32)) # eV

    PolarisHit(detno, x, y, z, edep)
end

struct PolarisSyncHeader
    n::Int64
end

function Base.read(src::IO, ::Type{PolarisSyncHeader})
    n = ntoh(read(src, UInt8))  # number of PolarisSyncValue entries to follow
    PolarisSyncHeader(n)
end

struct PolarisSyncValue
    x::UInt64
end

function Base.read(src::IO, ::Type{PolarisSyncValue})
    x = ntoh(read(src, Int64))
    PolarisSyncValue(x)
end

ext: jl


text:
"calculates the entropy of a given state"
entropy(state::InfiniteMPS) = [sum([-j^2*2*log(j) for j in entanglement_spectrum(state,i)]) for i in 1:length(state)]

"
given a thermal state, you can map it to an mps by fusing the physical legs together
to prepare a gibbs ensemble, you need to evolve this state with H working on both legs
here we return the 'superhamiltonian' (H*id,id*H)
"
function splitham(ham::MPOHamiltonian)
    fusers = [isomorphism(fuse(p*p'),p*p') for p in ham.pspaces]

    idham = Array{Union{Missing,typeof(ham[1,1,1])},3}(missing,ham.period,ham.odim,ham.odim)
    hamid = Array{Union{Missing,typeof(ham[1,1,1])},3}(missing,ham.period,ham.odim,ham.odim)

    for i in 1:ham.period
        for (k,l) in keys(ham,i)
            idt = isomorphism(ham.pspaces[i],ham.pspaces[i])

            @tensor hamid[i,k,l][-1 -2; -3 -4] :=   ham[i,k,l][-1,12,-3,15]*
                                                    idt[16,13]*
                                                    fusers[i][-2,12,13]*
                                                    conj(fusers[i][-4,15,16])

            @tensor idham[i,k,l][-1 -2; -3 -4] :=   ham[i,k,l][-1,13,-3,16]*
                                                    idt[12,15]*
                                                    fusers[i][-2,12,16]*
                                                    conj(fusers[i][-4,15,13])
        end
    end

    return MPOHamiltonian(hamid),MPOHamiltonian(idham)
end

infinite_temperature(ham::MPOHamiltonian) = [permute(isomorphism(Matrix{eltype(ham[1,1,1])},oneunit(sp)*sp,oneunit(sp)*sp),(1,2,4),(3,)) for sp in ham.pspaces]

"calculates the galerkin error"
calc_galerkin(state::Union{InfiniteMPS,FiniteMPS,MPSComoving},loc,envs) = norm(leftnull(state.AC[loc])'*ac_prime(state.AC[loc], loc,state,envs))
calc_galerkin(state::Union{InfiniteMPS,FiniteMPS,MPSComoving}, envs) = maximum([calc_galerkin(state,loc,envs) for loc in 1:length(state)])
calc_galerkin(state::MPSMultiline, envs::PerMPOInfEnv) = maximum([norm(leftnull(state.AC[row+1,col])'*ac_prime(state.AC[row,col], row,col,state,envs)) for (row,col) in Iterators.product(1:size(state,1),1:size(state,2))][:])
calc_galerkin(state::MPSMultiline, envs::MixPerMPOInfEnv) = maximum([norm(leftnull(state.AC[row+1,col])'*ac_prime(envs.above.AC[row,col], row,col,state,envs)) for (row,col) in Iterators.product(1:size(state,1),1:size(state,2))][:])

"
Calculates the (partial) transfer spectrum
"
function transfer_spectrum(above::InfiniteMPS;below=above,tol=Defaults.tol,num_vals = 20,sector=first(sectors(oneunit(virtualspace(above,1)))))
    init = TensorMap(rand, eltype(eltype(above)), virtualspace(below,0)*ℂ[sector => 1],virtualspace(above,0))

    num_vals = min(dim(virtualspace(above,0)*virtualspace(below,0)),num_vals); # we can ask at most this many values

    eigenvals, eigenvecs,convhist = eigsolve(x->transfer_left(x, above.AL, below.AL) , init, num_vals, :LM, tol=tol)
    convhist.converged < num_vals && @warn "correlation length failed to converge $(convhist.normres)"

    return eigenvals
end

"
Returns the (full) entanglement spectrum at site I
"
function entanglement_spectrum(st::Union{InfiniteMPS,FiniteMPS,MPSComoving},site::Int=0)
    @assert site<=length(st)

    (U,S,V) = tsvd(st.CR[site]);
    diag(convert(Array,S))
end

function correlation_length(above::InfiniteMPS; tol_angle=0.1,below=above,kwargs...)
    #get the transfer spectrum
    spectrum = transfer_spectrum(above;below=above,kwargs...);

    correlation_length(spectrum,above;tol_angle=tol_angle,below=below,kwargs...)
end

function correlation_length(spectrum,above::InfiniteMPS; tol_angle=0.1,below=above,kwargs...)
    #we also define a correlation length between different states
    (above === below) && (spectrum = spectrum[2:end])

    best_angle = mod1(angle(spectrum[1]), 2*pi)
    ind_at_angle = findall(x->x<tol_angle || abs(x-2*pi)<tol_angle, mod1.(angle.(spectrum).-best_angle, 2*pi))
    spectrum_at_angle = spectrum[ind_at_angle]

    lambdas = -log.(abs.(spectrum_at_angle));

    corlength = 1/first(lambdas);

    gap = Inf;
    if length(lambdas) > 2
        gap = lambdas[2]-lambdas[1]
    end

    return corlength, gap, best_angle
end

function variance(state::InfiniteMPS,ham::MPOHamiltonian,envs = environments(state,ham))
    rescaled_ham = ham-expectation_value(state,ham,envs);
    real(sum(expectation_value(state,rescaled_ham*rescaled_ham)))
end

function variance(state::FiniteMPS,ham::MPOHamiltonian,envs = environments(state,ham))
    ham2 = ham*ham;
    real(sum(expectation_value(state,ham2)) - sum(expectation_value(state,ham,envs))^2)
end

function variance(state::MPSComoving,ham::MPOHamiltonian,envs = environments(state,ham))
    #tricky to define
    (ham2,nenvs) = squaredenvs(state,ham,envs);
    real(expectation_value(state,ham2,nenvs)[2] - expectation_value(state,ham,envs)[2]^2)
end

variance(state::FiniteQP,ham::MPOHamiltonian,args...) = variance(convert(FiniteMPS,state),ham);

function variance(state::InfiniteQP,ham::MPOHamiltonian,envs=environments(state,ham))
    # I remember there being an issue here @gertian?
    state.trivial || throw(ArgumentError("variance of domain wall excitations is not implemented"));

    rescaled_ham = ham - expectation_value(state.left_gs,ham);

    #I don't remember where the formula came from
    E_ex = dot(state,effective_excitation_hamiltonian(ham,state,envs));
    E_f = expectation_value(state.left_gs,rescaled_ham,0);

    ham2 = rescaled_ham*rescaled_ham

    real(dot(state,effective_excitation_hamiltonian(ham2,state))-2*(E_f+E_ex)*E_ex+E_ex^2)
end

"""
You can impose periodic boundary conditions on an mpo-hamiltonian (for a given size)
That creates a new mpo-hamiltonian with larger bond dimension
The interaction never wraps around multiple times
"""
function periodic_boundary_conditions(ham::MPOHamiltonian{S,T,E},len = ham.period) where {S,T,E}
    sanitycheck(ham) || throw(ArgumentError("invalid ham"))
    mod(len,ham.period) == 0 || throw(ArgumentError("$(len) is not a multiple of unitcell"))

    fusers = PeriodicArray(map(1:len) do loc
        map(Iterators.product(ham.domspaces[len+1,:],ham.domspaces[loc,:],ham.domspaces[loc,:])) do (v1,v2,v3)
            isomorphism(fuse(v1'*v2*v3),v1'*v2*v3)
        end
    end)

    #a -> what virtual space did I "lend" in the beginning?
    #b -> what progress have I made in the lower layer?
    #c -> what progress have I made in the upper layer?
    χ = ham.odim;
    χ´ = Int((χ-1)*χ*(χ+1)/2+1);

    function indmap(a,b,c)
        Int((χ-a)*(χ-a+1)/2+(c-1)*χ*(χ+1)/2+(b-a)+1)
    end

    #do the bulk
    bulk = PeriodicArray(convert(Array{Union{T,E},3},fill(zero(E),ham.period,χ´,χ´)));

    for loc in 1:ham.period,
        (j,k) in keys(ham,loc)

        #apply (j,k) above
        l = ham.odim
        for i in 2:ham.odim

            k <= i && i<=l || continue

            f1 = fusers[loc][i,l,j]
            f2 = fusers[loc+1][i,l,k]
            @tensor bulk[loc,indmap(i,l,j),indmap(i,l,k)][-1 -2;-3 -4]:=ham[loc,j,k][1,-2,2,-4]*f1[-1,3,4,1]*conj(f2[-3,3,4,2])
        end

        #apply (j,k) below
        i = 1
        for l in 2:(ham.odim-1)

            l > 1 && l >= i && l<=j || continue

            f1 = fusers[loc][l,j,i];
            f2 = fusers[loc+1][l,k,i];

            @tensor bulk[loc,indmap(l,j,i),indmap(l,k,i)][-1 -2;-3 -4]:=ham[loc,j,k][1,-2,2,-4]*f1[-1,3,1,4]*conj(f2[-3,3,2,4])
        end
    end

    # make the starter
    starter = convert(Array{Union{T,E},2},fill(zero(E),χ´,χ´));
    for (j,k) in keys(ham,1)

        #apply (j,k) above
        if j == 1
            f1 = fusers[1][end,end,1];
            f2 = fusers[1][end,end,k];

            @tensor starter[1,indmap(ham.odim,ham.odim,k)][-1 -2;-3 -4]:=
            ham[1,j,k][1,-2,2,-4]*f1[-1,3,4,1]*conj(f2[-3,3,4,2])
        end

        #apply (j,k) below
        if j > 1 && j < ham.odim

            f1 = fusers[1][j,j,1];
            f2 = fusers[1][j,k,1];

            @tensor starter[1,indmap(j,k,1)][-1 -2;-3 -4]:=ham[1,j,k][1,-2,2,-4]*conj(f2[-3,1,2,-1])
        end
    end
    starter[1,1] = one(E);
    starter[end,end] = one(E);

    # make the ender
    ender = convert(Array{Union{T,E},2},fill(zero(E),χ´,χ´));
    for (j,k) in keys(ham,ham.period)

        if k > 1
            f1 = fusers[end][k,ham.odim,j]
            @tensor ender[indmap(k,ham.odim,j),end][-1 -2;-3 -4]:=ham[ham.period,j,k][1,-2,2,-4]*f1[-1,2,-3,1]
        end
    end
    ender[1,1] = one(E);
    ender[end,end] = one(E);

    # fill in the entire ham
    nos = convert(Array{Union{T,E},3},fill(zero(E),len,χ´,χ´))
    nos[1,:,:] = starter[:,:];
    nos[end,:,:] = ender[:,:];

    for i in 2:len-1
        nos[i,:,:] = bulk[i,:,:];
    end

    return MPOHamiltonian(nos)
end

#impose periodic boundary conditions on a normal mpo
function periodic_boundary_conditions(mpo::InfiniteMPO{O},len = length(mpo)) where O
    mod(len,length(mpo)) == 0 || throw(ArgumentError("len not a multiple of unitcell"))

    output = PeriodicArray{O,1}(undef,len);

    sp = space(mpo[1],1)';
    utleg = Tensor(ones,oneunit(sp));

    #do the bulk
    for j in 2:len-1
        f1 = isomorphism(fuse(sp*space(mpo[j],1)),sp*space(mpo[j],1))
        f2 = isomorphism(fuse(sp*space(mpo[j],3)'),sp*space(mpo[j],3)')

        @tensor output[j][-1 -2;-3 -4] := mpo[j][1,-2,2,-4]*f1[-1,3,1]*conj(f2[-3,3,2])
    end

    #do the left
    f2 = isomorphism(fuse(sp*space(mpo[1],3)'),sp*space(mpo[1],3)')
    @tensor output[1][-1 -2;-3 -4] := mpo[1][1,-2,2,-4]*conj(f2[-3,1,2])*utleg[-1]

    #do the right
    f2 = isomorphism(fuse(sp*space(mpo[len],1)),sp*space(mpo[len],1))
    @tensor output[end][-1 -2;-3 -4] := mpo[len][1,-2,2,-4]*f2[-1,2,1]*conj(utleg[-3])

    InfiniteMPO(output)
end


text:
include("../imports.jl")

@testset ExtendedTestSet "`StopOnNanLoss`" begin
    learner = testlearner(coeff = NaN, Recorder(), StopOnNaNLoss())
    # Epoch will be cancelled
    @test_nowarn fit!(learner, 1)
    @test learner.cbstate.history.epochs == 0
end


text:
###
### types
###

# TODO: should the various bits of free text / comments / timestamps be in this struct?
# Data can look like:
# 0,   100.00          / PSS/E-30.3    WED, SEP 15 2021  21:04
#    SE SNAPSHOT 09-15-2021 PEAK CASE 18:00
#    FULL COPY OF ETC.
"""
    $TYPEDEF

# Fields
$TYPEDFIELDS
"""
struct CaseID <: Tables.AbstractRow
    """
    IC Change code:
    0 - base case (i.e., clear the working case before adding data to it).
    1 - add data to the working case.
    """
    ic::Int
    "System base MVA."
    sbase::Float64
end

CaseID() = CaseID(0, 100.0)

Tables.columnnames(::CaseID) = fieldnames(CaseID)
Tables.getcolumn(cid::CaseID, i::Int) = getfield(cid, i)
Tables.getcolumn(cid::CaseID, nm::Symbol) = getfield(cid, nm)

# So all tabular data records (buses, loads, ...) can be handled the same.
abstract type Records end

# Create a instance of a `Records` subtype, with all fields (assumed to be Vector)
# containing `nrow` elements (initialised to undefined values, as they'll be overwritten).
(::Type{R})(nrow) where {R <: Records} = R(map(T -> T(undef, nrow), fieldtypes(R))...)

# Store data in column table so conversion to DataFrame efficient.
Tables.istable(::Type{<:Records}) = true
Tables.columnaccess(::Type{<:Records}) = true
Tables.columns(x::Records) = x
Tables.getcolumn(x::Records, i::Int) = getfield(x, i)
Tables.columnnames(x::R) where {R <: Records} = fieldnames(R)
Tables.schema(x::R) where {R <: Records} = Tables.Schema(fieldnames(R), fieldtypes(R))

"""
    $TYPEDEF

Each network bus to be represented in PSS/E is introduced through a bus data record.
Each bus data record includes not only data for the basic bus properties but also includes
information on an optionally connected shunt admittance to ground. That admittance can
represent a shunt capacitor or a shunt reactor (both with or without a real component) or a
shunt resistor. It must not represent line connected admittance, loads, line charging or
transformer magnetizing impedance, all of which are entered in other data categories.

# Fields
$TYPEDFIELDS
"""
struct Buses <: Records
    "Bus number (1 to 999997)."
    i::Vector{Int}
    """
    Alphanumeric identifier assigned to bus "I".
    The name may be up to twelve characters and must be enclosed in single quotes.
    NAME may contain any combination of blanks, uppercase letters, numbers and special characters, but the first character must not be a minus sign.
    """
    name::Vector{InlineString15}
    "Bus base voltage; entered in kV."
    basekv::Vector{Float64}
    """
    Bus type code:
    1 - load bus or other bus without any generator boundary condition.
    2 - generator or plant bus either regulating voltage or with a fixed reactive power (Mvar).
    A generator that reaches its reactive power limit will no longer control voltage but rather hold reactive power at its limit.
    3 - swing bus or slack bus.
    It has no power or reactive limits and regulates voltage at a fixed reference angle.
    4 - disconnected or isolated bus.
    """
    ide::Vector{Int}
    """
    Active component of shunt admittance to ground; entered in MW at one per unit voltage.
    GL should not include any resistive admittance load, which is entered as part of load data.
    """
    gl::Vector{Float64}
    """
    Reactive component of shunt admittance to ground; entered in Mvar at one per unit voltage.
    BL should not include any reactive impedance load, which is entered as part of load data;
    line charging and line connected shunts, which are entered as part of non-transformer branch data;
    or transformer magnetizing admittance, which is entered as part of transformer data.
    BL is positive for a capacitor, and negative for a reactor or an inductive load.
    """
    bl::Vector{Float64}
    "Area number. 1 through the maximum number of areas at the current size level."
    area::Vector{Int}
    "Zone number. 1 through the maximum number of zones at the current size level."
    zone::Vector{Int}
    "Bus voltage magnitude; entered in pu."
    vm::Vector{Float64}
    "Bus voltage phase angle; entered in degrees."
    va::Vector{Float64}
    "Owner number. 1 through the maximum number of owners at the current size level."
    owner::Vector{Int}
end

"""
    $TYPEDEF

Each network bus at which a load is to be represented must be specified in at least one load
data record. If multiple loads are to be represented at a bus, they must be individually
identified in a load data record for the bus with a different load identifier.
Each load at a bus can be a mixture of loads with different characteristics.

# Fields
$TYPEDFIELDS
"""
struct Loads <: Records
    "Buses number, or extended buses name enclosed in single quotes."
    i::Vector{Int}
    """
    One- or two-character uppercase non blank alphanumeric load identifier used to distinguish among multiple loads at bus "I".
    It is recommended that, at buses for which a single load is present, the load be designated as having the load identifier '1'.
    """
    id::Vector{InlineString3}  # TODO: confirm 3 is enough in practice, when whitespace can be present
    "Initial load status of one for in-service and zero for out-of-service."
    status::Vector{Bool}
    "Area to which the load is assigned (1 through the maximum number of areas at the current size level)."
    area::Vector{Int}
    "Zone to which the load is assigned (1 through the maximum number of zones at the current size level)."
    zone::Vector{Float64}
    "Active power component of constant MVA load; entered in MW."
    pl::Vector{Float64}
    "Reactive power component of constant MVA load; entered in Mvar."
    ql::Vector{Float64}
    "Active power component of constant current load; entered in MW at one per unit voltage."
    ip::Vector{Float64}
    "Reactive power component of constant current load; entered in Mvar at one per unit voltage."
    iq::Vector{Float64}
    "Active power component of constant admittance load; entered in MW at one per unit voltage."
    yp::Vector{Float64}
    """
    Reactive power component of constant admittance load; entered in Mvar at one per unit voltage.
    YQ is a negative quantity for an inductive load and positive for a capacitive load.
    """
    yq::Vector{Float64}
    "Owner to which the load is assigned (1 through the maximum number of owners at the current size level)."
    owner::Vector{Int}
end

"""
    $TYPEDEF

Each network bus to be represented as a generator or plant bus in PSS/E must be specified
in a generator data record. In particular, each bus specified in the bus data input with a
type code of two (2) or three (3) must have a generator data record entered for it.

# Fields
$TYPEDFIELDS
"""
struct Generators <: Records
    "Bus number, or extended bus name enclosed in single quotes."
    i::Vector{Int}
    """
    One- or two-character uppercase non blank alphanumeric machine identifier used to distinguish among multiple machines at bus "I".
    It is recommended that, at buses for which a single machine is present, the machine be designated as having the machine identifier ’1’.
    ID = ’1’ by default.
    """
    id::Vector{InlineString3}  # TODO: confirm 3 is enough in practice, when whitespace can be present
    "Generator active power output; entered in MW. PG = 0.0 by default."
    pg::Vector{Float64}
    """
    Generator reactive power output; entered in Mvar.
    QG needs to be entered only if the case, as read in, is to be treated as a solved case.
    QG = 0.0 by default.
    """
    qg::Vector{Float64}
    """
    Maximum generator reactive power output; entered in Mvar.
    For fixed output gen- erators (i.e., nonregulating), QT must be equal to the fixed Mvar output.
    QT = 9999.0 by default.
    """
    qt::Vector{Float64}
    """
    Minimum generator reactive power output; entered in Mvar.
    For fixed output generators, QB must be equal to the fixed Mvar output.
    QB = -9999.0 by default.
    """
    qb::Vector{Float64}
    "Regulated voltage setpoint; entered in pu. VS = 1.0 by default."
    vs::Vector{Float64}
    """
    Bus number, or extended bus name enclosed in single quotes,
    of a remote type 1 or 2 bus whose voltage is to be regulated by this plant to the value specified by VS.
    If bus IREG is other than a type 1 or 2 bus, bus "I" regulates its own voltage to the value specified by VS.
    IREG is entered as zero if the plant is to regulate its own voltage and must be zero for a type three (swing) bus.
    IREG = 0 by default.
    """
    ireg::Vector{Int}
    """
    Total MVA base of the units represented by this machine; entered in MVA.
    This quantity is not needed in normal power flow and equivalent construction work,
    but is required for switching studies, fault analysis, and dynamic simulation.
    MBASE = system base MVA by default.
    """
    mbase::Vector{Float64}
    """
    Complex machine impedance, ZSORCE; entered in pu on MBASE base.
    This data is not needed in normal power flow and equivalent construction work,
    but is required for switching studies, fault analysis, and dynamic simulation.
    For dynamic simulation, this impedance must be set equal to the unsaturated subtransient impedance for those generators to be modeled by subtransient level machine models,
    and to unsaturated transient impedance for those to be modeled by classical or transient level models.
    For short-circuit studies, the saturated subtransient or transient impedance should be used.
    ZR = 0.0 by default.
    """
    zr::Vector{Float64}
    "See `zr`. ZX = 1.0 by default."
    zx::Vector{Float64}
    """
    Step-up transformer impedance, XTRAN; entered in pu on MBASE base.
    XTRAN should be entered as zero if the step-up transformer is explicitly modeled as a network branch and bus "I" is the terminal bus.
    RT+jXT = 0.0 by default.
    """
    rt::Vector{Float64}
    "See `rt`. RT+jXT = 0.0 by default."
    xt::Vector{Float64}
    """
    Step-up transformer off-nominal turns ratio; entered in pu.
    GTAP is used only if XTRAN is nonzero.
    GTAP = 1.0 by default.
    """
    gtap::Vector{Float64}
    """
    Initial machine status of one for in-service and zero for out-of-service.
    STAT = 1 by default.
    """
    stat::Vector{Bool}
    """
    Percent of the total Mvar required to hold the voltage at the bus controlled by this bus "I" that are to be contributed by the generation at bus "I";
    RMPCT must be positive.
    RMPCT is needed if IREG specifies a valid remote bus and there is more than one local or remote voltage controlling device
    (plant, switched shunt, FACTS device shunt element, or VSC dc line converter) controlling the voltage at bus IREG to a setpoint.
    RMPCT is needed also if bus "I" itself is being controlled locally or remotely by one or more other setpoint mode voltage controlling devices.
    RMPCT = 100.0 by default.
    """
    rmpct::Vector{Float64}
    "Maximum generator active power output; entered in MW. PT = 9999.0 by default."
    pt::Vector{Float64}
    "Minimum generator active power output; entered in MW. PB = -9999.0 by default."
    pb::Vector{Float64}
    """
    Owner number; (1 through the maximum number of owners at the current size level).
    Each machine may have up to four owners.
    By default, O1 is the owner to which bus "I" is assigned and O2, O3, and O4 are zero.
    """
    oi::Vector{Int}
    """
    Fraction of total ownership assigned to owner Oi; each Fi must be positive.
    The Fi values are normalized such that they sum to 1.0 before they are placed in the working case.
    By default, each Fi is 1.0.
    """
    fi::Vector{Float64}
end

"""
    $TYPEDEF

In PSS/E, the basic transmission line model is an Equivalent Pi connected between network buses.

Data for shunt equipment units, such as reactors, which are connected to and switched with the line,
are entered in the same data record.

!!! note "Shunts connected to buses"
    To represent shunts connected to buses, that shunt data should be entered in the bus data record.

!!! note "Transformers"
    Branches to be modeled as transformers are not specified in this data category;
    rather, they are specified in the transformer data category.

# Fields
$TYPEDFIELDS
"""
struct Branches <: Records
    "Branch \"from bus\" number, or extended bus name enclosed in single quotes."
    i::Vector{Int}
    """
    Branch "to bus" number, or extended bus name enclosed in single quotes.
    "J" is entered as a negative number, or with a minus sign before the first character of the extended bus name,
    to designate it as the metered end; otherwise, bus "I" is assumed to be the metered end.
    """
    j::Vector{Int}
    """
    One- or two-character uppercase nonblank alphanumeric branch circuit identifier;
    the first character of CKT must not be an ampersand "&".
    It is recommended that single circuit branches be designated as having the circuit identifier '1'.
    CKT = '1' by default.
    """
    ckt::Vector{InlineString3}
    "Branch resistance; entered in pu. A value of R must be entered for each branch."
    r::Vector{Float64}
    "Branch reactance; entered in pu. A nonzero value of X must be entered for each branch."
    x::Vector{Float64}
    "Total branch charging susceptance; entered in pu. B = 0.0 by default."
    b::Vector{Float64}
    """
    First loading rating; entered in MVA.
    If RATEA is set to 0.0, the default value, this branch will not be included in any examination of circuit loading.

    Ratings are entered as:
    ``MVA_{rated} = sqrt(3) × E_{base} × I_{rated} × 10^{-6}`` where:
    - ``E_{base}`` is the base line-to-line voltage in volts of the buses to which the terminal of the branch is connected.
    - ``I_{rated}`` is the branch rated phase current in amperes.
    """
    rate_a::Vector{Float64}
    "Second loading rating; entered in MVA. RATEB = 0.0 by default."
    rate_b::Vector{Float64}
    "Third loading rating; entered in MVA. RATEC = 0.0 by default."
    rate_c::Vector{Float64}
    """
    Complex admittance of the line shunt at the bus "I" end of the branch; entered in pu.
    BI is negative for a line connected reactor and positive for line connected capacitor.
    GI + jBI = 0.0 by default.
    """
    gi::Vector{Float64}
    """
    Complex admittance of the line shunt at the bus "I" end of the branch; entered in pu.
    BI is negative for a line connected reactor and positive for line connected capacitor.
    GI + jBI = 0.0 by default.
    """
    bi::Vector{Float64}
    """
    Complex admittance of the line shunt at the bus "J" end of the branch; entered in pu.
    BJ is negative for a line connected reactor and positive for line connected capacitor.
    GJ + jBJ = 0.0 by default.
    """
    gj::Vector{Float64}
    """
    Complex admittance of the line shunt at the bus "J" end of the branch; entered in pu.
    BJ is negative for a line connected reactor and positive for line connected capacitor.
    GJ + jBJ = 0.0 by default.
    """
    bj::Vector{Float64}
    """
    Initial branch status where 1 designates in-service and 0 designates out-of-service.
    ST = 1 by default.
    """
    st::Vector{Bool}
    "Line length; entered in user-selected units. LEN = 0.0 by default."
    len::Vector{Float64}
    """
    Owner number; 1 through the maximum number of owners at the current size level.
    Each branch may have up to four owners.
    By default, O1 is the owner to which bus "I" is assigned and O2, O3, and O4 are zero.
    """
    oi::Vector{Int}
    """
    Fraction of total ownership assigned to owner Oi; each Fi must be positive.
    The Fi values are normalized such that they sum to 1.0 before they are placed in the working case.
    By default, each Fi is 1.0.
    """
    fi::Vector{Float64}
end

"""
    Network

Representation of a power network.

The PSS/E data format comprises 16 data categories of network and equipment
elements, each of which requires a particular type of data.

Similarly, a `Network` stores the data from each category in its own dedicated structure.

Currently supported are:
1. [`CaseID`](@ref)
1. [`Buses`](@ref)
1. [`Loads`](@ref)
1. [`Generators`](@ref)
1. [`Branches`](@ref)

`CaseID` data is a single row (in the Tables.jl-sense).
You can access it like `network.caseid` and interact with it like a `NamedTuple`,
or even convert it to a `NamedTuple` with `NamedTuple(caseid)`.

All other records (buses, loads, etc.) can be accessed also via the fields, for example
`network.buses`, and each is returned as lightweight table structure (again, in the Tables.jl-sense).
That is, all structures implement the Tables.jl interface, so can be passed to any valid
sink, such as a `DataFrame` like `DataFrame(network.buses)`.

For more info on working with tables see [Tables.jl](https://tables.juliadata.org/), and for
common table operations see [TableOperations.jl](https://github.com/JuliaData/TableOperations.jl).

# Fields
$TYPEDFIELDS
"""
struct Network
    "Case identification data."
    caseid::CaseID
    "Bus records."
    buses::Buses
    "Load records."
    loads::Loads
    "Generator records."
    generators::Generators
    "Non-transformer Branch records."
    branches::Branches
end


text:
# Generate IV regression data
begin
      n = 3000;
      sig = [1 0.5 0.5;
             0 0.5 0 ;
             0 0  1];

      mu = [1,9,0];
      mn = MvNormal(mu,sig * sig')
      W = rand(mn,n)'
      x = W[:,1];z=W[:,2];ϵ=W[:,3];
      y = x * 3 + ϵ;
      β_OLS = inv(x'*x)*(x'*y);
      β_IV = inv(z'*x)*(z'*y);
      w = hcat(x,y,z);
end

function apply_row(func::Function, w::Array{Float64,2})
      δ = zeros(size(w,1));
      for (index, value) in enumerate(w[:,1])
           δ[index] = func(w[index,:]);
      end
      return(δ);
end

obj = θ->(begin
      η = (w[:,1] .* θ - w[:,2]) .* w[:,3];
      η'*η
      end)

using Optim
using Optim: converged, maximum, maximizer, minimizer, iterations #some extra functions
begin "Test Optim"
      result = optimize(obj,zeros(1),Newton(),autodiff = :forward)
      # @show converged(result) || error("Failed to converge in $(iterations(result)) iterations")
      @show xmin = result.minimizer
      @show result.minimum
end


text:
# Methods to resolve ambiguities with `Diagonal`
@inline *(tvec::TransposedVector, diagonal::Diagonal) = transpose(diagonal * transpose(tvec))
@inline *(::Diagonal, ::TransposedVector) = error("Cannot right-multiply matrix by transposed vector")

@inline A_mul_Bt(d::Diagonal, tvec::TransposedVector) = d*transpose(tvec)

@inline At_mul_B(tvec::TransposedVector, d::Diagonal) = error("Cannot left-multiply matrix by vector")

@inline A_mul_Bc(d::Diagonal, tvec::TransposedVector) = d*ctranspose(tvec)

@inline Ac_mul_B(tvec::TransposedVector, d::Diagonal) = error("Cannot left-multiply matrix by vector")

@inline \(::Diagonal, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline \(::Bidiagonal, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline \{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::TransposedVector{TB}) = error("Cannot left-divide matrix by transposed vector")

@inline At_ldiv_B(::Bidiagonal, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline At_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::TransposedVector{TB}) = error("Cannot left-divide matrix by transposed vector")

@inline Ac_ldiv_B(::Bidiagonal, ::TransposedVector) = error("Cannot left-divide matrix by transposed vector")
@inline Ac_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::TransposedVector{TB}) = error("Cannot left-divide matrix by transposed vector")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Lifecycle describes actions that the management system should take in response to container lifecycle events. For the PostStart and PreStop lifecycle handlers, management of the container blocks until the action is complete, unless the container process fails, in which case the handler is aborted.

    IoK8sApiCoreV1Lifecycle(;
        postStart=nothing,
        preStop=nothing,
    )

    - postStart::IoK8sApiCoreV1Handler : PostStart is called immediately after a container is created. If the handler fails, the container is terminated and restarted according to its restart policy. Other management of the container blocks until the hook completes. More info: https://kubernetes.io/docs/concepts/containers/container-lifecycle-hooks/#container-hooks
    - preStop::IoK8sApiCoreV1Handler : PreStop is called immediately before a container is terminated due to an API request or management event such as liveness/startup probe failure, preemption, resource contention, etc. The handler is not called if the container crashes or exits. The reason for termination is passed to the handler. The Pod's termination grace period countdown begins before the PreStop hooked is executed. Regardless of the outcome of the handler, the container will eventually terminate within the Pod's termination grace period. Other management of the container blocks until the hook completes or until the termination grace period is reached. More info: https://kubernetes.io/docs/concepts/containers/container-lifecycle-hooks/#container-hooks
"""
mutable struct IoK8sApiCoreV1Lifecycle <: SwaggerModel
    postStart::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Handler } # spec name: postStart
    preStop::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Handler } # spec name: preStop

    function IoK8sApiCoreV1Lifecycle(;postStart=nothing, preStop=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Lifecycle, Symbol("postStart"), postStart)
        setfield!(o, Symbol("postStart"), postStart)
        validate_property(IoK8sApiCoreV1Lifecycle, Symbol("preStop"), preStop)
        setfield!(o, Symbol("preStop"), preStop)
        o
    end
end # type IoK8sApiCoreV1Lifecycle

const _property_map_IoK8sApiCoreV1Lifecycle = Dict{Symbol,Symbol}(Symbol("postStart")=>Symbol("postStart"), Symbol("preStop")=>Symbol("preStop"))
const _property_types_IoK8sApiCoreV1Lifecycle = Dict{Symbol,String}(Symbol("postStart")=>"IoK8sApiCoreV1Handler", Symbol("preStop")=>"IoK8sApiCoreV1Handler")
Base.propertynames(::Type{ IoK8sApiCoreV1Lifecycle }) = collect(keys(_property_map_IoK8sApiCoreV1Lifecycle))
Swagger.property_type(::Type{ IoK8sApiCoreV1Lifecycle }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Lifecycle[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Lifecycle }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Lifecycle[property_name]

function check_required(o::IoK8sApiCoreV1Lifecycle)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Lifecycle }, name::Symbol, val)
end


text:
using ValidatedNumerics
using FactCheck

facts("Complex interval operations") do
    a = @interval 1im
    @fact a -->  Interval(0) + Interval(1)*im
    @fact a * a --> Interval(-1)
    @fact a + a --> Interval(2)*im
    @fact a - a --> 0
    @fact a / a --> 1
    @fact a^2 --> -1
end


text:
include(normpath(joinpath(@__DIR__, "..", "..", "src", "models", "EMOM", "EMOM.jl")))
using DataStructures
using NCDatasets
using MPI
using Formatting
using ArgParse, JSON

println("""
This program generates initial file (technically a restart file) for
IOM to start.
""")

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--init-profile-TEMP"
            help = "Initial ocean profile. It must contains: TEMP, SALT, HMXL."
            arg_type = String
            required = true

        "--init-profile-SALT"
            help = "Initial ocean profile. It must contains: TEMP, SALT, HMXL."
            arg_type = String
            required = true

        "--init-profile-HMXL"
            help = "Initial ocean profile. It must contains: TEMP, SALT, HMXL."
            arg_type = String
            required = true

        "--config-file"
            help = "config TOML file"
            arg_type = String
            required = true

    end

    return parse_args(s)
end

parsed = DataStructures.OrderedDict(parse_commandline())
JSON.print(parsed,4)

MPI.Init()

println("Processing data...")

using TOML
config = TOML.parsefile(parsed["config-file"])

domain_file = config["MODEL_CORE"]["domain_file"]

init_file  = config["MODEL_MISC"]["init_file"]
topo_file  = config["MODEL_CORE"]["topo_file"]

Nz = length(config["MODEL_CORE"]["z_w"]) - 1 # Layers used. Thickness ≈ 503m

Dataset(domain_file, "r") do ds
    global Nx = ds.dim["ni"]
    global Ny = ds.dim["nj"]
    global mask = ds["mask"][:]
end

Dataset(parsed["init-profile-TEMP"], "r") do ds
    global TEMP = permutedims(nomissing(ds["TEMP"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

Dataset(parsed["init-profile-SALT"], "r") do ds
    global SALT = permutedims(nomissing(ds["SALT"][:, :, 1:Nz, 1],  NaN), [3, 1, 2])
end

Dataset(parsed["init-profile-HMXL"], "r") do ds
    global HMXL  = nomissing(ds["HMXL"][:, :, 1],  0.0)
end

valid_idx = isfinite.(TEMP)

Dataset(topo_file, "r") do ds
    global Nz_bot = ds["Nz_bot"][:]
end

println("Check if all the valid grid will be assigned a non-NaN value.")
for i=1:Nx, j=1:Ny
    local Nz = Nz_bot[i, j]

    if !isfinite(HMXL[i, j])
        throw(ErrorException("Error: Grid ($i, $j) of HMXL is assigned an NaN."))
    end

    for k=1:Nz
        if !isfinite(TEMP[k, i, j])
            throw(ErrorException("Error: Grid ($k, $i, $j) of TEMP is assigned an NaN."))
        end

        if !isfinite(SALT[k, i, j])
            throw(ErrorException("Error: Grid ($k, $i, $j) of SALT is assigned an NaN."))
        end

    end
end
println("NaNs are consistent.")

println("Create a model to save initial condition.")

ev = EMOM.Env(config["MODEL_CORE"])
mb = EMOM.ModelBlock(ev; init_core=false)

mb.fi.sv[:TEMP][valid_idx] .= TEMP[valid_idx]
mb.fi.sv[:SALT][valid_idx] .= SALT[valid_idx]

println(size(mb.fi.HMXL))
println(size(HMXL))
mb.fi.HMXL[:] = HMXL

println(format("Output file: {}.", init_file))

EMOM.takeSnapshot(DateTimeNoLeap(1,1,1), mb, init_file)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using NonconvexPavito, LinearAlgebra, Test

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

@testset "First order - $first_order" for first_order in [true, false]
    options = PavitoIpoptCbcOptions(first_order = first_order)
    @testset "Simple constraints" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = PavitoIpoptCbcAlg()

        setinteger!(m, 1, true)
        r1 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r1.minimizer[1] - round(Int, r1.minimizer[1]) ≈ 0 atol = 1e-7

        setinteger!(m, 1, false)
        setinteger!(m, 2, true)
        r2 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r2.minimizer[2] - round(Int, r2.minimizer[2]) ≈ 0 atol = 1e-7
    end

    @testset "Block constraints" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, FunctionWrapper(x -> [g(x, 2, 0), g(x, -1, 1)], 2))

        alg = PavitoIpoptCbcAlg()

        setinteger!(m, 1, true)
        r1 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r1.minimizer[1] - round(Int, r1.minimizer[1]) ≈ 0 atol = 1e-7

        setinteger!(m, 1, false)
        setinteger!(m, 2, true)
        r2 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r2.minimizer[2] - round(Int, r2.minimizer[2]) ≈ 0 atol = 1e-7
    end

    @testset "Infinite bounds" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [Inf, Inf])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = PavitoIpoptCbcAlg()

        setinteger!(m, 1, true)
        r1 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r1.minimizer[1] - round(Int, r1.minimizer[1]) ≈ 0 atol = 1e-7

        setinteger!(m, 1, false)
        setinteger!(m, 2, true)
        r2 = optimize(m, alg, [1.234, 2.345], options = options)
        @test r2.minimizer[2] - round(Int, r2.minimizer[2]) ≈ 0 atol = 1e-7
    end
end


text:
##################
# Initialization #
##################

initialize() = @checked API.AMGX_initialize()
initialize_plugins() = @checked API.AMGX_initialize_plugins()
finalize() = @checked API.AMGX_finalize()
finalize_plugins() = @checked API.AMGX_finalize_plugins()

##############
# Versioning #
##############

function api_version()
    major, minor = Ref{Cint}(), Ref{Cint}()
    @checked API.AMGX_get_api_version(major, minor)
    VersionNumber(major[], minor[])
end

function build_info()
    ref_version = Ref{Cstring}()
    ref_date = Ref{Cstring}()
    ref_time = Ref{Cstring}()
    @checked API.AMGX_get_build_info_strings(ref_version, ref_date, ref_time)
    return unsafe_string(ref_version[]), unsafe_string(ref_date[]), unsafe_string(ref_time[])
end

function versioninfo(io::IO=stdout)
    version, date, time = build_info()
    println(io, "AMGX version $(version)")
    println(io, "Built on $(date), $(time)")
    print(io, "API version $(api_version())")
end

###########
# Pinning #
###########

function pin_memory(v::Vector)
    GC.@preserve v begin
        pin_memory(pointer(v), sizeof(v))
    end
end

function pin_memory(ptr::Ptr, n_bytes::Int)
    @checked API.AMGX_pin_memory(ptr, n_bytes)
end

function unpin_memory(v::Vector)
    GC.@preserve v begin
        unpin_memory(pointer(v))
    end
end

function unpin_memory(ptr::Ptr)
    @checked API.AMGX_unpin_memory(ptr)
end

##################
# Print callback #
##################

# This can be set before initializing the library
function register_print_callback(f)
    run_f(str::Cstring, _::Cint) = f(unsafe_string(str))
    f_cfunc = @cfunction($run_f, Cvoid, (Cstring, Cint))
    @checked API.AMGX_register_print_callback(f_cfunc)
end

##################
# Signal handler #
##################

install_signal_handler() = @checked API.AMGX_install_signal_handler()
reset_signal_handler() = @checked API.AMGX_reset_signal_handler()


text:
module test_julia_iterate_counter

using Test

mutable struct G
     counter
end

function Base.iterate(g::G)
    (g.counter, 1)
end

function Base.iterate(g::G, val)
    if g.counter > 2
        return nothing
    end
    g.counter += val
    (g.counter, 1)
end

g = G(0)
for x in g
end
@test g.counter == 3

end # module test_julia_iterate_counter

module test_julia_iterate_basic

using Test

struct A
end

function Base.iterate(::A)
    (:value, :next_state)
end

entered = []
function Base.iterate(::A, state)
    push!(entered, state)
    nothing
end

a = A()
stack = []
for x in a
    push!(stack, x)
end
@test [:value] == stack
@test [:next_state] == entered

function Base.length(::A)
    2
end
vals = collect(a)
@test :value == vals[1]
@test_throws UndefRefError vals[2]

end # module test_julia_iterate_basic

module test_julia_iterate_dict

using Test

struct G
    d::Dict
end

struct V
    val
end

function Base.iterate(g::G, state=firstindex(g.d.vals))
    y = iterate(g.d, state)
    y === nothing && return nothing
    (V(y[1]), y[2])
end

g = G(Dict(:a => 1))
@test first(g) == V(:a => 1)

end # module test_julia_iterate_dict


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

name = "liblinear"
version = v"2.30.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/cjlin1/liblinear/archive/v230.tar.gz", "9b57710078206d4dbbe75e9015d4cf7fabe4464013fe0e89b8a2fe40038f8f51"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd liblinear-230/
mkdir -p ${prefix}/bin
mkdir -p ${prefix}/lib
if [[ "${target}" == *-freebsd* ]] || [[ "${target}" == *-apple-* ]]; then
    CC=gcc
    CXX=g++
fi
make
make lib
cp train${exeext} ${bindir}
cp predict${exeext} ${bindir}
cp liblinear.${dlext} ${libdir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("liblinear", :liblinear),
    ExecutableProduct("train", :train),
    ExecutableProduct("predict", :predict)
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Gnuplot
U  = LinRange(0,10π, 80)
V = LinRange(0,2π, 20)
x = [(1-0.1*cos(v))*cos(u) for u in U, v in V]
y = [(1-0.1*cos(v))*sin(u) for u in U, v in V]
z = [0.1*(sin(v) + u/1.7 - 10) for u in U, v in V]
@gsp "set pm3d depthorder" :-
@gsp :- "set style fill transparent solid 0.7" "set view equal xyz"
@gsp :- x y z "w pm3d notit" "set xyplane -0.05" palette(:ice)
save(term="pngcairo size 600,400", output="tubes2.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
println("Hello")
display("Julia")

# %%
println("Hello"); flush(stdout)
display("Julia")

# %%
show("Hello")
display("Julia")

# %%
show("Hello"); flush(stdout)
display("Julia")

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using ElectronicStructure
using ElectronicStructure.Utility

using SpecialFunctions: erf

type PsPot_HGH
  itype::Int
  atsymb::String
  zval::Float64
  lloc::Int
  lmax::Int
  rloc::Float64
  rc::Array{Float64}   # indexed (l+1), l=0,1,2,3
  c::Array{Float64}    # indexed 1,2,3,4
  h::Array{Float64,3}  # originally indexed [0:3,1:3,1:3]
  k::Array{Float64,3}  # indexed [0:3,1:3,1:3]
  nprj::Array{Int}
  snprj::Int           # snprj = sum(nprj)
  lll::Array{Int}
  ipr::Array{Int}
end

# Constructor
function PsPot_HGH( itype::Int, atsymb::String, filename::String )

  @printf("\nFile = %s\n", filename)
  file = open( filename )

  comment = readline(file)

  lines = split( readline(file) )
  zval = parse( Float64, lines[1] )
  zion = parse( Float64, lines[2] )

  lines  = split( readline(file) )
  pspcod = parse( Int, lines[1] )
  pspxc  = parse( Int, lines[2] )
  lmax = parse( Int, lines[3] )
  lloc = parse( Int, lines[4] )
  mmax   = parse( Int, lines[5] )
  r2well = parse( Int, lines[6] )

  rc = zeros(Float64,4)
  c = zeros(Float64,4)
  h = zeros(Float64,4,3,3)
  k = zeros(Float64,4,3,3)

  lines = split( readline(file) )
  rloc = parse( Float64, lines[1] )
  c[1] = parse( Float64, lines[2] )
  c[2] = parse( Float64, lines[3] )
  c[3] = parse( Float64, lines[4] )
  c[4] = parse( Float64, lines[5] )

  const ANGMOM = ["s", "p", "d", "f"]

  l = 0  # s
  lines = split( readline(file) )
  #
  rc[1] = parse( lines[1] )
  #
  h[1,1,1] = parse( lines[2] )
  h[1,2,2] = parse( lines[3] )
  h[1,3,3] = parse( lines[4] )

  for l = 1:3
    lines = split( readline(file) )
    #
    rc[l+1] = parse( lines[1] )
    #
    h[l+1,1,1] = parse( lines[2] )
    h[l+1,2,2] = parse( lines[3] )
    h[l+1,3,3] = parse( lines[4] )
    #
    lines = split( readline(file) )
    k[l+1,1,1] = parse( lines[1] )
    k[l+1,2,2] = parse( lines[2] )
    k[l+1,3,3] = parse( lines[3] )
    #
  end

  # additional info
  lines = split(readline(file))
  nprj = zeros(Int,4)
  nprj[1] = parse( Int, lines[1] )
  nprj[2] = parse( Int, lines[2] )
  nprj[3] = parse( Int, lines[3] )
  nprj[4] = parse( Int, lines[4] )

  snprj  = sum(nprj)
  lll = zeros(Int,snprj)
  ipr = zeros(Int,snprj)

  iv = 0
  for il = 0:lmax
    for ii = 1:nprj[il+1]
      iv = iv + 1
      lll[iv] = il
      ipr[iv] = ii
    end
  end

  close(file)

  const M_HALF = 0.5
  const M_THREE = 3.0
  const M_FIVE = 5.0
  const M_ONE = 1.0

  # from Octopus code, see also the appendix of HGH paper

  h[0+1, 1, 2] = -M_HALF    * sqrt(M_THREE/M_FIVE) * h[0+1, 2, 2]
  h[0+1, 1, 3] =  M_HALF    * sqrt(M_FIVE/21.0)    * h[0+1, 3, 3]
  h[0+1, 2, 3] = -M_HALF    * sqrt(100.0/63.0)     * h[0+1, 3, 3]
  h[1+1, 1, 2] = -M_HALF    * sqrt(M_FIVE/7.0)     * h[1+1, 2, 2]
  h[1+1, 1, 3] =  M_ONE/6.0 * sqrt(35.0/11.0)      * h[1+1, 3, 3]
  h[1+1, 2, 3] = -M_ONE/6.0 * ( 14.0 / sqrt(11.0)) * h[1+1, 3, 3]
  h[2+1, 1, 2] = -M_HALF    * sqrt(7.0/9.0)        * h[2+1, 2, 2]
  h[2+1, 1, 3] =  M_HALF    * sqrt(63.0/143.0)     * h[2+1, 3, 3]
  h[2+1, 2, 3] = -M_HALF    * (18.0/sqrt(143.0))   * h[2+1, 3, 3]

  k[0+1, 1, 2] = -M_HALF    * sqrt(M_THREE/M_FIVE) * k[0+1, 2, 2]
  k[0+1, 1, 3] =  M_HALF    * sqrt(M_FIVE/21.0)    * k[0+1, 3, 3]
  k[0+1, 2, 3] = -M_HALF    * sqrt(100.0/63.0)     * k[0+1, 3, 3]
  k[1+1, 1, 2] = -M_HALF    * sqrt(M_FIVE/7.0)     * k[1+1, 2, 2]
  k[1+1, 1, 3] =  M_ONE/6.0 * sqrt(35.0/11.0)      * k[1+1, 3, 3]
  k[1+1, 2, 3] = -M_ONE/6.0 * (14.0 / sqrt(11.0))  * k[1+1, 3, 3]
  k[2+1, 1, 2] = -M_HALF    * sqrt(7.0/9.0)        * k[2+1, 2, 2]
  k[2+1, 1, 3] =  M_HALF    * sqrt(63.0/143.0)     * k[2+1, 3, 3]
  k[2+1, 2, 3] = -M_HALF    * (18.0/sqrt(143.0))   * k[2+1, 3, 3]

  # Parameters are symmetric.
  # be careful of index !!!
  for ik = 1:4
    for i = 1:3
      for j = i+1:3
        h[ik, j, i] = h[ik, i, j]
        k[ik, j, i] = k[ik, i, j]
      end
    end
  end

  psp = PsPot_HGH( itype, atsymb, zval, lloc, lmax, rloc,
                   rc, c, h, k, nprj, snprj, lll, ipr )

  return psp
end

# Display information about HGH pseudopotential
function info_PsPot_HGH( psp::PsPot_HGH )

  const ANGMOM = ["s", "p", "d", "f"]

  @printf("\nLocal pseudopotential info:\n\n")
  @printf("rloc: %f, c: %f, %f, %f, %f\n", psp.rloc, psp.c[1], psp.c[2], psp.c[3], psp.c[4])
  @printf("\n")
  @printf("Nonlocal pseudopotential info:\n\n")
  for i=1:4
    @printf("Angular momentum: %s, rc = %f\n", ANGMOM[i], psp.rc[i])
    @printf("h = \n")
    PrintMatrix( reshape(psp.h[i,:,:],(3,3) ) )
    @printf("\n")
  end

end

"""
Evaluate HGH local pseudopotential in R-space
"""
function HGH_eval_Vloc_R( psp, r::Array{Float64,2} )

  Npoints = size(r)[1]
  Vloc = zeros(Npoints)

  term1 = psp.c[1]
  for ip = 1:Npoints
    rrloc = r[ip]/psp.rloc
    for i = 2:4
      term1 = term1 + psp.c[i]*rrloc^(2.0*(i-1))
    end
    Vloc[ip] = -psp.zval/r[ip] * erf( rrloc/sqrt(2.0) ) + exp(-0.5*rrloc^2)*term1
  end
  return Vloc
end

"""
Evaluate HGH local pseudopotential in G-space
"""
function HGH_eval_Vloc_G( psp, G2, Ω )

  Ng = size(G2)[1]
  Vg = zeros(Ng)

  rloc = psp.rloc
  zval = psp.zval
  c1 = psp.c[1]
  c2 = psp.c[2]
  c3 = psp.c[3]
  c4 = psp.c[4]

  pre1 = -4*pi*zval/Ω
  pre2 = sqrt(8*pi^3)*rloc^3/Ω
  #
  for ig=2:Ng
    Gr = sqrt(G2[ig])*rloc
    expGr2 = exp(-0.5*Gr^2)
    Vg[ig] = pre1/G2[ig]*expGr2 + pre2*expGr2 * (c1 + c2*(3-Gr^2) +
             c3*(15 - 10*Gr^2 + Gr^4) + c4*(105 - 105*Gr^2 + 21*Gr^4 - Gr^6) )
  end
  # limiting value, with minus sign ?
  Vg[1] = 2*pi*zval*rloc^2 + (2*pi)^1.5 * rloc^3 * (c1 + 3.0*c2 + 15*c3 + 105*c4)

  return Vg
end

"""
Evaluate HGH projector function in R-space.
"""
function HGH_eval_proj_R( psp, l, i, r::Float64 )
  x = sqrt( gamma( l + (4*i-1)/2.0 ) )
  if l==0 & i==1
    rr = 1.0
  else
    rr = r^(l + 2*(i-1))
  end
  fprj = sqrt(2) * rr * exp(-r^2/(2*psp.rc[l+1])) /
         ( p.rc[l+1]^(l + (4*i-1)/2) * x )
  return fprj
end

"""
Evaluate HGH projector function in G-space.
"""
function HGH_eval_proj_G( psp, l, iproj, Gm, Ω )

  # NOTICE that G is the magnitudes of G-vectors
  Ng = size(Gm)[1]

  Vprj = zeros(Ng)

  rrl = psp.rc[l+1]

  # s-channel
  if l == 0

    if iproj==1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = exp( -0.5*Gr2 )
      end

    elseif iproj==2

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = 2.0/sqrt(15.0) * exp( -0.5*Gr2 ) * ( 3.0 - Gr2 )
      end

    elseif iproj==3

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (4.0/3.0)/sqrt(105.0) * exp( -0.5*Gr2 ) * (15.0 - 10.*Gr2 + Gr2^2)
       end

    end  # if iproj

  # p-channel
  elseif l == 1

    if iproj == 1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (1.0/sqrt(3.0)) * exp(-0.5*Gr2) * Gm[ig]
      end

    elseif iproj == 2

      for ig = 1:Ng
        Gr2 = (Gm[ig]*rrl)^2
        Vprj[ig] = (2./sqrt(105.)) * exp(-0.5*Gr2) * Gm[ig]*(5. - Gr2)
      end

    elseif iproj == 3

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl)^2
        Vprj[ig] = (4./3.)/sqrt(1155.) * exp(-0.5*Gr2) * Gm[ig] * (35. - 14.*Gr2 + Gr2^2)
      end

    end # if iproj

  # d-channel
  elseif l == 2

    if iproj == 1

      for ig = 1:Ng
        Gr2 = ( Gm[ig]*rrl )^2
        Vprj[ig] = (1.0/sqrt(15.0)) * exp(-0.5*Gr2) * Gm[ig]^2
      end

    elseif iproj == 2

      for ig = 1:Ng
        Gr2 = (Gm[ig]*rrl)^2
        Vprj[ig] = (2./3.)/sqrt(105.) * exp(-0.5*Gr2) * Gm[ig]^2 * (7.-Gr2)
      end

    end # if iproj

  # f-channel
  elseif l == 3

    for ig = 1:Ng
      Gr2 = ( Gm[ig]*rrl )^2
      Vprj[ig] = Gm[ig]^3 * exp(-0.5*Gr2)
    end

  end  # if l

  pre =  4 * pi^(5./4.) * sqrt( 2.^(l+1) * rrl^(2*l+3) / Ω )

  Vprj[:] = pre * Vprj[:]

  return Vprj

end


text:
using Dates

include("../structures.jl")

function loadPalette(theme)::Palette
    local palette = PlotThemes.palette(theme) # get_color_palette(:auto, plot_color(:white), 17)
    return Palette(palette, 1)
end

function nextColor!(p::Palette)
    p.index = (p.index + 1) % length(p.palette) + 1
    return p.palette[p.index]
end

function defaultStyles()::Dict{String,Any}
    return Dict{String,Any}(
        "color" => nothing,
        "fillcolor" => nothing,
        "fillalpha" => 0.2,
        "fillrange" => nothing,
        "offset" => 0.0,
    )
end

using Formatting
global Y_FORMATTER = yi -> begin
    if yi > 100
        yi = round(yi)
    end
    replace(format(round(yi, digits=3), commas = true), "," => " ")
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
export A2CGAELearner

using Flux

"""
    A2CGAELearner(;kwargs...)
# Keyword arguments
- `approximator`, an [`ActorCritic`](@ref) based [`NeuralNetworkApproximator`](@ref)
- `γ::Float32`, reward discount rate.
- 'λ::Float32', lambda for GAE-lambda
- `actor_loss_weight::Float32`
- `critic_loss_weight::Float32`
- `entropy_loss_weight::Float32`
"""
Base.@kwdef mutable struct A2CGAELearner{A<:ActorCritic} <: AbstractLearner
    approximator::A
    γ::Float32
    λ::Float32
    max_grad_norm::Union{Nothing,Float32} = nothing
    norm::Float32 = 0.f0
    actor_loss_weight::Float32
    critic_loss_weight::Float32
    entropy_loss_weight::Float32
    actor_loss::Float32 = 0.f0
    critic_loss::Float32 = 0.f0
    entropy_loss::Float32 = 0.f0
    loss::Float32 = 0.f0
end

(learner::A2CGAELearner)(obs::BatchObs) =
    learner.approximator.actor(send_to_device(
        device(learner.approximator),
        get_state(obs),
    )) |> send_to_host

function RLBase.update!(learner::A2CGAELearner, t::AbstractTrajectory)
    isfull(t) || return

    states = get_trace(t, :state)
    actions = get_trace(t, :action)
    rewards = get_trace(t, :reward)
    terminals = get_trace(t, :terminal)
    rollout = t[:state]

    AC = learner.approximator
    γ = learner.γ
    λ = learner.λ
    w₁ = learner.actor_loss_weight
    w₂ = learner.critic_loss_weight
    w₃ = learner.entropy_loss_weight

    states = send_to_device(device(AC), states)
    rollout = flatten_batch(rollout)
    rollout = send_to_device(device(AC), rollout)

    states_flattened = flatten_batch(states) # (state_size..., n_thread * update_step)
    actions = flatten_batch(actions)
    actions = CartesianIndex.(actions, 1:length(actions))

    rollout_values = AC.critic(rollout)
    rollout_values = send_to_host(rollout_values)
    rollout_values = reshape(
        rollout_values,
        size(states, ndims(states) - 1),
        size(states, ndims(states)) + 1,
    )
    advantages = generalized_advantage_estimation(
        rewards,
        rollout_values,
        γ,
        λ;
        dims = 2,
        terminal = terminals,
    )

    gains = advantages + select_last_dim(rollout_values, 1:(nframes(rollout_values)-1))
    gains = send_to_device(device(AC), gains)
    advantages = flatten_batch(advantages)
    advantages = send_to_device(device(AC), advantages)

    ps = Flux.params(AC)
    gs = gradient(ps) do
        logits = AC.actor(states_flattened)
        probs = softmax(logits)
        log_probs = logsoftmax(logits)
        log_probs_select = log_probs[actions]
        values = AC.critic(states_flattened)
        advantage = vec(gains) .- vec(values)
        actor_loss = -mean(log_probs_select .* advantages)
        critic_loss = mean(advantage .^ 2)
        entropy_loss = -sum(probs .* log_probs) * 1 // size(probs, 2)
        loss = w₁ * actor_loss + w₂ * critic_loss - w₃ * entropy_loss
        ignore() do
            learner.actor_loss = actor_loss
            learner.critic_loss = critic_loss
            learner.entropy_loss = entropy_loss
            learner.loss = loss
        end
        loss
    end

    if !isnothing(learner.max_grad_norm)
        learner.norm = clip_by_global_norm!(gs, ps, learner.max_grad_norm)
    end

    update!(AC, gs)
end

function (agent::Agent{<:QBasedPolicy{<:A2CGAELearner},<:CircularCompactSARTSATrajectory})(
    ::Training{PreActStage},
    obs,
)
    action = agent.policy(obs)
    state = get_state(obs)
    push!(agent.trajectory; state = state, action = action)
    update!(agent.policy, agent.trajectory)

    # the main difference is we'd like to flush the buffer after each update!
    if isfull(agent.trajectory)
        empty!(agent.trajectory)
        push!(agent.trajectory; state = state, action = action)
    end

    action
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using ExprSearch
using ExprSearch.CE

using SymbolicRegression
using RLESUtils, MemPools, RandChannels
using DepthAwarePCFGs
using LinearDerivTrees
using AbstractTrees

const MAXDEPTH = 10
const NUM_SAMPLES = 20
const RANDCHANNEL_WIDTH = 500

problem = Symbolic()
cfg = get_grammar(problem)
dpcfg = DepthAwarePCFG(cfg)
tree_params = LDTParams(cfg, typemax(Int64))
samples = [LinearDerivTree(tree_params; nodepool=MemPool(DerivTreeNode, 20, 500))
    for i=1:NUM_SAMPLES]
rc = RandChannel(NUM_SAMPLES, RANDCHANNEL_WIDTH)
wrc_vec = [WrappedRandChannel(rc, 0) for i=1:Threads.nthreads()]

ExprSearch.CE.ramped_rand!(wrc_vec, samples, dpcfg, MAXDEPTH)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
#=
Implement integer exponentiation. That is, implement the pow(x, y) function, where x and y are integers and returns x^y.

Do this faster than the naive method of repeated multiplication.

For example, pow(2, 10) should return 1024.
=#

function pow(x::Int, y::Int)
    if y == 0
        return 1
    elseif y % 2 == 0
        return pow(x, Int(y/2)) * pow(x, Int(y/2))
    else
        return x * pow(x,Int((y-1)/2)) * pow(x,Int((y-1)/2))
    end
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# host-level indexing

export allowscalar, @allowscalar, @disallowscalar, assertscalar

# mechanism to disallow scalar operations

@enum ScalarIndexing ScalarAllowed ScalarWarn ScalarWarned ScalarDisallowed

"""
    allowscalar(allow=true, warn=true)
    allowscalar(allow=true, warn=true) do end

Configure whether scalar indexing is allowed depending on the value of `allow`.

If allowed, `warn` can be set to throw a single warning instead. Calling this function will
reset the state of the warning, and throw a new warning on subsequent scalar iteration.

For temporary changes, use the do-block version, or [`@allowscalar`](@ref).
"""
function allowscalar(allow::Bool=true, warn::Bool=true)
    val = if allow && !warn
        ScalarAllowed
    elseif allow
        ScalarWarn
    else
        ScalarDisallowed
    end

    task_local_storage(:ScalarIndexing, val)
    return
end

@doc (@doc allowscalar) ->
function allowscalar(f::Base.Callable, allow::Bool=true, warn::Bool=false)
    val = if allow && !warn
        ScalarAllowed
    elseif allow
        ScalarWarn
    else
        ScalarDisallowed
    end

    task_local_storage(f, :ScalarIndexing, val)
end

"""
    assertscalar(op::String)

Assert that a certain operation `op` performs scalar indexing. If this is not allowed, an
error will be thrown ([`allowscalar`](@ref)).
"""
function assertscalar(op = "operation")
    val = get!(task_local_storage(), :ScalarIndexing) do
        if haskey(ENV, "JULIA_GPU_ALLOWSCALAR")
            parse(Bool, ENV["JULIA_GPU_ALLOWSCALAR"]) ? ScalarAllowed : ScalarDisallowed
        else
            ScalarWarn
        end
    end
    if val == ScalarDisallowed
        error("$op is disallowed")
    elseif val == ScalarWarn
        @warn "Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`"
        task_local_storage(:ScalarIndexing, ScalarWarned)
    end
    return
end

"""
    @allowscalar ex...
    @disallowscalar ex...
    allowscalar(::Function, ...)

Temporarily allow or disallow scalar iteration.

Note that this functionality is intended for functionality that is known and allowed to use
scalar iteration (or not), i.e., there is no option to throw a warning. Only use this on
fine-grained expressions.
"""
macro allowscalar(ex)
    quote
        task_local_storage(:ScalarIndexing, ScalarAllowed) do
            $(esc(ex))
        end
    end
end

@doc (@doc @allowscalar) ->
macro disallowscalar(ex)
    quote
        task_local_storage(:ScalarIndexing, ScalarDisallowed) do
            $(esc(ex))
        end
    end
end

# basic indexing with integers

Base.IndexStyle(::Type{<:AbstractGPUArray}) = Base.IndexLinear()

function Base.getindex(xs::AbstractGPUArray{T}, I::Integer...) where T
    assertscalar("scalar getindex")
    i = Base._to_linear_index(xs, I...)
    x = Array{T}(undef, 1)
    copyto!(x, 1, xs, i, 1)
    return x[1]
end

function Base.setindex!(xs::AbstractGPUArray{T}, v::T, I::Integer...) where T
    assertscalar("scalar setindex!")
    i = Base._to_linear_index(xs, I...)
    x = T[v]
    copyto!(xs, i, x, 1, 1)
    return xs
end

Base.setindex!(xs::AbstractGPUArray, v, I::Integer...) =
    setindex!(xs, convert(eltype(xs), v), I...)

# basic indexing with cartesian indices

Base.@propagate_inbounds Base.getindex(A::AbstractGPUArray, I::Union{Integer, CartesianIndex}...) =
    A[Base.to_indices(A, I)...]
Base.@propagate_inbounds Base.setindex!(A::AbstractGPUArray, v, I::Union{Integer, CartesianIndex}...) =
    (A[Base.to_indices(A, I)...] = v; A)

# generalized multidimensional indexing

Base.getindex(A::AbstractGPUArray, I...) = _getindex(A, to_indices(A, I)...)

function _getindex(src::AbstractGPUArray, Is...)
    shape = Base.index_shape(Is...)
    dest = similar(src, shape)
    any(isempty, Is) && return dest # indexing with empty array
    idims = map(length, Is)

    AT = typeof(src).name.wrapper
    # NOTE: we are pretty liberal here supporting non-GPU indices...
    gpu_call(getindex_kernel, dest, src, idims, adapt(AT, Is)...)
    return dest
end

@generated function getindex_kernel(ctx::AbstractKernelContext, dest, src, idims,
                                    Is::Vararg{<:Any,N}) where {N}
    quote
        i = @linearidx dest
        is = @inbounds CartesianIndices(idims)[i]
        @nexprs $N i -> I_i = @inbounds(Is[i][is[i]])
        val = @ncall $N getindex src i -> I_i
        @inbounds dest[i] = val
        return
    end
end

Base.setindex!(A::AbstractGPUArray, v, I...) = _setindex!(A, v, to_indices(A, I)...)

function _setindex!(dest::AbstractGPUArray, src, Is...)
    isempty(Is) && return dest
    idims = length.(Is)
    len = prod(idims)
    len==0 && return dest

    AT = typeof(dest).name.wrapper
    # NOTE: we are pretty liberal here supporting non-GPU sources and indices...
    gpu_call(setindex_kernel, dest, adapt(AT, src), idims, len, adapt(AT, Is)...;
             total_threads=len)
    return dest
end

@generated function setindex_kernel(ctx::AbstractKernelContext, dest, src, idims, len,
                                    Is::Vararg{<:Any,N}) where {N}
    quote
        i = linear_index(ctx)
        i > len && return
        is = @inbounds CartesianIndices(idims)[i]
        @nexprs $N i -> I_i = @inbounds(Is[i][is[i]])
        @ncall $N setindex! dest src[i] i -> I_i
        return
    end
end


text:
using Documenter

makedocs(
    sitename = "🎢 Control Theory",
    format = Documenter.LaTeX(),
    authors = "Joe Carpinelli",
    pages = [
        "Frontmatter" => [
            "Welcome!" => "index.md",
            "Chapter 1: Introduction" => "Chapter 1: Introduction.md"
        ],
        "Topic 1: Dynamics" => [
        "Chapter 2: General Dynamics" => "Chapter 2: General Dynamics.md",
        "Chapter 3: Flight Dynamics" => "Chapter 3: Flight Dynamics.md",
        "Chapter 4: Linear Dynamics" => "Chapter 4: Linear Dynamics.md",
        "Chapter 5: Equilibrium Points" => "Chapter 5: Equilibrium Points.md",
        "Chapter 6: Linearization" => "Chapter 6: Linearization.md"
        ],
        "Topic 2: Systems" => [],
        "Topic 3: Controls" => [],
        "Topic 4: Analysis" => []
    ]
)

max_stars_repo_name: cadojo/controls

max_issues_repo_name: cadojo/exploring-control-theory

max_forks_repo_name: cadojo/exploring-control-theory


text:
# Autogenerated wrapper script for FGlT_jll for i686-linux-musl
export libfglt

using CompilerSupportLibraries_jll
using cilkrts_jll
JLLWrappers.@generate_wrapper_header("FGlT")
JLLWrappers.@declare_library_product(libfglt, "libfglt.so.0")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, cilkrts_jll)
    JLLWrappers.@init_library_product(
        libfglt,
        "lib/libfglt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "EmbedOnehot" begin
    using Flux: onehot, onecold, data
    multi(x, n) = collect(Iterators.repeated(x, n))

    v = Vocabulary(["a", "b", "c", "d", "e"], "?")

    before_enc = ["a", "a", "c" ,"d", "e","x", "sad"]
    before_unk_enc = ["a", "a", "c" ,"d", "e","?", "?"]
    after_enc = [2,2,4,5,6,1,1]
    multi_after_enc = hcat(multi(after_enc, 5)...)

    e = Embed(10, v)

    @test data(e(after_enc)) == data(hcat(map(i->e.embedding[:,i], after_enc)...))
    @test data(e(before_enc)) == data(hcat(map(x->e.embedding[:, v(x)], before_unk_enc)...))
    @test onehot(e, before_enc) == onehotarray(length(v), v(before_unk_enc))
    @test onehot(e, after_enc) == onehotarray(length(v), after_enc)

    @test onehot(e, multi(before_enc,5)) == onehotarray(length(v), v(multi(before_unk_enc,5)))
    @test onehot(e, multi_after_enc) == onehotarray(length(v), multi_after_enc)

    @test onecold(e, onehot(e, multi(before_enc,5))) == multi(before_unk_enc,5)
    @test onecold(e, onehot(e, multi_after_enc)) == multi(before_unk_enc,5)
    @test onecold(e, onehot(e, before_enc)) == before_unk_enc
    @test onecold(e, onehot(e, after_enc)) == before_unk_enc
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FLINT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FLINT")
JLLWrappers.@generate_main_file("FLINT", UUID("e134572f-a0d5-539d-bddf-3cad8db41a82"))
end  # module FLINT_jll


text:
#=
The tree structure used for all the ClaDS functions
=#

"""
A phylogeny object for the module ClaDS. It is represented as a branch with a given length and optional attributes, and its daughter trees.

- `offsprings::Array{Tree,1}`: the two daughter trees.
- `branch_length::Float64`: the length of the branch
- `attributes::Array{T,1} where {T<:Number}`: soem attributes of the branch. Used to store the speciation rates.
- `n_nodes::Int64`: the number of nodes in the tree (internal nodes + tips)
- `extant::Bool`: does the tree have any extant species?
- `label::String`: if the tree is a tip (ie n_nodes == 1), contains the name of the corresponding species.
"""
struct Tree
    offsprings::Array{Tree,1}
    branch_length::Float64
    attributes::Array{T,1} where {T<:Number}
    n_nodes::Int64
    extant::Bool
    label::String
end

#=
Constructors (12 methods with various default options)
=#

# specify all
function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::T, n_nodes::Int64, extant::Bool, label::String) where {T<:Number}
    Tree(offsprings, branch_length, [attributes], n_nodes, extant, label)
end

# all except tree labels, set to empty strings
function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::T, n_nodes::Int64, extant::Bool) where {T<:Number}
    Tree(offsprings, branch_length, attributes, n_nodes, extant, "")
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::Array{T,1}, n_nodes::Int64, extant::Bool) where {T<:Number}
    Tree(offsprings, branch_length, attributes, n_nodes, extant, "")
end

# without extant information (set to true)
function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::Array{T,1}, n_nodes::Int64, label::String) where {T<:Number}
    Tree(offsprings, branch_length, attributes, n_nodes, true, label)
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::T, n_nodes::Int64, label::String) where {T<:Number}
    Tree(offsprings, branch_length, [attributes], n_nodes, true, label)
end

# without tip labels (set to empty strings) nor extant information (set to true)
function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::Array{T,1}, n_nodes::Int64) where {T<:Number}
    Tree(offsprings, branch_length, attributes, n_nodes, true)
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes::T, n_nodes::Int64) where {T<:Number}
    Tree(offsprings, branch_length, [attributes], n_nodes, true)
end

# whithout specification of the node number (set to the sum of the node number of the offspring +1)
function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes)
    n_nodes = length(branch_length)
    extant = false
    for i in offsprings
        n_nodes += i.n_nodes
        extant = extant || i.extant
    end
    extant = extant || n_nodes == 1
    Tree(offsprings, branch_length, attributes, n_nodes, extant)
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes, label::String)
    n_nodes = length(branch_length)
    extant = false
    for i in offsprings
        n_nodes += i.n_nodes
        extant = extant || i.extant
    end
    extant = extant || n_nodes == 1
    Tree(offsprings, branch_length, attributes, n_nodes, extant, label)
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes, extant::Bool)
    n_nodes = length(branch_length)
    for i in offsprings
        n_nodes += i.n_nodes
    end
    Tree(offsprings, branch_length, attributes, n_nodes, extant)
end

function Tree(offsprings::Array{Tree,1}, branch_length::Float64, attributes, extant::Bool, label::String)
    n_nodes = length(branch_length)
    for i in offsprings
        n_nodes += i.n_nodes
    end
    Tree(offsprings, branch_length, attributes, n_nodes, extant, label)
end

# empty tree
function Tree()
    Tree(Array{Tree,1}(undef,0), 0., Array{Float64,1}(undef,0), 0, false)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# Tests for Radial DEA Models
@testset "RadialDEAModel" begin

    ## Test Radial DEA Models with FLS Book data
    # Test against results with R
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Input oriented CRS
    deaio = dea(X, Y, orient = :Input, rts = :CRS)

    @test nobs(deaio) == 11
    @test ninputs(deaio) == 2
    @test noutputs(deaio) == 1
    @test efficiency(deaio) ≈ [1.0000000000;
                               0.6222896791;
                               0.8198562444;
                               1.0000000000;
                               0.3103709311;
                               0.5555555556;
                               1.0000000000;
                               0.7576690896;
                               0.8201058201;
                               0.4905660377;
                               1.0000000000]
    @test convert(Matrix, peers(deaio)) ≈
    [1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.4249783174  0  0 0.10928013877  0  0   0   0;
     1.134321653  0  0 0.4380053908  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.2573807721  0  0 0.04844814534  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.33333333333  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0 1.0348650979  0  0 0.11457435013  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.14814814815  0  0   0   0;
     0.000000000  0  0 0.4905660377  0  0 0.49056603774  0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   1.000000000]
     @test slacks(deaio, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                4.444444444   0;
                                0.000000000   0;
                                0.000000000   0;
                                1.640211640   0;
                                0.000000000   0;
                                0.000000000   4]
    @test slacks(deaio, :Y) ≈ zeros(11)

    # Otuput oriented CRS
    deaoo = dea(X, Y, orient = :Output, rts = :CRS)

    @test nobs(deaoo) == 11
    @test ninputs(deaoo) == 2
    @test noutputs(deaoo) == 1
    @test efficiency(deaoo) ≈ [1.0000000000;
                               1.606968641;
                               1.219726027;
                               1.0000000000;
                               3.221951220;
                               1.800000000;
                               1.0000000000;
                               1.319837398;
                               1.219354839;
                               2.038461538;
                               1.0000000000]
    @test convert(Matrix, peers(deaoo)) ≈
    [1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.6829268293  0  0 0.1756097561   0  0   0   0;
     1.383561644  0  0 0.5342465753  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0 0.8292682927  0  0 0.1560975610   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 0.6000000000   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0 1.3658536585  0  0 0.1512195122   0  0   0   0;
     0.000000000  0  0 0.0000000000  0  0 1.4000000000   0  0   0   0;
     0.000000000  0  0 1.0000000000  0  0 1.0000000000   0  0   0   0;
     1.000000000  0  0 0.0000000000  0  0 0.00000000000  0  0   0   0]
     @test slacks(deaoo, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                8   0;
                                0.000000000   0;
                                0.000000000   0;
                                2   0;
                                0.000000000   0;
                                0.000000000   4]
    @test slacks(deaoo, :Y) ≈ zeros(11)

    # Input oriented VRS
    deaiovrs = dea(X, Y, orient = :Input, rts = :VRS)

    @test nobs(deaiovrs) == 11
    @test ninputs(deaiovrs) == 2
    @test noutputs(deaiovrs) == 1
    @test efficiency(deaiovrs) ≈ [1.0000000000;
                                  0.8699861687;
                                  1.0000000000;
                                  1.0000000000;
                                  0.7116402116;
                                  1.0000000000;
                                  1.0000000000;
                                  1.0000000000;
                                  1.0000000000;
                                  0.4931209269;
                                  1.0000000000]
    @test convert(Matrix, peers(deaiovrs)) ≈
    [1.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.52558782849  0  0 0.0000000000  0  0.2842323651  0.1901798064   0  0   0   0;
     0.000000000    0  1 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.000000000    0  0 1.0000000000  0  0.00000000000 0.00000000000  0  0   0   0;
     0.56613756614  0  0 0.0000000000  0  0.4338624339  0.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  1.00000000000 0.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 1.00000000000  0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  1  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  1   0   0;
     0.03711078928  0  0 0.4433381608  0  0.00000000000 0.5195510500   0  0   0   0;
     0.000000000    0  0 0.0000000000  0  0.00000000000 0.00000000000  0  0   0   1.000000000]
     @test slacks(deaiovrs, :X) ≈ [0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0   0;
                                   0.000000000   0;
                                   0.000000000   0;
                                   0   0;
                                   0.000000000   0;
                                   0.000000000   4]
    @test slacks(deaiovrs, :Y) ≈ [0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  2.698412698;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000;
                                  0.000000000]

    # Output oriented VRS
    deaoovrs = dea(X, Y, orient = :Output, rts = :VRS)

    @test nobs(deaoovrs) == 11
    @test ninputs(deaoovrs) == 2
    @test noutputs(deaoovrs) == 1
    @test efficiency(deaoovrs) ≈ [1.0000000000;
                                  1.507518797;
                                  1.0000000000;
                                  1.0000000000;
                                  3.203947368;
                                  1.000000000;
                                  1.0000000000;
                                  1.000000000;
                                  1.000000000;
                                  1.192307692;
                                  1.0000000000]
    @test convert(Matrix, peers(deaoovrs)) ≈
    [1.000000000   0  0 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.38157894737 0  0 0.1710526316  0  0 0.4473684211   0  0   0   0;
     0.000000000   0  1 0.0000000000  0  0 0.00000000000  0  0   0   0;
     0.000000000   0  0 1.0000000000  0  0 0.00000000000  0  0   0   0;
     0.03947368421 0  0 0.7763157895  0  0 0.1842105263   0  0   0   0;
     0.000000000   0  0 0.0000000000  0  1 0.00000000000  0  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 1.00000000000  0  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  1  0   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  0  1   0   0;
     0.000000000   0  0 0.0000000000  0  0 0.00000000000  0  1   0   0;
     1.000000000   0  0 0.0000000000  0  0 0.00000000000  0  0   0   0]
    @test slacks(deaoovrs, :X) ≈ [0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                0.000000000   0;
                                5   11;
                                0.000000000   4]
    @test slacks(deaoovrs, :Y) ≈ zeros(11) atol=1e-15

    # Test no slacks
    deaionoslack = dea(X, Y, slack = false)
    @test efficiency(deaionoslack) == efficiency(deaio)
    @test isempty(slacks(deaionoslack, :X)) == 1
    @test isempty(slacks(deaionoslack, :Y)) == 1

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaio_ref_eff = zeros(size(X, 1))
    deaoo_ref_eff = zeros(size(X, 1))

    deaiovrs_ref_eff = zeros(size(X, 1))
    deaoovrs_ref_eff = zeros(size(X, 1))

    deaiovrs_ref_slackX = zeros(size(X))
    deaiovrs_ref_slackY = zeros(size(Y))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deaio_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Input, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaoo_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Output, rts = :CRS, Xref = Xref, Yref = Yref))[1]

        deaiovrs_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref))[1]
        deaoovrs_ref_eff[i] = efficiency(dea(Xeval, Yeval, orient = :Output, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaiovrs_ref_slackX[i,:] = slacks(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :X)
        deaiovrs_ref_slackY[i,:] = slacks(dea(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref), :Y)
    end

    @test deaio_ref_eff ≈ efficiency(deaio)
    @test deaoo_ref_eff ≈ efficiency(deaoo)

    @test deaiovrs_ref_eff ≈ efficiency(deaiovrs)
    @test deaoovrs_ref_eff ≈ efficiency(deaoovrs)

    @test deaiovrs_ref_slackX ≈ slacks(deaiovrs, :X) atol=1e-14
    @test deaiovrs_ref_slackY ≈ slacks(deaiovrs, :Y) atol=1e-15

    # Print
    show(IOBuffer(), deaio)
    show(IOBuffer(), deaionoslack)

    # Test errors
    @test_throws ErrorException dea([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws ErrorException dea([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws ErrorException dea([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws ErrorException dea([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ErrorException dea([1; 2; 3], [4; 5; 6], orient = :Error) # Invalid orientation
    @test_throws ErrorException dea([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale

end


text:
using ConvexTests, Tulip

@info "Starting Tulip tests"
do_tests("Tulip", () -> Tulip.Optimizer(); exclude = [r"mip", r"exp", r"socp", r"sdp"], description="Tests run with default parameters in type `Float64`.", variant="Float64", first = true, last = false)

do_tests("Tulip", () -> Tulip.Optimizer{BigFloat}(); exclude = [r"mip", r"exp", r"socp", r"sdp"], description="Tests run with default parameters in type `BigFloat`.", first = false, last = true, T=BigFloat, variant="BigFloat")


text:
import PyPlot

function plot1D(gp::GP; clim::Tuple{Float64, Float64}=(minimum(gp.X), maximum(gp.X)), CI::Float64=1.96, res::Int=1000, fname::AbstractString="")

    sx = (clim[2]-clim[1])/(res-1)
    x=collect(clim[1]:sx:clim[2])
    mu, Sigma = predict(gp, x)
    conf = CI*sqrt(Sigma)
    u = mu + conf
    l = mu - conf

    fig = PyPlot.figure("Gaussian Process Regression")
    ax = fig[:add_subplot](1,1,1)
    ax[:set_xlim](xmin=clim[1], xmax=clim[2])
    ax[:set_ylim](ymin=minimum(l), ymax=maximum(u))

    ax[:plot](x,mu,label="Mean function")       #plot the mean function
    ax[:scatter](gp.X,gp.y,marker="o",label="Observations",color="black")  #plot the observations
    ax[:plot](x,l,label="Confidence Region",color="red")            # plot the confidence region
    ax[:plot](x,u,color="red")
    ax[:fill_between](x,u,l, facecolor="blue", alpha=0.25)
    PyPlot.legend(loc = "best", fontsize=15);
    if fname != ""
        PyPlot.savefig(fname)
    end
end

function plot2D(gp::GP; clim::Tuple{Float64, Float64, Float64, Float64} = (minimum(gp.X[1,:]), maximum(gp.X[1,:]), minimum(gp.X[2,:]), maximum(gp.X[2,:])), res::Int=50, fname::AbstractString="")

    x = linspace(clim[1], clim[2], res)
    y = linspace(clim[3], clim[4], res)
    xgrid = repmat(x', res, 1 )
    ygrid = repmat(y, 1, res )

    mu = predict(gp,[vec(xgrid)';vec(ygrid)'])[1]
    zgrid  = reshape(mu,res,res)

    fig = PyPlot.figure("Contour Plot of Gaussian Process Regression")
    ax = fig[:add_subplot](1,1,1)
    ax[:set_xlim](xmin=clim[1], xmax=clim[2])
    ax[:set_ylim](ymin=clim[3], ymax=clim[4])
    ax[:contour](xgrid,ygrid,zgrid)
    ax[:scatter](gp.X[1,:],gp.X[2,:],marker="o",label="Observations",color="black")  #plot the observations
    if fname != ""
        PyPlot.savefig(fname)
    end
end

function PyPlot.plot(gp::GP; kwargs...)
    d, n = size(gp.X)
    PyPlot.close()
    if d>2
        error("Only 1D and 2D plots are permitted")
    elseif d==1
        plot1D(gp; kwargs...)
    elseif d==2
        plot2D(gp; kwargs...)
    end
end


text:
using LinearAlgebra
@testset "degopt_poly" begin
    x1a = [3; 4.0]
    x1b = [0.1; 0.3]
    x2a = [0.1; 0.1; -0.1]
    x2b = [0.3; 0; -0.1]
    x = [(x1a, x1b), (x2a, x2b)]

    z = [1.0; -0.1; 0.1; 0.01]
    (graph, crefs) = graph_degopt(x, z)

    x0 = 3.1415
    ev1 = eval_graph(graph, x0, comporder = get_topo_order_degopt(2))

    # Manual evaluation
    β1 = (x1a[1] + x1a[2] * x0) * (x1b[1] + x1b[2] * x0)
    β2 =
        (x2a[1] + x2a[2] * x0 + x2a[3] * β1) *
        (x2b[1] + x2b[2] * x0 + x2b[3] * β1)
    ev2 = z[1] + z[2] * x0 + z[3] * β1 + z[4] * β2
    @test ev1 ≈ ev2

    # Also check that the crefs are correct

    # Modify the graph: set almost all coeffs to one
    T = Float64
    t = ones(T, size(crefs, 1))
    t[2] = 2.0
    set_coeffs!(graph, t, crefs)

    # Create a new graph object graph2: with coeffs almost all zero
    x1a = [1; 2.0]
    x1b = [1; 1.0]
    x2a = [1.0; 1.0; 1.0]
    x2b = [1.0; 1.0; 1.0]
    z = ones(T, 4)
    x = [(x1a, x1b), (x2a, x2b)]
    (graph2, crefs2) = graph_degopt(x, z)

    x1 = 0.3
    @test eval_graph(graph, x1, comporder = get_topo_order_degopt(2)) ==
          eval_graph(graph2, x1)

    for k = 1:7
        x = Vector{Tuple{Vector{T},Vector{T}}}(undef, k)
        for i = 1:k
            x[i] = (collect(1.0:i+1), collect(-1.0:-1:-(i + 1)))
        end
        z = collect(range(4, 5, length = k + 2))
        (graph, cref) = graph_degopt(x, z)

        (xx, zz) = get_degopt_crefs(k)
        for i = 1:k
            for N = 1:2
                @test get_coeffs(graph, xx[i][N]) == x[i][N]
            end
        end
        @test get_coeffs(graph, zz) == z
    end

    # Check the inteface with all ones
    (graph, crefs) = graph_degopt(3)
    x = randn(3)
    a = (1 .+ x) .^ 2
    b = (1 .+ x .+ a) .^ 2
    c = (1 .+ x .+ a .+ b) .^ 2
    g = 1 .+ x .+ a .+ b .+ c
    @test eval_graph(graph, x) ≈ g
end


text:
using Documenter
include("../master/_core.jl")

makedocs()

run(`pandoc build/index.md --latex-engine=xelatex -o build/index.pdf`)

deploydocs(repo = "github.com/285714/ncm.git")