text:
@testset "multipole" begin
    #let's check that rotation works properly
    k0 = 0.1
    kin = 100k0
    λ0 = 2π/k0
    θ_i = 0.0

    N = 1000
    P = 5
    shapes = [squircle(0.2λ0, N)]
    centers = [-1.5λ0 0; 1.5λ0 0]
    Ns = size(centers,1)
    φs = [0.0; 0.23]
    ids = ones(Int, Ns)
    sp = ScatteringProblem(shapes, ids, centers, φs)
    points = λ0*[5 0; -2 3]

    verify_min_distance(shapes, centers, ids, points)
    β1, σ1 = solve_particle_scattering(k0, kin, P, sp, PlaneWave(θ_i);
        get_inner = true, verbose = false)
    Ez1 = scattered_field_multipole(k0, β1, centers, points)

    #now rotate everything by some angle and compare
    θ_r = 1.23#2π*rand()
    #notation is transposed due to structure of centers
    shapes2 = [squircle(0.2λ0, N); CircleParams(15)] #no effect,for extra code coverage
    centers2 = centers*[cos(θ_r) sin(θ_r); -sin(θ_r) cos(θ_r)]
    φs2 = θ_r .+ φs
    sp2 = ScatteringProblem(shapes2, ids, centers2, φs2)
    points2 = points*[cos(θ_r) sin(θ_r); -sin(θ_r) cos(θ_r)]

    β2, σ2 = solve_particle_scattering(k0, kin, P, sp2, PlaneWave(θ_i + θ_r);
        get_inner = true, verbose = true)
    Ez2 = scattered_field_multipole(k0, β2, centers2, points2)

    @test Ez1 ≈ Ez2
    β2_r = copy(β2)
    for ic = 1:Ns, p = -P:P
        β2_r[(ic-1)*(2P+1) + p + P + 1] *= exp(1.0im*θ_r*p)
    end
    @test β1 ≈ β2_r
    @test σ1 ≈ σ2
end

@testset "circle" begin
    k0 = 1
    kin = 0.1
    R = 0.2
    P = 10
    S, gamma = ParticleScattering.circleScatteringMatrix(k0, kin, R, P; gamma = true)
    center = [0 2R]
    points = center .+ [0 0.5R; 0.1R -0.2R]
    E1 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points)
    E2 = ParticleScattering.innerFieldCircle(kin, gamma, center[1,:], points[1,:])
    @test E1[1] == E2
end


text:
include("GeneralMethods.jl")
include("SuperOperator.jl")
include("Ising.jl")

#using SymPy

###########################
###Core of the algorithm###
###########################
"""
    joinDensityMatrix!(extendedRho, matRho, TracesSmall, invMatJoin, pseudoIdentity)

Build rho_0:l+1 from rho_0:l-1, rho_1:l and rho_2:l+1
extendedRho stocks the larger DMs
matRho are the current DMs
TracesSmall are the trace matrix for matRho
invMatJoin is the pseudo-inverse to find the marginal
pseudoIdentity is used to minimize the distance with the previous matrix
"""
function joinDensityMatrix!(extendedRho, matRho, TracesSmall, invMatJoin, pseudoIdentity)
    for l in 2:length(matRho)-1
        joinDensityMatAux!(extendedRho[l],  matRho[l-1], matRho[l], matRho[l+1], invMatJoin, pseudoIdentity)
    end
    joinDensityMatAux!(extendedRho[1],  leftFusionIdentity(reduceDensityMatrixRight(matRho[1], TracesSmall)), matRho[1], matRho[2], invMatJoin, pseudoIdentity)
    joinDensityMatAux!(extendedRho[end],  matRho[end-1], matRho[end], rightFusionIdentity(reduceDensityMatrixLeft(matRho[end], TracesSmall)), invMatJoin, pseudoIdentity)
end

"""
    joinDensityMatAux!(rhoGuess, rhoLeft, rhoMiddle,  rhoRight, invMatJoin, pseudoIdentity)

Build rhoGuess from rhoLeft, rhoMiddle and rhoRight
rhoGuess is both the initial guess and the returned matrix
rhoLeft is the initial left matrix
rhoMiddle is the initial middle matrix, that we are approximating
rhoRight is the initial right matrix
invMatJoin is the pseudo-inverse to find the marginal
pseudoIdentity is used to minimize the distance with the previous matrix
"""
function joinDensityMatAux!(rhoGuess, rhoLeft, rhoMiddle,  rhoRight, invMatJoin, pseudoIdentity)
    len=getLength(length(rhoGuess));
    ###Fast method
    temp=transpose(At_mul_B(vcat(rhoLeft, rhoMiddle, rhoRight), invMatJoin))
    temp-=transpose(At_mul_B(rhoGuess, pseudoIdentity))
    temp/=2
    tempconj=reshape(transpose(conj.(reshape(temp, (2^(len), 2^(len))))), length(temp))  ###A ameliorer
    for j=1:length(rhoGuess)
        rhoGuess[j]+=temp[j]+tempconj[j]
    end
    return 0
end

"""
    evolveMatrix!(extendedRho, ULeft, tULeft)

Let the extended matrix evolve under Uleft and its transconjugate. ToDO: improve
"""
function evolveMatrix!(extendedRho,temp,  ULeft, tULeft)
    lenmat=size(ULeft)[1];
    for l in 1:length(extendedRho)
        A_mul_B!(temp,ULeft, reshape(view(extendedRho[l], :), (lenmat, lenmat)))
        A_mul_B!(reshape(view(extendedRho[l], :), (lenmat, lenmat)), temp, tULeft)
    end
end

"""
    massReduce!(matRho, timeEvolvedRho, Traces)

Update matRho from the evolved larger DMs
matRho are the current DMs to be updated
timeEvolvedRho are the time evolved larger DMs
Traces are the trace matrix for extendedRho
"""
function massReduce!(matRho, timeEvolvedRho, Traces)
    len=getLength(length(matRho[1]));
    for l=1:length(matRho)
        matRho[l]=buildReducedDensityMatrix(timeEvolvedRho[l], 1, len, Traces)
        normalizeDM!(matRho[l])
    end
end

"""
    extendZone!(matRho,  indexmin, extendedRho, Traces, TracesSmall)

Check if one needs to extend the considered range
matRho are the current DMs
indexmin is the lefttmost site in matRho
extendedRho are the current extendedDM (odd or even)
Traces are the trace matrix for extendedRho
TracesSmall are the trace matrix for matRho
"""
function extendZone!(matRho,  indexmin, extendedRho, Traces, TracesSmall)
    len=getLength(length(matRho[1]));
    ##Check if left is entangled
    mini=indexmin;
    mi=mutualInformation(matRho[1], 0, 1, len, TracesSmall)
    if mi>1e-7
        unshift!(matRho, buildReducedDensityMatrix(extendedRho[1], 0, len, Traces))
        unshift!(extendedRho, leftFusionIdentity(reduceDensityMatrixRight(extendedRho[1], Traces)))
        mini-=1
    end
    ##Same for right extremity
    mi=mutualInformation(matRho[end], 0, len-1, len, TracesSmall)
    if mi>1e-7
        push!(matRho, buildReducedDensityMatrix(extendedRho[end], 2, len, Traces))
        push!(extendedRho, rightFusionIdentity(reduceDensityMatrixLeft(extendedRho[end], Traces)))
    end
    return mini
end

###################
###Main function###
###################
function main(hx::Float64, hz::Float64, T::Float64; maxLen=3, dt=0.0001)
    #Define Pauli Matrices
    sx=[0 1; 1 0]; sz=-[1  0 ; 0 -1]; sy=[0 -1im; 1im 0]
    tim=dt:dt:T
    #Define Hamiltonian
    localHamiltonianSingleSite=buildSingleSiteHamiltonian(hx, hz, maxLen+2)  #### Build local Hamiltonian i:i+1 for i=0:maxLen-1 for maxLen+1 sites matrices.
    localHamiltonianMultiSite=buildMultiSiteHamiltonian(hx, hz, maxLen+2)  #### Build local Hamiltonian i:i+1 for i=0:maxLen-1 for maxLen+1 sites matrices.
    ###Evolution operator
    ULeft=expm(full(-1im*dt/2*localHamiltonianSingleSite))*expm(full(-1im*dt*localHamiltonianMultiSite))*expm(full(-1im*dt/2*localHamiltonianSingleSite))
    tULeft=ULeft'
    #StockedOperator
    extendedRho=Array{Array{Complex128, 1}}((maxLen));
    for l in 1:length(extendedRho)
        extendedRho[l]=reshape(eye(2^(maxLen+2), 2^(maxLen+2))/2^(maxLen+2),  2^(maxLen*2+4))
    end
    tempMatrix=reshape(similar(extendedRho[1]), (size(ULeft)[1],size(ULeft)[1]))
    ###Setting up the global trace operators
    Traces=initializeTraceMatrices(maxLen+2);
    TracesSmall=initializeTraceMatrices(maxLen);
    MatJoin=transpose(sparse(vcat(Traces[1, 3]', Traces[2,2]', Traces[3,1]')));
    invMatJoin=sparse(triming.(pinv(full(MatJoin))));
    pseudoIdentity=sparse(triming.(MatJoin*invMatJoin));
    ###Observables###
    if length(tim)<=1000
        tim2=0:10*dt:T;
    else
        tim2=0:100*dt:T;
    end
    szObservable=zeros((length(tim)+99)÷100+1);
    sxObservable=zeros((length(tim)+99)÷100+1);
    #initializeState
    matRho=map(x->reshape(x, length(x)), generateInitialState_z(maxLen));
    minMatRho=-maxLen+1;
    #timeEvolution
    tic()
    cnt=0; totaltime=0; cnt2=2;
    szObservable[1]=real(trace(sz*reshape(buildReducedDensityMatrix(matRho[1-minMatRho], 0, TracesSmall), (2, 2))))
    sxObservable[1]=real(trace(sx*reshape(buildReducedDensityMatrix(matRho[1-minMatRho], 0, TracesSmall), (2, 2))))
    for j in 1:length(tim)
        ####Convention: apply Odd then Even
        if mod(j, (length(tim)+99)÷100)==0
            cnt+=1
            ctime=toq()
            println(string(cnt)*"% of simulation done. Time for the last %: "*string(round(ctime, 2))*"sec.")
            totaltime+=ctime
            println("Total elapsed time: "*string(round(totaltime, 2))*"sec.")
            println((minMatRho, minMatRho+length(matRho)+maxLen-2))
            tic()
        end
        if tim[j]==tim2[cnt2]
            szObservable[cnt2]=real(trace(sz*reshape(buildReducedDensityMatrix(matRho[1-minMatRho], 0, TracesSmall), (2, 2))))
            sxObservable[cnt2]=real(trace(sx*reshape(buildReducedDensityMatrix(matRho[1-minMatRho], 0, TracesSmall), (2, 2))))
            cnt2+=1
        end
        joinDensityMatrix!(extendedRho, matRho, TracesSmall, invMatJoin, pseudoIdentity);   #Build the larger density Matrix
        evolveMatrix!(extendedRho, tempMatrix, ULeft, tULeft);   #Evolve the larger densityMatrix
        massReduce!(matRho, extendedRho, Traces) #Compute the update densityMatrix
        minMatRho=extendZone!(matRho, minMatRho, extendedRho, Traces, TracesSmall)
        ##Measure sz
    end
    toc()
    return tim2, szObservable, sxObservable
end

###ToDo: remplacer isposdef par issemiposdef
###ToDo: optimize matrix product in joinDensity ? Have an explicit function compiled probably.

ext: jl

max_stars_repo_name: LoicHerviou/DensitySimulation

max_issues_repo_name: LoicHerviou/DensitySimulation

max_forks_repo_name: LoicHerviou/DensitySimulation


text:
function Source!(F,U,CG,Param)
nz=Param.Grid.nz;
if strcmp(Param.Thermo,"Energy")
else
  @views Pres=Pressure(U[:,:,Param.ThPos],U[:,:,Param.RhoPos],U[:,:,Param.RhoPos],Param);
end
sigma=Pres/Param.p0;
@views T=Pres./(Param.Rd*U[:,:,Param.RhoPos]);

matlab_max2(X) = max.(X, 0)
height_factor = matlab_max2((sigma .- Param.sigma_b) ./ (1-Param.sigma_b));

temp = (Param.T_equator .-
    Param.DeltaT_y*repmat(sin.(Param.latN).^2,1,nz) .-
    Param.DeltaTh_z*log.(sigma).*repmat(cos.(Param.latN).^2,1,nz)) .*
    sigma.^(Param.Rd/Param.Cpd)

matlab_max(X) = max.(X, Param.T_min)

@views DeltaRhoT=(Param.k_a .+ (Param.k_s - Param.k_a) .* height_factor .*
  repmat(cos.(Param.latN).^4,1,nz)) .* U[:,:,Param.RhoPos] .*
  (T .- matlab_max(temp));
@views F[:,:,Param.uPos:Param.vPos] .-= Param.k_f*height_factor.*U[:,:,Param.uPos:Param.vPos];
if strcmp(Param.Thermo,"Energy")
else
@views  F[:,:,Param.ThPos] .-= DeltaRhoT.*(Param.p0./Pres).^Param.kappa;
end
end


text:
include("pm_io/matpower.jl")
include("pm_io/common.jl")
include("pm_io/pti.jl")
include("pm_io/psse.jl")
include("pm_io/data.jl")
include("pm_io/multiconductor.jl")


text:
using Test

input = joinpath(@__DIR__, "input")
tiles = map(t -> split(t, "\n"), split(read(input, String), "\n\n"))

mutable struct Tile
    id::Int
    raw::Array{Char, 2}
end

function copy(tile::Tile)
    new = copy(tile.raw)
    return Tile(tile.id, new)
end

t = Tile(1, ['a' 'b' 'c'; 'd' 'e' 'f'])

function size(tile::Tile)
    return Base.size(tile.raw)
end

@test size(t) == (2, 3)

function top(tile::Tile)
    (y, x) = size(tile)
    return tile.raw[range(1, length=x, step=y)]
end

@test top(t) == ['a', 'b', 'c']

function right(tile::Tile)
    (y, x) = size(tile)
    return tile.raw[range((x * y) - y + 1, length=y, step=1)]
end

@test right(t) == ['c', 'f']

function bottom(tile::Tile)
    (y, x) = size(tile)
    return tile.raw[range(y, length=x, step=y)]
end

@test bottom(t) == ['d', 'e', 'f']

function left(tile::Tile)
    (y, x) = size(tile)
    return tile.raw[range(1, length=y, step=1)]
end

@test left(t) == ['a', 'd']

function flip!(tile)
    new = tile.raw[:,end:-1:1]
    tile.raw = new
end

function rotate!(tile)
    new = rotr90(tile.raw)
    tile.raw = new
end

function parse_tiles(tiles)
    return map(tile -> parse_tile(tile), tiles)
end

function parse_tile(tile)
    m = match(r"(\d+)", first(tile))
    id = parse(Int, m.captures[1])
    raw = reduce(vcat, permutedims.(collect.(tile[2:end])))
    return Tile(id, raw)
end

function find_edges(tile)
    edges = []
    push!(edges, top(tile))
    push!(edges, right(tile))
    push!(edges, bottom(tile))
    push!(edges, left(tile))
    push!(edges, reverse(top(tile)))
    push!(edges, reverse(right(tile)))
    push!(edges, reverse(bottom(tile)))
    push!(edges, reverse(left(tile)))
    return edges
end

function orient!(tile, matches)
    for flipped in (false, true), i in 1:4
        matches(tile) && return
        flipped && flip!(tile)
        flipped && matches(tile) && return
        flipped && flip!(tile)
        rotate!(tile)
    end
end

function find_corners(tiles)
    corners = []
    ec = count_edges(tiles)

    for tile in tiles
        c = sum([length(ec[top(tile)]), length(ec[bottom(tile)]), length(ec[left(tile)]), length(ec[right(tile)])])
        c == 6 && push!(corners, tile)
    end
    return corners
end

function count_edges(tiles)
    acc = Dict{Vector{Char}, Vector{Int}}()
    for tile in tiles
        edges = find_edges(tile)
        for edge in edges
            haskey(acc, edge) ? acc[edge] = append!(acc[edge], tile.id) : acc[edge] = [tile.id]
        end
    end
    return acc
end

function find_origin(tiles)
    edges = count_edges(tiles)
    corner = first(find_corners(tiles))
    for _ in 1:4
        t = edges[top(corner)]
        l = edges[left(corner)]
        length(t) == 1 && length(l) == 1 && return corner
        rotate!(corner)
    end
end

function trim(tile)
    raw = tile.raw
    (x, y) = size(tile)
    return raw[2:(x - 1),2:(y - 1)]
end

function assemble(tiles)
    origin = find_origin(deepcopy(tiles))
    edges = count_edges(deepcopy(tiles))

    rows = []
    row = [origin]

    rowstart = origin
    current = origin

    while !isnothing(current)
        match = edges[right(current)]
        if length(match) == 2
            next_id = only(filter(id -> id != current.id, match))
            next_idx = findfirst(t -> t.id == next_id, tiles)
            next = tiles[next_idx]
            matcher = (t) -> left(t) == right(current)
            orient!(next, matcher)
            push!(row, next)
            current = next
        else
            push!(rows, row)
            match = edges[bottom(rowstart)]
            if length(match) == 2
                next_id = only(filter(id -> id != rowstart.id, match))
                next_idx = findfirst(t -> t.id == next_id, tiles)
                next = tiles[next_idx]
                orient!(next, (t) -> top(t) == bottom(rowstart))
                current = next
                rowstart = current
                row = [rowstart]
            else
                current = nothing
            end
        end
    end

    return vcat(map(row -> hcat(map(tile -> trim(tile), row)...), rows)...)
end

function find_seamonsters(assembled)
    # seamonster pattern
    indices = [
        CartesianIndex(19, 1),
        CartesianIndex(1, 2),
        CartesianIndex(6, 2),
        CartesianIndex(7, 2),
        CartesianIndex(12, 2),
        CartesianIndex(13, 2),
        CartesianIndex(18, 2),
        CartesianIndex(19, 2),
        CartesianIndex(20, 2),
        CartesianIndex(2, 3),
        CartesianIndex(5, 3),
        CartesianIndex(8, 3),
        CartesianIndex(11, 3),
        CartesianIndex(14, 3),
        CartesianIndex(17, 3)
    ]

    (xr, yr) = Base.size(assembled) .- (20, 3)

    counts = []

    for flipped in (false, true), i in 1:4
        flipped ? assembled = assembled[end:-1:1,:] : nothing
        count = 0
        idx = CartesianIndex(0,0)
        for j in 1:yr
            (x, y) = Tuple(idx)
            idx = CartesianIndex(0, y + 1)
            for i in 1:xr
                idx += CartesianIndex(1, 0)
                all(i -> assembled[idx + i] == '#', indices) ? count += 1 : nothing
            end
        end
        push!(counts, count)
        flipped ? assembled = assembled[end:-1:1,:] : nothing
        assembled = rotr90(assembled)
    end
    return counts
end

function check(puzzle)
    acc = 0
    for i in eachindex(puzzle)
        puzzle[i] == '#' ? acc += 1 : nothing
    end
    return acc
end

tiles = parse_tiles(tiles)
corners = find_corners(tiles)

assembled = assemble(tiles)
counts = find_seamonsters(assembled)

p1 = prod(map(tile -> tile.id, corners))
p2 = check(assembled) - (maximum(counts) * 15)

println("-----------------------------------------------------------------------")
println("jurassic jigsaw -- part one :: $p1")
println("jurassic jigsaw -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 174206308298779)
@assert(p2 == 2409)


text:
function phi(x, ::Real, ::Product)
	return asin(x)
end

function phi(x, b::Real, ::Union{ChebyshevSparse, ChebyshevTower})
	return b*acos(x)
end


text:
"""
    JuliaEcosystem.jl

    This module is used to generate a text file with
        name, version, dependency, repository
    information for all registered installable packages
"""
module JuliaEcosystem
    using Base.Iterators: flatten
    using CSV: CSV, File
    using HTTP: request
    using Pkg: METADATA_compatible_uuid
    using Pkg.Operations: load_package_data_raw, deps_graph, simplify_graph!, resolve
    using Pkg.TOML: parsefile
    using Pkg.Types: Context, Fixed, Requires, UUID, uuid_julia, VersionRange, VersionSpec

    export PkgMetaData

    # Record information about the registry
    """
        All package names in the General registry
    """
    const pkgs =
        readdir.(joinpath.(homedir(), ".julia/registries/General", string.('A':'Z'))) |>
        flatten |>
        collect |>
        (x -> filter!(x -> ~any(x ∈ ["julia", ".DS_Store"]), x))
    const deps = Dict{UUID,Dict{VersionRange,Dict{String,UUID}}}()
    const compat = Dict{UUID,Dict{VersionRange,Dict{String,VersionSpec}}}()
    const uuid_to_name =
        Dict(UUID("ade2ca70-3891-5945-98fb-dc099432e06a")=>"Dates",
             UUID("7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee")=>"FileWatching",
             UUID("ea8e919c-243c-51af-8825-aaa63cd721ce")=>"SHA",
             UUID("6462fe0b-24de-5631-8697-dd941f90decc")=>"Sockets",
             UUID("8bb1440f-4735-579b-a4ab-409b98df4dab")=>"DelimitedFiles",
             UUID("9a3f8284-a2c9-5f02-9a11-845980a1fd5c")=>"Random",
             UUID("8ba89e20-285c-5b6f-9357-94700520ee1b")=>"Distributed",
             UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e")=>"LinearAlgebra",
             UUID("a63ad114-7e13-5084-954f-fe012c677804")=>"Mmap",
             UUID("2f01184e-e22b-5df5-ae63-d93ebab69eaf")=>"SparseArrays",
             UUID("cf7118a7-6976-5b1a-9a39-7adc72f591a4")=>"UUIDs",
             UUID("56ddb016-857b-54e1-b83d-db4d58db5568")=>"Logging",
             UUID("10745b16-79ce-11e8-11f9-7d13ad32a3b2")=>"Statistics",
             UUID("8bf52ea8-c179-5cab-976a-9e18b702a9bc")=>"CRC32c",
             UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f")=>"Pkg",
             UUID("9e88b42a-f829-5b0c-bbe9-9e923198166b")=>"Serialization",
             UUID("d6f4376e-aef5-505a-96c1-9c027394607a")=>"Markdown",
             UUID("4607b0f0-06f3-5cda-b6b1-a6196a1729e9")=>"SuiteSparse",
             UUID("3fa0cd96-eef1-5676-8a61-b3b8758bbffb")=>"REPL",
             UUID("8f399da3-3557-5675-b5ff-fb832c97cbdb")=>"Libdl",
             UUID("2a0f44e3-6c83-55bd-87e4-b1978d98bd5f")=>"Base64",
             UUID("9abbd945-dff8-562f-b5e8-e1ebf5ef1b79")=>"Profile",
             UUID("b77e0a4c-d291-57a0-90e8-8db25a27a240")=>"InteractiveUtils",
             UUID("1a1011a3-84de-559e-8e89-a11a2f7dc383")=>"SharedArrays",
             UUID("de0858da-6303-5e67-8744-51eddeeeb8d7")=>"Printf",
             UUID("9fa8497b-333b-5362-9e8d-4d0656e87820")=>"Future",
             UUID("76f85450-5226-5b5a-8eaa-529ad045b433")=>"LibGit2",
             UUID("8dfed614-e22c-5e08-85e1-65c5234f0b40")=>"Test",
             UUID("4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5")=>"Unicode",
             uuid_julia=>"julia")
    const versions = Dict{UUID,Set{VersionNumber}}()
    for pkg ∈ pkgs
        dir = joinpath(homedir(), ".julia/registries/General", uppercase(pkg[1:1]), pkg)
        uuid = UUID(parsefile(joinpath(dir, "Package.toml"))["uuid"])
        uuid_to_name[uuid] = pkg
        versions[uuid] = Set(VersionNumber.(keys(parsefile(joinpath(dir, "Versions.toml")))))
        deps[uuid] = load_package_data_raw(UUID, joinpath(dir, "Deps.toml"))
        compat[uuid] = load_package_data_raw(VersionSpec, joinpath(dir, "Compat.toml"))
    end
    """
        PkgMetaData(pkg::AbstractString)

        Given a pkg, it returns the name, version, dependencies, and repository if installable
    """
    struct PkgMetaData
        name::String
        uuid::UUID
        repo::String
        version::VersionNumber
        dependencies::Vector{String}
        installable::Bool
        function PkgMetaData(pkg::AbstractString)
            # First get the package information
            dir = joinpath(homedir(), ".julia/registries/General", pkg[1:1], pkg)
            package = parsefile(joinpath(dir, "Package.toml"))
            name = package["name"]
            repo = replace(package["repo"], r"\.git$" => "")
            uuid = UUID(package["uuid"])
            # Verify that it can be installed for Julia 1 (current LTS release)
            try
                graph = deps_graph(Context(),
                                   uuid_to_name,
                                   Requires(uuid => VersionSpec()),
                                   Dict(uuid_julia => Fixed(v"1")))
                simplify_graph!(graph)
                solved = resolve(graph)
                version = solved[uuid]
                dependencies = get.(Ref(uuid_to_name),
                                    filter(!isequal(uuid), keys(solved)),
                                    nothing) |>
                               sort!
                # Check whether the repository exists or redirects
                # If it isn't accessible it is not installable
                # If it redirects, updates to the current resource location
                response = request("GET", repo)
                if !isa(response.request.parent, Nothing)
                    # It redirected so replace with final address
                    repo = header(response.request.parent, "Location")
                end
                new(name, uuid, repo, version, dependencies, true)
            catch
                new(name, uuid, repo, VersionNumber(0), Vector{String}(), false)
            end
        end
    end

end


text:
module Compute

export demand!, demandD1only!, demandOnly
export setupIndex!, jacobS!, passThrough
using TypeDef, Helper, Tools
using Statistics, LinearAlgebra, LoopVectorization

function share!(
  s,
  p,
  num,
  μ̃,
  αᵢ::Vector{Float64},
  active::BitArray{1},
  dim::dimType)

  for i=1:dim.person
    den = 1.0;
    for j=1:dim.product
      num[j] = (active[j]==true) ? exp(μ̃[j,i]+αᵢ[i]*p[j]) : 0
      den += num[j]
    end
    for j=1:dim.product
      s[j] = num[j]/den
    end
  end
end

function demand!(
  sim::simType,
  p::Vector{Float64},
  dim::dimType,
  mat::matType,
  par::parType,
  indx::indxType)

  s,num,D,ms,∂s∂p,∂S∂p,∂²S∂p²,∂D∂p,∂²D∂p² = alias(dim,sim,true)

  @inbounds for r=1:dim.region
    for k=1:dim.product
      for j=1:dim.product
        ∂S∂p[j,k] = 0
        for i=1:dim.product
          ∂²S∂p²[i,j,k] = 0
        end
      end
    end
    αᵢ,μ̃,pop,active = aliasRegion!(num,indx,r,dim,mat,par,sim)

    for n=1:dim.person
      jacobS!(∂s∂p,∂S∂p,s,view(ms,:,r),view(μ̃,:,n),αᵢ[n],p,num,active,dim)

      # ∂²S∂p²[i,j,k] = ∂sᵢ / (∂p_j ∂p_k)
      @avx for ki in eachindex(active)
        k = active[ki]
        for ji in eachindex(active)
          j = active[ji]
          for ii in eachindex(active)
            i = active[ii]
            sp = αᵢ[n]*(∂s∂p[i,k]*(i==j)-s[j]*∂s∂p[i,k]-s[i]*∂s∂p[j,k])
            ∂²S∂p²[i,j,k] += sp
          end
        end
      end
    end # end of n loop

    for k in active
      ms[k,r] /= dim.person
      D[k]    += ms[k,r] * pop
      for j in active
        ∂S∂p[j,k] /= dim.person
        ∂D∂p[j,k] += pop * ∂S∂p[j,k]
        for i in active
          ∂²S∂p²[i,j,k] /= dim.person
          ∂²D∂p²[i,j,k] += pop * ∂²S∂p²[i,j,k]
        end
      end
    end
  end # end of r loop
end

function demandD1only!(
  sim::simType,
  p::Vector{Float64},
  dim::dimType,
  mat::matType,
  par::parType,
  indx::indxType)

  s,num,D,ms,∂s∂p,∂S∂p,∂D∂p = alias(dim,sim)

  @inbounds for r=1:dim.region
    fill!(∂S∂p,0)
    # fastFill!(∂S∂p,dim.product,dim.product)
    αᵢ,μ̃,pop,active = aliasRegion!(num,indx,r,dim,mat,par,sim)

    for n=1:dim.person
      jacobS!(∂s∂p,∂S∂p,s,view(ms,:,r),view(μ̃,:,n),αᵢ[n],p,num,active,dim)
    end # end of n loop

    for k in active
      ms[k,r] /= dim.person
      D[k]  += ms[k,r] * pop
      for j in active
        ∂S∂p[j,k] /= dim.person
        ∂D∂p[j,k] += pop * ∂S∂p[j,k]
      end
    end
  end # end of g loop
end

function demandOnly(
  sim::simType,
  p::Vector{Float64},
  dim::dimType,
  mat::matType,
  par::parType,
  indx::indxType)

  s,num,D,ms,∂s∂p,∂S∂p,∂D∂p = alias(dim,sim)

  csplusTotal = 0.0
  @inbounds for r=1:dim.region
    αᵢ,μ̃,pop,active = aliasRegion!(num,indx,r,dim,mat,par,sim)

    csplus = 0.0
    for n=1:dim.person
      den = 1.0
      @avx for ji in eachindex(active)
        j = active[ji]
        num[j] = exp(μ̃[j,n]+αᵢ[n]*p[j])
        den += num[j]
      end
      csplus += -3log(den) * pop / αᵢ[n]
      @avx for ji in eachindex(active)
        j = active[ji]
        s[j]     = num[j]/den
        ms[j,r] += s[j]
      end
    end # end of n loop
    csplus /= dim.person
    csplusTotal += csplus

    for k in active
      ms[k,r] /= dim.person
      D[k]    += ms[k,r] * pop
    end
  end # end of g loop

  return copy(D), csplusTotal
end

function jacobS!(
  ∂s∂p,
  ∂S∂p,
  s::Vector{Float64},
  ms,
  μ̃,
  αᵢ::Float64,
  p::Vector{Float64},
  num::Vector{Float64},
  active,
  dim::dimType)

  den = 1.0
  @avx for ji in eachindex(active)
    j = active[ji]
    num[j] = exp(μ̃[j]+αᵢ*p[j])
    den += num[j]
  end

  @avx for ji in eachindex(active)
    j = active[ji]
    s[j]   = num[j]/den
    ms[j] += s[j]
  end
  @avx for ji in eachindex(active)
    j = active[ji]
    for ii in eachindex(active)
      i = active[ii]
      ∂s∂p[i,j]  = αᵢ*(s[i]*(i==j) - s[i]*s[j])
      ∂S∂p[i,j] += ∂s∂p[i,j]
    end
  end

end

function passThrough(
  markupR,
  markupW,
  markupV,
  dim::dimType,
  mat::matType,
  sim::simType,
  indx::indxType,
  mask::maskType)

  A = zeros(dim.product,dim.product)
  B = zeros(dim.product,dim.mvno)
  Djpi = sim.∂D∂p' .* mask.Djpi
  any = indx.anyActive

  for f=1:4
    Lf = (mat.firm[:,2].==f) .& any
    Lfid = findall(Lf)
    a1 = sim.∂D∂p[Lf,:]
    for l in Lfid
      a1 = a1 + markupR[l] * sim.∂²D∂p²[Lf,:,l]
    end

    if f<=3
      L0 = mat.mvno[:,f] .& any
      L0id = findall(L0)

      a2 = sim.∂D∂p[L0,:]
      for l in L0id
        a1 = a1 + markupW[l] * sim.∂²D∂p²[Lf,:,l]
        a2 = a2 + markupV[l] * sim.∂²D∂p²[L0,:,l]
      end
      A[Lf.|L0,:] = A[Lf.|L0,:] + [a1;a2] + Djpi[Lf.|L0,:]
      B[Lf.|L0,indx.mvnp[:,f]] = [-sim.∂D∂p[L0,Lf]'; sim.∂D∂p[L0,L0]]
    else
      A[Lf,:] = A[Lf,:] + a1 + Djpi[Lf,:]
    end
  end
  x = zeros(dim.product,dim.mvno)
  # x[any,:] = inv(A[any,any]) * B[any,:]
  x[any,:] = A[any,any] \ B[any,:]

  return x
end

end # end of module

max_stars_repo_name: yutec/FightBrand

max_issues_repo_name: yutec/FightBrand

max_forks_repo_name: yutec/FightBrand


text:
using Arblib, Test, LinearAlgebra

@testset "Arblib" begin
    include("arb_types-test.jl")
    include("types-test.jl")
    include("arbcall-test.jl")
    include("precision-test.jl")
    include("setters-test.jl")
    include("constructors-test.jl")
    include("predicates-test.jl")
    include("show-test.jl")
    include("examples.jl")
    include("arithmetic.jl")
    include("rand.jl")
    include("vector.jl")
    include("matrix.jl")
    include("eigen.jl")
    include("calc_integrate.jl")
    include("ref-test.jl")
    include("poly.jl")
    include("series.jl")
end


text:
module RollingMill

const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)

    module DataModel

    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)

    include("rollingmill_model_pb.jl")

    end

end


text:
module ConstrainedProblems

import ..OptimizationProblem
import ..ConstraintData

examples = Dict{AbstractString, OptimizationProblem}()

hs9_obj(x::AbstractVector) = sin(π*x[1]/12) * cos(π*x[2]/16)
hs9_c!(c::AbstractVector, x::AbstractVector) = (c[1] = 4*x[1]-3*x[2]; c)
hs9_h!(h, x, λ) = h

function hs9_obj_g!(g::AbstractVector, x::AbstractVector)
    g[1] = π/12 * cos(π*x[1]/12) * cos(π*x[2]/16)
    g[2] = -π/16 * sin(π*x[1]/12) * sin(π*x[2]/16)
    g
end
function hs9_obj_h!(h::AbstractMatrix, x::AbstractVector)
    v = hs9_obj(x)
    h[1,1] = -π^2*v/144
    h[2,2] = -π^2*v/256
    h[1,2] = h[2,1] = -π^2 * cos(π*x[1]/12) * sin(π*x[2]/16) / 192
    h
end

function hs9_jacobian!(J, x)
    J[1,1] = 4
    J[1,2] = -3
    J
end

# TODO: IPNewtons  gets stuck when using x0 = [0,0].
#       Check with Tim if this also happened before?
examples["HS9"] = OptimizationProblem("HS9",
                                      hs9_obj,
                                      hs9_obj_g!,
                                      nothing,
                                      hs9_obj_h!,
                                      ConstraintData(hs9_c!, hs9_jacobian!, hs9_h!,
                                                     [], [], [0.0], [0.0]),
                                      [-1.0,2.0],#[0.0, 0.0],
                                      [-3.0,-4.0],#[[12k-3, 16k-4] for k in (0, 1, -1)], # any integer k will do...
                                      hs9_obj([-3.0,-4.0]),
                                      true,
                                      true)

# Hock and Schittkowski problem number 39.
#
#   Source:
#   Problem 39 in
#   W. Hock and K. Schittkowski,
#   Test examples for nonlinear programming codes,
#   Lectures Notes in Economics and Mathematical Systems 187,
#   Springer Verlag, Heidelberg, 1981.
#
#   classification LOR2-AN-4-2
#
#
# Original source:
# Example 8.7 in
# A. Miele, J.L Tiete, and A.V. Levy
# COMPARISON OF SEVERAL GRADIENT ALGORITHMS FOR MATHEMATICAL PROGRAMMING PROBLEMS, 1972
#

hs39_obj(x) = -x[1]
function hs39_obj_g!(g::AbstractVector{T}, x) where T
    g[1] = -T(1)
    g[2:4] .= T(0)
    g
end
function hs39_obj_h!(h::AbstractMatrix{T}, x) where T
    h .= T(0)
    h
end

function hs39_c!(c, x)
    c[1] = x[2] - x[1]^3 - x[3]^2
    c[2] = x[1]^2 - x[2] - x[4]^2
    c
end

function hs39_jacobian!(J::AbstractMatrix{T}, x) where T
    J[1,1] = -3x[1]^2
    J[1,2] = T(1)
    J[1,3] = -2x[3]

    J[2,1] = 2x[1]
    J[2,2] = -T(1)
    J[2,4] = -2x[4]

    J
end

function hs39_h!(h::AbstractMatrix{T}, x, λ) where T
    h[1,1] += -λ[1]*6*x[1]
    h[3,3] += -λ[1]*2

    h[1,1] += λ[2]*2
    h[4,4] += -λ[2]*2

    h
end

examples["HS39"] = OptimizationProblem("HS39",
                                       hs39_obj,
                                       hs39_obj_g!,
                                       nothing,
                                       hs39_obj_h!,
                                       ConstraintData(hs39_c!, hs39_jacobian!, hs39_h!,
                                                      [], [], [0.0,0.0], [0.0,0.0]),
                                       fill(2.0,4),
                                       [1.0, 1.0, 0, 0],
                                       hs39_obj([1.0, 1.0, 0, 0]),
                                       true,
                                       true)

end  # module


text:
@testset "deprecations" begin
    @testset "propagate" begin
        struct GCN{A<:AbstractMatrix, B, F} <: GNNLayer
            weight::A
            bias::B
            σ::F
        end

        Flux.@functor GCN # allow collecting params, gpu movement, etc...

        function GCN(ch::Pair{Int,Int}, σ=identity)
            in, out = ch
            W = Flux.glorot_uniform(out, in)
            b = zeros(Float32, out)
            GCN(W, b, σ)
        end

        GraphNeuralNetworks.compute_message(l::GCN, xi, xj, e) = xj

        function (l::GCN)(g::GNNGraph, x::AbstractMatrix{T}) where T
            x, _ = propagate(l, g, +, x)
            return l.σ.(l.weight * x .+ l.bias)
        end

        function new_forward(l, g, x)
            x = propagate(copy_xj, g, +, xj=x)
            return l.σ.(l.weight * x .+ l.bias)
        end

        g = GNNGraph(random_regular_graph(10, 4), ndata=randn(3, 10))
        l = GCN(3 => 5, tanh)
        @test l(g, g.ndata.x) ≈ new_forward(l, g, g.ndata.x)
    end
end


text:
using SafeTestsets

@safetestset "Utilities" begin include("utils.jl") end
@safetestset "Index Reduction & SCC" begin include("index_reduction.jl") end
@safetestset "Tearing" begin include("tearing.jl") end
@safetestset "Bareiss" begin include("bareiss.jl") end


text:
using julia_try2_package
using Test

@testset "julia_try2_package.jl" begin
    # Write your tests here.
end


text: const BPF_F_FAST_STACK_CMP = UInt64(1 << 9)


text:
#===============================================================================

            Utility struct for saving most recently sampled paths.

===============================================================================#

mutable struct CyclicCounter
    first::Int64
    last::Int64
    i::Int64

    function CyclicCounter(n; first=1, start_from=first)
        new(first, n, start_from)
    end
end

function next!(cc::CyclicCounter)
    cc.i = mod1(cc.i+1, cc.last)
    cc.i < cc.first && (cc.i = cc.first)
end

(cc::CyclicCounter)() = cc.i
function (cc::CyclicCounter)(::Val{:next})
    i = cc.i
    next!(cc)
    i
end

struct PathSavingBuffer{TX}
    XX::TX
    iter::CyclicCounter

    function PathSavingBuffer(XX, N)
        XX_container = map(1:length(XX)) do i
            trajectory(
                glue_containers( map(x->x.t, XX[i]) ),
                glue_containers( map(x->x.x, XX[i]) ),
            )
        end
        XX_containers = [deepcopy(XX_container) for _ in 1:N]

        new{typeof(XX_containers)}(XX_containers, CyclicCounter(N))
    end
end

function glue_containers(xs)
    glued_xs = collect(Iterators.flatten(map(x->x[1:end-1], xs)))
    append!(glued_xs, [xs[end][end]])
    glued_xs
end

save!(buf::PathSavingBuffer, XX) = _copy_path!(buf.XX[buf.iter(_NEXT)], XX)

function _copy_path!(_to, _from)
    for i in eachindex(_to, _from)
        offset = 1
        for j in eachindex(_from[i])
            N = length(_from[i][j].x)
            view(_to[i].x, offset:offset+N-2) .= view(_from[i][j].x, 1:N-1)
            offset += N-1
        end
        _to[i].x[end] = _from[i][end].x[end]
    end
end


text:
# Autogenerated wrapper script for SCOTCH_jll for i686-linux-musl
export libscotch, libscotcherr, libscotcherrexit, libscotchmetis

using Zlib_jll
JLLWrappers.@generate_wrapper_header("SCOTCH")
JLLWrappers.@declare_library_product(libscotch, "libscotch.so")
JLLWrappers.@declare_library_product(libscotcherr, "libscotcherr.so")
JLLWrappers.@declare_library_product(libscotcherrexit, "libscotcherrexit.so")
JLLWrappers.@declare_library_product(libscotchmetis, "libscotchmetis.so")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libscotch,
        "lib/libscotch.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscotcherr,
        "lib/libscotcherr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscotcherrexit,
        "lib/libscotcherrexit.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscotchmetis,
        "lib/libscotchmetis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@par function cross!(outx::AbstractArray{T},outy::AbstractArray{T},outz::AbstractArray{T},
                     ux::AbstractArray{T},uy::AbstractArray{T},uz::AbstractArray{T},
                     vx::AbstractArray{T},vy::AbstractArray{T},vz::AbstractArray{T},
                     s::@par(AbstractSimulation)) where {T<:Complex}
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                outx[i,j,k] = uy[i,j,k]*vz[i,j,k] - uz[i,j,k]*vy[i,j,k]
                outy[i,j,k] = uz[i,j,k]*vx[i,j,k] - ux[i,j,k]*vz[i,j,k]
                outz[i,j,k] = ux[i,j,k]*vy[i,j,k] - uy[i,j,k]*vx[i,j,k]
            end
    end
end

@par function crossk!(outx,outy,outz,vx,vy,vz,s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                outx[i,j,k] = im*(KY[j]*vz[i,j,k] - KZ[k]*vy[i,j,k])
                outy[i,j,k] = im*(KZ[k]*vx[i,j,k] - KX[i]*vz[i,j,k])
                outz[i,j,k] = im*(KX[i]*vy[i,j,k] - KY[j]*vx[i,j,k])
            end
    end
end

function ccross!(out::VectorField,u::VectorField,v::VectorField,s::AbstractSimulation)
    cross!(out.cx,out.cy,out.cz,u.cx,u.cy,u.cz,v.cx,v.cy,v.cz,s)
    return out
end

function curl!(out::VectorField,u::VectorField,s::AbstractSimulation)
    crossk!(out.cx,out.cy,out.cz,u.cx,u.cy,u.cz,s)
    return out
end

@par function grad!(out::VectorField,f::AbstractArray{<:Complex,3},s::@par(AbstractSimulation))
    outx = out.c.x
    outy = out.c.y
    outz = out.c.z
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                outx[i,j,k] = f[i,j,k]*im*KX[i]
                outy[i,j,k] = f[i,j,k]*im*KY[j]
                outz[i,j,k] = f[i,j,k]*im*KZ[k]
            end
    end
    #dx!(out.cx,f,s)
    #dy!(out.cy,f,s)
    #dz!(out.cz,f,s)
    #dealias!(out,s)
end

@par function dx!(out::AbstractArray{<:Complex,3},f::AbstractArray{<:Complex,3},s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                out[i,j,k] = f[i,j,k]*im*KX[i]
            end
    end
end

@par function dy!(out::AbstractArray{<:Complex,3},f::AbstractArray{<:Complex,3},s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                out[i,j,k] = f[i,j,k]*im*KY[j]
            end
    end
end

@par function dz!(out::AbstractArray{<:Complex,3},f::AbstractArray{<:Complex,3},s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                out[i,j,k] = f[i,j,k]*im*KZ[k]
            end
    end
end

@par function div!(out::AbstractArray{Complex{Float64},3},ux,uy,uz,w,mdρdz::Real, s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
               # out[i,j,k] = mim*(kx[i]*ux[i,j,k] + ky[j]*uy[i,j,k] + kz[k]*uz[i,j,k]) + mdœÅdz*w[i,j,k]
                out[i,j,k] = muladd(im,muladd(KX[i], ux[i,j,k], muladd(KY[j], uy[i,j,k], KZ[k]*uz[i,j,k])), mdρdz*w[i,j,k])
            end
    end
end

div!(out::AbstractArray{<:Complex,3},u::VectorField,s) = div!(out,u.cx,u.cy,u.cz,s)

@par function div!(out::AbstractArray{Complex{Float64},3},ux,uy,uz, s::@par(AbstractSimulation))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
            @inbounds @msimd for i in XRANGE
                #out[i,j,k] = -im*(kx[i]*ux[i,j,k] + ky[j]*uy[i,j,k] + kz[k]*uz[i,j,k])
                out[i,j,k] = im*muladd(KX[i], ux[i,j,k], muladd(KY[j], uy[i,j,k], KZ[k]*uz[i,j,k]))
            end
    end
end


text:
#=======================================================================
Convenience methods related to construction of pipeline::StanfordCoreNLP.

Usage: pline = pipeline(annotators::String...)

Example: pline = pipeline("pos", "dcoref")

Uses `check_annotators` to validate the annotators in the argument of pipeline,
i.e., it checks that each annotator is a key in ANNOTATORS. ANNOTATORS must be
located in ./data/annotators.jl.

Uses `annotator_list` to build a single, comma-separated string of the
annotators passed to `pipeline` and their dependencies. For example,
annotator_list("pos") returns string "tokenize, ssplit, pos". Annotators
listed in the returned string are unique and ordered: each
annotator is to the right of its dependencies.
=======================================================================#

#Construct a CoreNLP pipeline object suitable for `annotate!` calls.
#=
function pipeline(annotators::String...)
    if check_annotators(annotators...)
        annList = annotator_list(annotators...)
        println("annList = ", annList)  #diagnostic
        return StanfordCoreNLP(Dict("annotations" => annList))
    else
        println("Invalid annotator in `pipeline` call.")
        return 1
    end
end

function pipeline(annotators::String...)
    if check_annotators(annotators...)
        annList = annotator_list(annotators...)
        println("annList = ", annList)  #diagnostic
        jprops = JProperties()
        jcall(jprops, "setProperty", JObject, (JString, JString), "annotators", annList)
        return StanfordCoreNLP(JStanfordCoreNLP((JProperties,), jprops))
    else
        println("Invalid annotator in `pipeline` call.")
        return 1
    end
end

#Returns true if all arguments are valid annotators, false otherwise.
include("./data/annotators.jl")
function check_annotators(annotators...)
    valid = true
        for a in annotators
            valid = valid && in(a, keys(ANNOTATORS))
        end
    return valid
end

#Accepts an unordered list of annotators and returns an ordered, comma-separated
#string of annotators and their dependencies. In the returned string, each
# annotator is to the right of its dependencies. For example,
#annotatorList("pos", "ner") returns "tokenize, ssplit, pos, lemma, ner".
#The return value is intended as the argument of Pipeline(...)
function annotator_list(annotators...)
    ##Push annotators onto a list. Make unique and sort.
    list = String[]
    for a in annotators
        for b in ANNOTATORS[a]
            push!(list, b)
        end
        push!(list, a)
    end
    list = unique(list)

    ##Construct Pipeline() call argument by stringing annotators together.
    arg = ""
    for a in list
       arg = Base.string(arg, a, ",")
    end
#    arg = arg[1:length(arg)-1] #Drop the trailing comma.
    pop!(arg)  #Drop the trailing comma.
end

#Annotates doc.jann according to pipeline. Note, doc is an Annotation
#whereas doc.jann is a JAnnotation.
function annotate!(pipeline::JavaCoreNLP.StanfordCoreNLP, doc::Annotation)   #??
    jcall(pipeline.jpipeline, "annotate", Void, (JAnnotation,), doc.jann)
end
=#


text:
## Types ##

struct ModelState{T<:AbstractArray, N<:Integer}
    Nx::N
    Ny::N
    dtype::Type

    # State variables
    state::OrderedDict{Symbol, T}

    # time increments
    inc::OrderedDict{Symbol, T}

    # diagnostic variables
    diagnostic::OrderedDict{Symbol, T}

    # inner constructor
    function ModelState{T,N}(Nx::N, Ny::N) where {T<:AbstractArray, N<:Integer}
        d = OrderedDict{Symbol, T}()
        new{T, N}(Nx, Ny, eltype(T), d, copy(d), copy(d))
    end
end
ModelState{T}(Nx::N, Ny::N) where {T<:AbstractArray, N<:Integer} = ModelState{T, N}(Nx, Ny)

function add_variable!(s::MS, reg::Symbol, var::Symbol) where {MS<:ModelState}
    getproperty(s, reg)[var] = zeros(eltype(valtype(getproperty(s, reg))), s.Nx, s.Ny)
end

function add_variable!(s::MS, eq::PrognosticEquation) where {MS<:ModelState}
    add_variable!(s, :state, eq.var)
    add_variable!(s, :inc, eq.var)
    add_variable!(s, :inc, np1_symbol(eq.var))
end

add_variable!(s::MS, eq::DiagnosticEquation) where {MS<:ModelState} = add_variable!(s, :diagnostic, eq.var)

@inline np1_symbol(s::Symbol) = Symbol(s, "_next")

struct Model{S<:ModelState}
    # model state
    s::S

    # set of equations
    equations::Tuple

    # scratch workspace
    scratch::Vector

    function Model(state::S, eqs::Tuple) where {S<:ModelState}
        for eq in eqs
            add_variable!(state, eq)
        end
        new{S}(state, eqs, [])
    end
end

## Model mechanics ##

function create_scratch(m::M) where {M<:Model}
    template_var = first(keys(m.s.state))
    push!(m.scratch, similar(m.s.state[template_var]))
end
n_scratch(m::M) where M<:Model = length(m.scratch)

function guarantee_n_scratch!(m::M, n::I) where {M<:Model, I<:Integer}
    if n_scratch(m) >= n
        return
    end
    for i in 1:n - n_scratch(m)
        create_scratch(m)
    end
end

@inline function step!(m::M) where {M<:Model}
    # setup workspace
    guarantee_n_scratch!(m, 1)
    scratch = m.scratch[1]

    # compute diagnostic variables
    for eq in m.equations
        is_diagnostic(eq) ? step!(eq, m.s, scratch::T) : continue
    end

    # prognostic variables increment
    for eq in m.equations
        is_prognostic(eq) ? step!(eq, m.s, scratch::T) : continue
    end
end

@inline step!(eq::DiagnosticEquation, s::S, scratch::T) where {E<:Equation, S<:ModelState, T<:AbstractArray} =
    step!(eq.var, eq.terms, s, :diagnostic, scratch)

@inline step!(eq::PrognosticEquation, s::S, scratch::T) where {S<:ModelState, T<:AbstractArray} =
    step!(eq.var, eq.terms, s, :inc, scratch)

@inline function step!(var::Symbol, terms::Tuple, s::S, reg::Symbol, scratch::T) where {S<:ModelState, T<:AbstractArray}
    if length(terms) == 0
        return
    end
    v_arr = getproperty(s, reg)[var]
    scratch .= 0
    for (i, term!) in enumerate(terms)
        if i == 1
            term!(v_arr, s)
        else
            term!(scratch, s)
            v_arr .+= scratch
        end
    end
end

@inline function advance!(m::M) where M<:Model
    map(advance!, x -> advance!(eq, m.s), m.equations)
end
@inline function advance!(eq::PrognosticEquation, s::S) where {S<:ModelState}
    s.state[eq.var] .= s.inc[eq.next];
end
@inline function advance!(eq::E, s::S) where {E<:Equation, S<:ModelState} end

## Model iteration (the full workflow of one iteration)

function iterate!(m::M, dt::R, nt::N) where {M<:Model, R<:Real, N<:Integer}
    for i in 1:nt
        iterate!(m, dt)
    end
end
@inline function iterate!(m::M, dt::R) where {M<:Model, R<:Real}
    step!(m)
    integrate!(m, dt)
    advance!(m)
end
@inline function iterate!(eq::E, s::S, scratch::T, dt::N) where {E<:Equation, S<:ModelState, T<:AbstractArray, N<:Real}
    step!(eq, s, scratch)
    integrate!(eq, s, dt)
    advance!(eq, s)
end
@inline function iterate_heaps!(m::M, dt::R) where {M<:Model, R<:Real}
    guarantee_n_scratch!(m, 1)
    scratch = m.scratch[1]

    for eq in m.equations
        iterate!(eq, m.s, scratch, dt)
    end
end

function iterate_heaps!(m::M, dt::R, nt::N) where {M<:Model, R<:Real, N<:Integer}
    for i in 1:nt
        iterate_heaps!(m, dt)
    end
end


text:
push!(LOAD_PATH, "../src/")
ENV["GKSwstype"] = "100"
using Documenter
using GEMPIC
using Literate
using Plots

# generate examples
# output = joinpath(@__DIR__, "src", "generated")
# examples = String[]
# push!(examples, "strong_landau_damping.jl")
#
# for example in examples
#     jl_file = joinpath(@__DIR__, "examples", example)
#     Literate.markdown(jl_file, output, documenter=true)
#     #Literate.notebook(EXAMPLE, OUTPUT)
#     #Literate.script(EXAMPLE, OUTPUT)
# end

makedocs(;
    sitename="GEMPIC",
    doctest=true,
    authors="Julia Vlasov",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        mathengine=MathJax(
            Dict(
                :TeX =>
                    Dict(:equationNumbers => Dict(:autoNumber => "AMS"), :Macros => Dict()),
            ),
        ),
    ),
    modules=[GEMPIC],
    pages=[
        "Documentation" => [
            "index.md",
            "mesh.md",
            "low_level_bsplines.md",
            "splinepp.md",
            "particle_mesh_coupling.md",
            "distributions.md",
            "particle_group.md",
            "particle_sampling.md",
            "maxwell_solver.md",
            "Splitting" => [
                "hamiltonian_splitting.md",
                "hamiltonian_splitting_spin.md",
                "hamiltonian_splitting_boris.md",
            ],
            "diagnostics.md",
        ],
        "Examples" => ["strong_landau_damping.md"],
        # "Scalar Spin Vlasov-Maxwell" => "scalar_spin_vlasov_maxwell.md",
        "Contents" => "contents.md",
    ],
)

deploydocs(; repo="github.com/JuliaVlasov/GEMPIC.jl.git")


text:
using MathOptInterface
const MOI = MathOptInterface
const AFF = MOI.ScalarAffineFunction{Cdouble}
const EQ = MOI.EqualTo{Cdouble}
const AFFEQ = MOI.ConstraintIndex{AFF,EQ}

mutable struct Optimizer <: MOI.AbstractOptimizer
    objective_constant::Cdouble
    objective_sign::Int
    blockdims::Vector{CSDP_INT}
    varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
    num_entries::Dict{Tuple{Int, Int}, Int}
    b::Vector{Cdouble}
    C::blockmatrix
    problem::Union{Nothing, LoadingProblem}
    X::blockmatrix
    y::Union{Nothing, Vector{Cdouble}}
    Z::blockmatrix
    status::CSDP_INT
    pobj::Cdouble
    dobj::Cdouble
    solve_time::Float64
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
        optimizer = new(
            zero(Cdouble), 1, CSDP_INT[], Tuple{Int, Int, Int}[],
            Dict{Tuple{Int, Int}, Int}(), Cdouble[],
            blockmatrix(), nothing, blockmatrix(), nothing, blockmatrix(),
            -1, NaN, NaN, NaN, false, Dict{Symbol, Any}())
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawOptimizerAttribute(String(key)), value)
        end
        # May need to call `free_loaded_prob` and `free_loading_prob`.
        finalizer(MOI.empty!, optimizer)
        return optimizer
    end
end

varmap(optimizer::Optimizer, vi::MOI.VariableIndex) = optimizer.varmap[vi.value]

function MOI.supports(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    return Symbol(param.name) in ALLOWED_OPTIONS
end
function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value)
    if !(param.name isa String)
        Base.depwarn(
            "passing `$(param.name)` to `MOI.RawOptimizerAttribute` as type " *
            "`$(typeof(param.name))` is deprecated. Use a string instead.",
            Symbol("MOI.set")
        )
    end
    if !MOI.supports(optimizer, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    optimizer.options[Symbol(param.name)] = value
end
function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO: This gives a poor error message if the name of the parameter is invalid.
    if !(param.name isa String)
        Base.depwarn(
            "passing `$(param.name)` to `MOI.RawOptimizerAttribute` as type " *
            "`$(typeof(param.name))` is deprecated. Use a string instead.",
            Symbol("MOI.set")
        )
    end
    return optimizer.options[Symbol(param.name)]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.get(::Optimizer, ::MOI.SolverName) = "CSDP"
# See table "Return codes for easy_sdp() and CSDP" in `doc/csdpuser.pdf`.
const RAW_STATUS = [
    "Problem solved to optimality.",
    "Problem is primal infeasible.",
    "Problem is dual infeasible.",
    "Problem solved to near optimality.",
    "Maximum iterations reached.",
    "Stuck at edge of primal feasibility.",
    "Stuck at edge of dual feasibility.",
    "Lack of progress.",
    "X, Z, or O is singular.",
    "NaN or Inf values encountered.",
    "Program stopped by signal (SIXCPU, SIGTERM, etc.)"]

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return RAW_STATUS[optimizer.status + 1]
end
function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec)
    return optimizer.solve_time
end

function MOI.is_empty(optimizer::Optimizer)
    return iszero(optimizer.objective_constant) &&
        isone(optimizer.objective_sign) &&
        isempty(optimizer.blockdims) &&
        isempty(optimizer.varmap) &&
        isempty(optimizer.num_entries) &&
        isempty(optimizer.b) &&
        iszero(optimizer.C.nblocks) &&
        optimizer.C.blocks == C_NULL &&
        optimizer.problem === nothing
end

function MOI.empty!(optimizer::Optimizer)
    optimizer.objective_constant = zero(Cdouble)
    optimizer.objective_sign = 1
    empty!(optimizer.blockdims)
    empty!(optimizer.varmap)
    empty!(optimizer.num_entries)
    empty!(optimizer.b)
    if optimizer.problem !== nothing
        if optimizer.y !== nothing
            free_loaded_prob(optimizer.problem, optimizer.X, optimizer.y, optimizer.Z)
        end
        free_loading_prob(optimizer.problem)
    end
    optimizer.problem = nothing
    optimizer.C.nblocks = 0
    optimizer.C.blocks = C_NULL
    optimizer.X.nblocks = 0
    optimizer.X.blocks = C_NULL
    optimizer.y = nothing
    optimizer.Z.nblocks = 0
    optimizer.Z.blocks = C_NULL
    optimizer.status = -1
    optimizer.pobj = 0.0
    optimizer.dobj = 0.0
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{AFF}})
    return true
end

MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false

const SupportedSets = Union{MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}
MOI.supports_add_constrained_variables(::Optimizer, ::Type{<:SupportedSets}) = true
function MOI.supports_constraint(::Optimizer, ::Type{AFF}, ::Type{EQ})
    return true
end

function _new_block(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.blockdims, -MOI.dimension(set))
    blk = length(optimizer.blockdims)
    for i in 1:MOI.dimension(set)
        push!(optimizer.varmap, (blk, i, i))
    end
end

function _new_block(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.blockdims, set.side_dimension)
    blk = length(optimizer.blockdims)
    for i in 1:set.side_dimension
        for j in 1:i
            push!(optimizer.varmap, (blk, i, j))
        end
    end
end

function _add_constrained_variables(optimizer::Optimizer, set::SupportedSets)
    offset = length(optimizer.varmap)
    _new_block(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function _error(start, stop)
    error(start, ". Use `MOI.instantiate(CSDP.Optimizer, with_bridge_type = Float64)` ", stop)
end

function constrain_variables_on_creation(
    dest::MOI.ModelLike,
    src::MOI.ModelLike,
    index_map::MOI.Utilities.IndexMap,
    ::Type{S},
) where {S<:MOI.AbstractVectorSet}
    for ci_src in
        MOI.get(src, MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}())
        f_src = MOI.get(src, MOI.ConstraintFunction(), ci_src)
        if !allunique(f_src.variables)
            _error("Cannot copy constraint `$(ci_src)` as variables constrained on creation because there are duplicate variables in the function `$(f_src)`",
                   "to bridge this by creating slack variables.")
        elseif any(vi -> haskey(index_map, vi), f_src.variables)
            _error("Cannot copy constraint `$(ci_src)` as variables constrained on creation because some variables of the function `$(f_src)` are in another constraint as well.",
                   "to bridge constraints having the same variables by creating slack variables.")
        else
            set = MOI.get(src, MOI.ConstraintSet(), ci_src)::S
            vis_dest, ci_dest = _add_constrained_variables(dest, set)
            index_map[ci_src] = ci_dest
            for (vi_src, vi_dest) in zip(f_src.variables, vis_dest)
                index_map[vi_src] = vi_dest
            end
        end
    end
end

function count_entry(optimizer::Optimizer, con_idx::Integer, blk::Integer)
    key = (con_idx, blk)
    optimizer.num_entries[key] = get(optimizer.num_entries, key, 0) + 1
end

# Loads objective coefficient α * vi
function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    # in SDP format, it is max and in MPB Conic format it is min
    coef = optimizer.objective_sign * α
    if i != j
        coef /= 2
    end
    addentry(optimizer.problem, 0, blk, i, j, coef, true)
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    MOI.empty!(dest)
    index_map = MOI.Utilities.IndexMap()

    # Step 1) Compute the dimensions of what needs to be allocated
    constrain_variables_on_creation(
        dest,
        src,
        index_map,
        MOI.Nonnegatives,
    )
    constrain_variables_on_creation(
        dest,
        src,
        index_map,
        MOI.PositiveSemidefiniteConeTriangle,
    )
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    if length(vis_src) < length(index_map.var_map)
        _error("Free variables are not supported by CSDP",
               "to bridge free variables into `x - y` where `x` and `y` are nonnegative.")
    end
    cis_src = MOI.get(src, MOI.ListOfConstraintIndices{AFF,EQ}())
    dest.b = Vector{Cdouble}(undef, length(cis_src))
    funcs = Vector{AFF}(undef, length(cis_src))
    for (k, ci_src) in enumerate(cis_src)
        funcs[k] = MOI.get(src, MOI.CanonicalConstraintFunction(), ci_src)
        set = MOI.get(src, MOI.ConstraintSet(), ci_src)
        if isempty(funcs[k].terms)
            throw(ArgumentError("Empty constraint $cis_src: $(funcs[k])-in-$set. Not supported by CSDP."))
        end
        if !iszero(MOI.constant(funcs[k]))
            throw(MOI.ScalarFunctionConstantNotZero{
                Cdouble, AFF, EQ}(
                    MOI.constant(funcs[k])))
        end
        for t in funcs[k].terms
            if !iszero(t.coefficient)
                blk, _, _ = varmap(dest, index_map[t.variable])
                count_entry(dest, k, blk)
            end
        end
        dest.b[k] = MOI.constant(set)
        index_map[ci_src] = AFFEQ(k)
    end

    # Step 2) Allocate CSDP datastructures
    dummy = isempty(dest.b)
    if dummy
        # See https://github.com/coin-or/Csdp/issues/2
        dest.b = [one(Cdouble)]
        dest.blockdims = [dest.blockdims; CSDP_INT(-1)]
        count_entry(dest, 1, length(dest.blockdims))
    end
    dest.C.nblocks = length(dest.blockdims)
    num_entries = zeros(CSDP_INT, length(dest.b), length(dest.blockdims))
    for (key, value) in dest.num_entries
        num_entries[key...] = value
    end
    dest.problem = allocate_loading_prob(Ref(dest.C), dest.blockdims, length(dest.b), num_entries, 3)
    if dummy
        # See https://github.com/coin-or/Csdp/issues/2
        duplicate = addentry(dest.problem, 1, length(dest.blockdims), 1, 1, 1.0, true)
        @assert !duplicate
    end

    # Step 3) Load data in the datastructures
    for k in eachindex(funcs)
        setconstant(dest.problem, k, dest.b[k])
        for term in funcs[k].terms
            if !iszero(term.coefficient)
                blk, i, j = varmap(dest, index_map[term.variable])
                coef = term.coefficient
                if i != j
                    coef /= 2
                end
                duplicate = addentry(dest.problem, k, blk, i, j, coef, true)
                @assert !duplicate
            end
        end
    end

    # Throw error for variable attributes
    MOI.Utilities.pass_attributes(dest, src, index_map, vis_src)
    # Throw error for constraint attributes
    MOI.Utilities.pass_attributes(dest, src, index_map, cis_src)

    # Pass objective attributes and throw error for other ones
    model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
    for attr in model_attributes
        if attr != MOI.ObjectiveSense() && attr != MOI.ObjectiveFunction{AFF}()
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    # We make sure to set `objective_sign` first before setting the objective
    if MOI.ObjectiveSense() in model_attributes
        sense = MOI.get(src, MOI.ObjectiveSense())
        dest.objective_sign = sense == MOI.MIN_SENSE ? -1 : 1
    end
    if MOI.ObjectiveFunction{AFF}() in model_attributes
        func = MOI.get(src, MOI.ObjectiveFunction{AFF}())
        obj = MOI.Utilities.canonical(func)
        dest.objective_constant = obj.constant
        for term in obj.terms
            if !iszero(term.coefficient)
                load_objective_term!(dest, term.coefficient, index_map[term.variable])
            end
        end
    end
    return index_map
end

function MOI.optimize!(optimizer::Optimizer)
    write_prob(optimizer)

    start_time = time()
    optimizer.y = loaded_initsoln(optimizer.problem, length(optimizer.b), Ref(optimizer.X), Ref(optimizer.Z))

    options = optimizer.options
    if optimizer.silent
        options = copy(options)
        options[:printlevel] = 0
    end

    optimizer.status, optimizer.pobj, optimizer.dobj = loaded_sdp(
        optimizer.problem, optimizer.objective_sign * optimizer.objective_constant,
        Ref(optimizer.X), optimizer.y, Ref(optimizer.Z), options,
    )
    optimizer.solve_time = time() - start_time
    return
end

function MOI.get(m::Optimizer, ::MOI.TerminationStatus)
    status = m.status
    if status == -1
        return MOI.OPTIMIZE_NOT_CALLED
    elseif status == 0
        return MOI.OPTIMAL
    elseif status == 1
        return MOI.INFEASIBLE
    elseif status == 2
        return MOI.DUAL_INFEASIBLE
    elseif status == 3
        return MOI.ALMOST_OPTIMAL
    elseif status == 4
        return MOI.ITERATION_LIMIT
    elseif 5 <= status <= 7
        return MOI.SLOW_PROGRESS
    elseif 8 <= status <= 9
        return MOI.NUMERICAL_ERROR
    else
        error("Internal library error: status=$status")
    end
end

function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = m.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 1
        return MOI.INFEASIBLE_POINT
    elseif status == 2
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == 3
        return MOI.NEARLY_FEASIBLE_POINT
    elseif 4 <= status <= 9
        return MOI.UNKNOWN_RESULT_STATUS
    else
        error("Internal library error: status=$status")
    end
end

function MOI.get(m::Optimizer, attr::MOI.DualStatus)
    if attr.result_index > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    status = m.status
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif status == 1
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == 2
        return MOI.INFEASIBLE_POINT
    elseif status == 3
        return MOI.NEARLY_FEASIBLE_POINT
    elseif 4 <= status <= 9
        return MOI.UNKNOWN_RESULT_STATUS
    else
        error("Internal library error: status=$status")
    end
end

MOI.get(m::Optimizer, ::MOI.ResultCount) = 1
function MOI.get(m::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objective_sign * m.pobj
end
function MOI.get(m::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return m.objective_sign * m.dobj
end
struct PrimalSolutionMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::PrimalSolutionMatrix) = true
MOI.get(optimizer::Optimizer, ::PrimalSolutionMatrix) = optimizer.X

struct DualSolutionVector <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSolutionVector) = true
MOI.get(optimizer::Optimizer, ::DualSolutionVector) = optimizer.y

struct DualSlackMatrix <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::DualSlackMatrix) = true
MOI.get(optimizer::Optimizer, ::DualSlackMatrix) = optimizer.Z

function block(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    return optimizer.varmap[ci.value][1]
end
function dimension(optimizer::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    blockdim = optimizer.blockdims[block(optimizer, ci)]
    if blockdim < 0
        return -blockdim
    else
        return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(blockdim))
    end
end
function vectorize_block(M, blk::Integer, s::Type{MOI.Nonnegatives})
    return diag(block(M, blk))
end
function vectorize_block(M::AbstractMatrix{Cdouble}, blk::Integer, s::Type{MOI.PositiveSemidefiniteConeTriangle}) where T
    B = block(M, blk)
    d = LinearAlgebra.checksquare(B)
    n = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))
    v = Vector{Cdouble}(undef, n)
    k = 0
    for j in 1:d
        for i in 1:j
            k += 1
            v[k] = B[i, j]
        end
    end
    @assert k == n
    return v
end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(optimizer, attr)
    blk, i, j = varmap(optimizer, vi)
    return block(MOI.get(optimizer, PrimalSolutionMatrix()), blk)[i, j]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, PrimalSolutionMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return optimizer.b[ci.value]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
    MOI.check_result_index_bounds(optimizer, attr)
    return vectorize_block(MOI.get(optimizer, DualSlackMatrix()), block(optimizer, ci), S)
end
function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::AFFEQ)
    MOI.check_result_index_bounds(optimizer, attr)
    return -MOI.get(optimizer, DualSolutionVector())[ci.value]
end


text:
module ReactiveMPModelsAutoregressiveTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra, StableRNGs

## Model definition
## -------------------------------------------- ##
@model [ default_factorisation = MeanField() ] function ar_model(n, order)

    x = datavar(Vector{Float64}, n)
    y = datavar(Float64, n)

    γ ~ GammaShapeRate(1.0, 1.0)
    θ ~ MvNormalMeanPrecision(zeros(order), diageye(order))

    for i in 1:n
        y[i] ~ NormalMeanPrecision(dot(x[i], θ), γ)
    end

    return x, y, θ, γ
end

@model [ default_factorisation = MeanField() ] function lar_model(T::Type{ Multivariate }, n, order, c, stype, τ)

    # Parameter priors
    γ  ~ GammaShapeRate(1.0, 1.0)
    θ  ~ MvNormalMeanPrecision(zeros(order), diageye(order))

    # We create a sequence of random variables for hidden states
    x = randomvar(n)
    # As well a sequence of observartions
    y = datavar(Float64, n)

    ct = constvar(c)
    # We assume observation noise to be known
    cτ = constvar(τ)

    # Prior for first state
    x0 ~ MvNormalMeanPrecision(zeros(order), diageye(order))

    x_prev = x0

    # AR process requires extra meta information
    meta = ARMeta(Multivariate, order, stype)

    for i in 1:n
        # Autoregressive node uses structured factorisation assumption between states
        x[i] ~ AR(x_prev, θ, γ) where { q = q(y, x)q(γ)q(θ), meta = meta }
        y[i] ~ NormalMeanPrecision(dot(ct, x[i]), cτ)
        x_prev = x[i]
    end

    return x, y, θ, γ
end
@model [ default_factorisation = MeanField() ] function lar_model(T::Type{Univariate}, n, order, c, stype, τ)

    # Parameter priors
    γ  ~ GammaShapeRate(1.0, 1.0)
    θ  ~ NormalMeanPrecision(0.0, 1.0)

    # We create a sequence of random variables for hidden states
    x = randomvar(n)
    # As well a sequence of observartions
    y = datavar(Float64, n)

    ct = constvar(c)
    # We assume observation noise to be known
    cτ = constvar(τ)

    # Prior for first state
    x0 ~ NormalMeanPrecision(0.0, 1.0)

    x_prev = x0

    # AR process requires extra meta information
    meta = ARMeta(Univariate, order, stype)

    for i in 1:n
        x[i] ~ AR(x_prev, θ, γ) where { q = q(y, x)q(γ)q(θ), meta = meta }
        y[i] ~ NormalMeanPrecision(ct * x[i], cτ)
        x_prev = x[i]
    end

    return x, y, θ, γ
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function ar_inference(inputs, outputs, order, niter)
    n = length(outputs)

    model, (x, y, θ, γ) = ar_model(n, order, options = (limit_stack_depth = 500, ))

    γ_buffer = keep(Marginal)
    θ_buffer = keep(Marginal)
    fe       = keep(Float64)

    γ_sub = subscribe!(getmarginal(γ), γ_buffer)
    θ_sub = subscribe!(getmarginal(θ), θ_buffer)
    f_sub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    setmarginal!(γ, GammaShapeRate(1.0, 1.0))

    for _ in 1:niter
        update!(x, inputs)
        update!(y, outputs)
    end

    unsubscribe!((γ_sub, θ_sub, f_sub))

    return γ_buffer, θ_buffer, fe
end

function lar_init_marginals!(::Type{ Multivariate }, order, γ, θ)
    setmarginal!(γ, GammaShapeRate(1.0, 1.0))
    setmarginal!(θ, MvNormalMeanPrecision(zeros(order), diageye(order)))
end

function lar_init_marginals!(::Type{ Univariate }, order, γ, θ)
    setmarginal!(γ, GammaShapeRate(1.0, 1.0))
    setmarginal!(θ, NormalMeanPrecision(0.0, 1.0))
end

function lar_inference(data, order, artype, stype, niter, τ)

    # We build a full graph based on nber of observatios
    n = length(data)

    # Depending on the order of AR process `c` is
    # either a nber or a vector
    c = ReactiveMP.ar_unit(artype, order)

    # Note that to run inference for huge model it might be necessary to pass extra
    # options = (limit_stack_depth = 100,) to limit stack depth during recursive inference procedure
    model, (x, y, θ, γ) = lar_model(artype, n, order, c, stype, τ)

    # We are going to keep `γ` and `θ` estimates for all VMP iterations
    # But `buffer` only last posterior estimates for a sequence of hidden states `x`
    # We also will keep Bethe Free Energy in `fe`
    γ_buffer = keep(Marginal)
    θ_buffer = keep(Marginal)
    x_buffer = buffer(Marginal, n)
    fe       = keep(Float64)

    γsub  = subscribe!(getmarginal(γ), γ_buffer)
    θsub  = subscribe!(getmarginal(θ), θ_buffer)
    xsub  = subscribe!(getmarginals(x), x_buffer)
    fesub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    lar_init_marginals!(artype, order, γ, θ)

    # We update data several times to perform several VMP iterations
    for i in 1:niter
        update!(y, data)
    end

    # It is important to unsubscribe from running observables
    unsubscribe!((γsub, θsub, xsub, fesub))

    return γ_buffer, θ_buffer, x_buffer, fe
end

@testset "Autoregressive model" begin

    @testset "Autoregressive model" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        function ar_ssm(series, order)
            inputs = [reverse!(series[1:order])]
            outputs = [series[order + 1]]
            for x in series[order+2:end]
                push!(inputs, vcat(outputs[end], inputs[end])[1:end-1])
                push!(outputs, x)
            end
            return inputs, outputs
        end
        rng  = StableRNG(1234)
        series          = randn(rng, 1_000)
        ## -------------------------------------------- ##
        ## Inference execution and test inference results
        for order in 1:5
            inputs, outputs = ar_ssm(series, order)
            γ_buffer, θ_buffer, fe = ar_inference(inputs, outputs, order, 15)
            @test length(γ_buffer) === 15
            @test length(θ_buffer) === 15
            @test length(fe) === 15
            @test last(fe) < first(fe)
            @test all(filter(e -> abs(e) > 1e-3, diff(getvalues(fe))) .< 0)
        end
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        benchmark_output = joinpath(base_output, "ar_model_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output benchmarks
        inputs5, outputs5 = ar_ssm(series, 5)
        benchmark = @benchmark ar_inference($inputs5, $outputs5, 5, 15)#
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

    @testset "Latent autoregressive model" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        # The following coefficients correspond to stable poles
        coefs_ar_5 = [ 0.10699399235785655, -0.5237303489793305, 0.3068897071844715, -0.17232255282458891, 0.13323964347539288 ]

        function generate_lar_data(rng, n, θ, γ, τ)
            order        = length(θ)
            states       = Vector{Vector{Float64}}(undef, n + 3order)
            observations = Vector{Float64}(undef, n + 3order)

            γ_std = sqrt(inv(γ))
            τ_std = sqrt(inv(γ))

            states[1] = randn(rng, order)

            for i in 2:(n + 3order)
                states[i]       = vcat(rand(rng, Normal(dot(θ, states[i - 1]), γ_std)), states[i-1][1:end-1])
                observations[i] = rand(rng, Normal(states[i][1], τ_std))
            end

            return states[1+3order:end], observations[1+3order:end]
        end
        # Seed for reproducibility
        rng  = StableRNG(123)
        # Number of observations in synthetic dataset
        n = 500
        # AR process parameters
        real_γ = 5.0
        real_τ = 5.0
        real_θ = coefs_ar_5
        states, observations = generate_lar_data(rng, n, real_θ, real_γ, real_τ)
        ## -------------------------------------------- ##
        ## Inference execution

        # AR order 1
        γ, θ, xs, fe = lar_inference(observations, 1, Univariate, ARsafe(), 15, real_τ)
        @test length(xs) === n
        @test length(γ)  === 15
        @test length(θ)  === 15
        @test length(fe) === 15
        @test abs(last(fe) - 518.9182342) < 0.01
        @test all(filter(e -> abs(e) > 1e-3, diff(getvalues(fe))) .< 0)

        for i in 1:4
            γ, θ, xs, fe = lar_inference(observations, i, Multivariate, ARsafe(), 15, real_τ)
            @test length(xs) === n
            @test length(γ)  === 15
            @test length(θ)  === 15
            @test length(fe) === 15
        end

        # AR order 5
        γ, θ, xs, fe = lar_inference(observations, length(real_θ), Multivariate, ARsafe(), 15, real_τ)

        ## -------------------------------------------- ##
        ## Test inference results
        @test length(xs) === n
        @test length(γ)  === 15
        @test length(θ)  === 15
        @test length(fe) === 15
        @test abs(last(fe) - 514.66086) < 0.01
        @test all(filter(e -> abs(e) > 1e-1, diff(getvalues(fe))) .< 0)
        @test (mean(last(γ)) - 3.0std(last(γ)) < real_γ < mean(last(γ)) + 3.0std(last(γ)))
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "lar_model_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "lar_model_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        p1 = plot(first.(states), label="Hidden state")
        p1 = scatter!(p1, observations, label="Observations")
        p1 = plot!(p1, first.(mean.(xs)), ribbon = sqrt.(first.(var.(xs))), label="Inferred states", legend = :bottomright)

        p2 = plot(mean.(γ), ribbon = std.(γ), label = "Inferred transition precision", legend = :bottomright)
        p2 = plot!([ real_γ ], seriestype = :hline, label = "Real transition precision")

        p3 = plot(getvalues(fe), label = "Bethe Free Energy")

        p = plot(p1, p2, p3, layout = @layout([ a; b c ]))
        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark lar_inference($observations, length($real_θ), Multivariate, ARsafe(), 15, $real_τ) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
using PortHamiltonian

# mixed finite elements (Golo)
 Nel = 15; alfa = 1; a = 0; b = 1;
 	p = PortHamiltonian.new_spectral_element_phs_golo(Nel,alfa,a,b)
p = finelem(20,2,0,1);
# pseudo-spectral discretization (Moulla)
#p = discrete_phs(10,0,e-1)
#wi = p.disc_data.flow.wi[:];
#xi = p.disc_data.flow.xi[:];
#Qi = diagm(wi .* (1 .*(1+xi)));
#p.Q = blkdiag(Qi, Qi)
p.GradHam = x-> p.Q *x

# pseudo-spectral "elements"
#p = PortHamiltonian.new_spectral_element_phs(1,40,0,1)

# pseudo-spectral elements mixed with golo alpha=1 elements:
#p = PortHamiltonian.new_block_spectral_element_phs(Nel,Npol,a,b)

Ns = size(p.J,1)
freq = 1.

function dynamics(t,x, xd)
	xd[:] = p.J * p.GradHam(x) + p.B[:,2] * sin(freq*t) - p.B[:,1] * output(t,x)
end
function output(t,x)
	return (p.B.' * p.GradHam(x) + p.D[:,2] * sin(freq*t))[1]
end

using Sundials
to = collect(linspace(0,10/freq,10000))
x0 = zeros(Ns);
res = Sundials.cvode(dynamics, x0, to)

y = zeros(length(to))
for i = 1:length(y)
	y[i] = output(to[i], res[i,:][:])
end

using PyPlot

figure()
#plot(to+1, sin(freq*to),to, y , to, (blkdiag(diagm(1./wi),diagm(1./wi)) * p.Q*res')[end/2+5,:]')
plot(to+1, -sin(freq*to),to, y)
figure()
plot(to, y+(sin(freq*(to-1)).*(to.>=1)  ))

max_stars_repo_name: flavioluiz/PortHamiltonian

max_issues_repo_name: flavioluiz/PortHamiltonian

max_forks_repo_name: flavioluiz/PortHamiltonian


text:
module PlotlyJSTest
using TestSetExtensions

using Base.Test

using PlotlyJS
const M = PlotlyJS

try
    @testset ExtendedTestSet "PlotlyJS Tests" begin
        @includetests ARGS
    end
catch
    exit(-1)
end

end


text:
function get_num_dens_spectrum(CR_N::Vector{T}, CR_S::Vector{T},
    CR_Cut::T, ρ::T;
    par::CRMomentumDistributionConfig,
    mode::Integer = 3) where {T}

    # transform the norm dependent on IO mode
    norm = transform_norm(CR_N, mode)

    bin_width = log10(par.pmax / par.pmin) / par.Nbins

    bounds = 10.0 .^ collect(log10(par.pmin):par.bin_width:log10(par.pmax))
    bound_up = Vector{Float64}(undef, par.Nbins)

    for i = 1:par.Nbins
        bound_up[i] = (bounds[i+1] > CR_Cut) ? CR_Cut : bounds[i+1]
    end

    bound_up[end] = (bounds[end] > CR_Cut) ? CR_Cut : bounds[end]

    bin_centers = Vector{Float64}(undef, par.Nbins)
    density = Vector{Float64}(undef, par.Nbins)

    @inbounds for i = 1:par.Nbins
        bin_centers[i] = 10.0^(0.5 * (log10(bounds[i]) + log10(bound_up[i])))
        density[i] = density_integral(bounds[i], bound_up[i], norm[i], CR_S[i], ρ)
    end

    return CR_NumSpectrum(bin_centers, density)
end


text:
using InterpretableModels: zerooneconvert, optimize

@testset "LIMScoringSystem via Lasso" begin

N, P = 1_000, 500
realcoefs = zeros(Int, P)
realcoefs[1] = 4
realcoefs[42] = -2
realcoefs[244] = -8
X = rand(N, P) .< 0.5
yraw = X * realcoefs
y = zerooneconvert.(yraw) # class is 0 or 1
models = optimize(LIMScoringSystem, LassoLIMOptimizer(), X, y)

@test models isa AbstractVector
@test models[1] isa LIMScoringSystem

end


text:
struct LambaEMConstantCache <: StochasticDiffEqConstantCache end
@cache struct LambaEMCache{uType,rateType,rateNoiseType,randType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  du2::rateType
  K::rateType
  tmp::uType
  L::rateNoiseType
  gtmp::rateNoiseType
  dW_cache::randType
end

alg_cache(alg::LambaEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = LambaEMConstantCache()

function alg_cache(alg::LambaEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  du1 = zero(rate_prototype); du2 = zero(rate_prototype)
  K = zero(rate_prototype); tmp = zero(u);
  L = zero(noise_rate_prototype)
  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    dW_cache = nothing
  else
    dW_cache = zero(ΔW)
  end
  LambaEMCache(u,uprev,du1,du2,K,tmp,L,gtmp,dW_cache)
end

struct LambaEulerHeunConstantCache <: StochasticDiffEqConstantCache end
@cache struct LambaEulerHeunCache{uType,rateType,rateNoiseType,randType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  du1::rateType
  du2::rateType
  K::rateType
  tmp::uType
  L::rateNoiseType
  gtmp::rateNoiseType
  dW_cache::randType
end

alg_cache(alg::LambaEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = LambaEulerHeunConstantCache()

function alg_cache(alg::LambaEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
  du1 = zero(rate_prototype); du2 = zero(rate_prototype)
  K = zero(rate_prototype); tmp = zero(u);
  L = zero(noise_rate_prototype)
  gtmp = zero(noise_rate_prototype)
  if is_diagonal_noise(prob)
    dW_cache = nothing
  else
    dW_cache = zero(ΔW)
  end
  LambaEulerHeunCache(u,uprev,du1,du2,K,tmp,L,gtmp,dW_cache)
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


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

@doc raw"""current status of a horizontal pod autoscaler

    IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus(;
        currentCPUUtilizationPercentage=nothing,
        currentReplicas=nothing,
        desiredReplicas=nothing,
        lastScaleTime=nothing,
        observedGeneration=nothing,
    )

    - currentCPUUtilizationPercentage::Int32 : current average CPU utilization over all pods, represented as a percentage of requested CPU, e.g. 70 means that an average pod is using now 70% of its requested CPU.
    - currentReplicas::Int32 : current number of replicas of pods managed by this autoscaler.
    - desiredReplicas::Int32 : desired number of replicas of pods managed by this autoscaler.
    - lastScaleTime::IoK8sApimachineryPkgApisMetaV1Time : last time the HorizontalPodAutoscaler scaled the number of pods; used by the autoscaler to control how often the number of pods is changed.
    - observedGeneration::Int64 : most recent generation observed by this autoscaler.
"""
mutable struct IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus <: SwaggerModel
    currentCPUUtilizationPercentage::Any # spec type: Union{ Nothing, Int32 } # spec name: currentCPUUtilizationPercentage
    currentReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: currentReplicas
    desiredReplicas::Any # spec type: Union{ Nothing, Int32 } # spec name: desiredReplicas
    lastScaleTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastScaleTime
    observedGeneration::Any # spec type: Union{ Nothing, Int64 } # spec name: observedGeneration

    function IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus(;currentCPUUtilizationPercentage=nothing, currentReplicas=nothing, desiredReplicas=nothing, lastScaleTime=nothing, observedGeneration=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("currentCPUUtilizationPercentage"), currentCPUUtilizationPercentage)
        setfield!(o, Symbol("currentCPUUtilizationPercentage"), currentCPUUtilizationPercentage)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("currentReplicas"), currentReplicas)
        setfield!(o, Symbol("currentReplicas"), currentReplicas)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("desiredReplicas"), desiredReplicas)
        setfield!(o, Symbol("desiredReplicas"), desiredReplicas)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("lastScaleTime"), lastScaleTime)
        setfield!(o, Symbol("lastScaleTime"), lastScaleTime)
        validate_property(IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus, Symbol("observedGeneration"), observedGeneration)
        setfield!(o, Symbol("observedGeneration"), observedGeneration)
        o
    end
end # type IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus

const _property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus = Dict{Symbol,Symbol}(Symbol("currentCPUUtilizationPercentage")=>Symbol("currentCPUUtilizationPercentage"), Symbol("currentReplicas")=>Symbol("currentReplicas"), Symbol("desiredReplicas")=>Symbol("desiredReplicas"), Symbol("lastScaleTime")=>Symbol("lastScaleTime"), Symbol("observedGeneration")=>Symbol("observedGeneration"))
const _property_types_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus = Dict{Symbol,String}(Symbol("currentCPUUtilizationPercentage")=>"Int32", Symbol("currentReplicas")=>"Int32", Symbol("desiredReplicas")=>"Int32", Symbol("lastScaleTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("observedGeneration")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV1HorizontalPodAutoscalerStatus)
    (getproperty(o, Symbol("currentReplicas")) === nothing) && (return false)
    (getproperty(o, Symbol("desiredReplicas")) === nothing) && (return false)
    true
end

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


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
using DataFrames
using Turing
using Plots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %%
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
sort!(d, :weight)  # to make the plots look nice

@model function height_log(weight, height)
    # β: A person with 100kg will be maybe 180cm. log10(100) = 2. Missing factor: 90
    β ~ Normal(90, 20)
    σ ~ Exponential(5)
    μ = f_log(weight, β)
    height ~ MvNormal(μ, σ)
end

f_log(weight, β) = β .* log10.(weight)

q = quap(height_log(d.weight, d.height))
post = DataFrame(rand(q.distr, 1000)', ["β", "σ"])

mu = f_log(d.weight, post.β')
mu_m, mu_lower, mu_upper = meanlowerupper(mu)

scatter(d.weight, d.height, alpha = 0.3, legend = false)
plot!(d.weight, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m))

# %%
# Some predictions from the inferred values.
estimmean, estimlower, estimupper = estimparam(post)
estimmean[1] * log10(2)
"""
I've left out an intercept again because it doesn't make any sense. If you want to have the
best fit you should definitely put it in though.
β: A person double your weight is expected to be on average 28 cm taller.
"""

# %%
sim = rand.(Normal.(mu, post.σ'))
sim_m, sim_lower, sim_upper = meanlowerupper(sim)

scatter(d.weight, d.height, alpha = 0.3, legend = false)
plot!(d.weight, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m))
plot!(d.weight, sim_lower, fillrange = sim_upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %%
"""
Uhhh, the model looks pretty screwed towards the beginning there...
"""


text:
module API

using ..VsRecorderBase: VsConfig,
    parse_frame!, parse_frame, parse_image,
    parse_stream, parse_video, parse_camera

export parse_frame!, parse_frame, parse_image,
    parse_stream, parse_video, parse_camera

end


text:
using Random
using Images

include("../../configs.jl")

mutable struct Image_Reader
    augment::Bool
    pre_scales::Array{Float64}
    function Image_Reader(augment::Bool=false)
        crop_ratios = [0.6, 0.8, 1.0]
        return new(augment, crop_ratios)
    end
end

# r is the reader object
function read_img(dir; len=nothing, boxes=nothing, r=nothing)

    img = convert(Array{Int64}, rawview(channelview(load(dir))))
    new_boxes = boxes === nothing ? nothing : deepcopy(boxes)
    c, h, w = size(img)

    if len === nothing
        len = max(w, h) # makes image square without changing its size, for evaluation
        div = Int(ceil(len / 64))
        len = 64 * div
    end

    if r !== nothing && r.augment
        img = convert(Array{Float32}, img) ./ 255
        roi = _find_roi(r, img)
        # horizontal flip
        if rand() < 0.5 img, new_boxes = flip_horizontal(img, new_boxes) end
        # color distortion
        img = distort_color(img)
        # cropping the picture with a scale
        img, new_boxes = crop_image(img, new_boxes, roi)
        # resizing to <len>
        img, new_boxes = resize_square_img(img, new_boxes, len, (roi[3] - roi[1]))
        img .*= 255
    else
        img = convert(Array{Float32}, img) ./ 255
        img = squaritize_img(img, len)
        # img, new_boxes = resize_square_img(img, new_boxes, len, maxlen)
        img .*= 255
    end

    if new_boxes !== nothing
        new_boxes[new_boxes .< 0] .= -1
        bbox_indices = getindex.(findall(new_boxes[1,:] .>= 0)) # images having bboxes
        new_boxes = new_boxes[:,bbox_indices]
    end

    img .-= avg_img
    return img, new_boxes
end

function _find_roi(r::Image_Reader, img)
    c, h, w = size(img)
    maxlen = max(h, w); minlen = min(h, w)
    scale = r.pre_scales[rand(1:length(r.pre_scales))]
    new_short = Int(ceil(minlen * scale))
    x1 = 1; y1 = 1
    if scale < 1
        x1 = rand(1:max(1, w-new_short))
        y1 = rand(1:max(1, h-new_short))
    end
    x2 = (new_short-1) + x1
    y2 = (new_short-1) + y1
    #region of interest is determined
    return [x1 y1 x2 y2]
end

function crop_image(img, new_boxes, roi)
    img = img[:, roi[2]:roi[4], roi[1]:roi[3]]
    box_size = size(new_boxes)

    for person in 1:box_size[2]
        # setting landmarks outside of the region to -1
        for lm in 5:2:box_size[1]-2
            # if not in the cropped area
            if (new_boxes[lm, person] > roi[3]) ||
                (new_boxes[lm, person] < roi[1]) ||
                (new_boxes[lm + 1, person] > roi[4]) ||
                (new_boxes[lm + 1, person] < roi[2])
                new_boxes[lm, person] = -1
                new_boxes[lm + 1, person] = -1
            # if in the cropped area
            else
                new_boxes[lm, person] -= roi[1]
                new_boxes[lm + 1, person] -= roi[2]
            end
        end
        if !all(>=(0), new_boxes[5:14,person])
            new_boxes[15, person] = -1.0
        end
        # setting bounding box coordinates

        # completely off the cropped area
        if (new_boxes[3,person] <= roi[1]) || (new_boxes[1,person] >= roi[3])
            new_boxes[1:4, person] .= -1
        elseif (new_boxes[4,person] <= roi[2]) || (new_boxes[2,person] >= roi[4])
            new_boxes[1:4, person] .= -1
        else
            # partial fits
            if new_boxes[3,person] > roi[3]
                new_boxes[3,person] = roi[3]
            end
            if new_boxes[4,person] > roi[4]
                new_boxes[4,person] = roi[4]
            end

            if 0 <= new_boxes[1,person] && new_boxes[1,person] < roi[1]
                new_boxes[1,person] = 0
            else
                new_boxes[1,person] -= roi[1]
            end
            new_boxes[3,person] -= roi[1]

            if 0 <= new_boxes[2,person] && new_boxes[2,person] < roi[2]
                new_boxes[2,person] = 0
            else
                new_boxes[2,person] -= roi[2]
            end
            new_boxes[4,person] -= roi[2]
        end
    end

    return img, new_boxes
end

function resize_square_img(img, new_boxes, new_len, old_len)
    img = imresize(img, (3, new_len, new_len))
    ratio = new_len / old_len
    if new_boxes !== nothing
        new_boxes = ratio .* new_boxes
    end
    return img, new_boxes
end

function squaritize_img(img, len)
    c, h, w = size(img)
    full_img = zeros(3, len, len)
    full_img[:,1:h,1:w] .= img
    return full_img
end

function flip_horizontal(img, new_boxes)
    img = reverse(img, dims=3)
    new_boxes[1:2:13,:] .= size(img, 3) .- new_boxes[1:2:13,:]
    temp = new_boxes[3,:]
    new_boxes[3,:] .= new_boxes[1,:]
    new_boxes[1,:] .= temp
    new_boxes[new_boxes .> size(img, 3)] .= -1
    return img, new_boxes
end

function distort_color(img)
    # brightness distortion
    if rand() < 0.5 img = _convert(img, beta=rand(-0.125:0.125)) end
    # contrast distortion
    if rand() < 0.5 img = _convert(img, alpha=rand(0.5:1.5)) end

    img = channelview(colorview(HSV, float.(img)))
    # saturation distortion
    if rand() < 0.5 img[2,:,:] = _convert(img[2,:,:], alpha=rand(0.5:1.5)) end
    # hue distortion
    if rand() < 0.5 img[1,:,:] = _convert(img[1,:,:], beta=rand(-(18/256):(18/256))) end

    img = channelview(colorview(RGB, colorview(HSV, img)))
    return img
end

function _convert(image; alpha=1, beta=0)
    image = image .* alpha .+ beta
    image[image .< 0] .= 0
    image[image .> 1] .= 1
    return image
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


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

mutable struct IoK8sApiCoreV1StorageOSVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalObjectReference } # spec name: secretRef
    volumeName::Any # spec type: Union{ Nothing, String } # spec name: volumeName
    volumeNamespace::Any # spec type: Union{ Nothing, String } # spec name: volumeNamespace

    function IoK8sApiCoreV1StorageOSVolumeSource(;fsType=nothing, readOnly=nothing, secretRef=nothing, volumeName=nothing, volumeNamespace=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1StorageOSVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1StorageOSVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1StorageOSVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1StorageOSVolumeSource, Symbol("volumeName"), volumeName)
        setfield!(o, Symbol("volumeName"), volumeName)
        validate_property(IoK8sApiCoreV1StorageOSVolumeSource, Symbol("volumeNamespace"), volumeNamespace)
        setfield!(o, Symbol("volumeNamespace"), volumeNamespace)
        o
    end
end # type IoK8sApiCoreV1StorageOSVolumeSource

const _property_map_IoK8sApiCoreV1StorageOSVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("volumeName")=>Symbol("volumeName"), Symbol("volumeNamespace")=>Symbol("volumeNamespace"))
const _property_types_IoK8sApiCoreV1StorageOSVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1LocalObjectReference", Symbol("volumeName")=>"String", Symbol("volumeNamespace")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1StorageOSVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1StorageOSVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1StorageOSVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1StorageOSVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1StorageOSVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1StorageOSVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1StorageOSVolumeSource)
    true
end

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