text:
const MapableTypes = Union{AbstractString, Symbol}
struct NamedCategorical{N<:MapableTypes} <: DiscreteUnivariateDistribution
    cat::Categorical
    map::CategoricalDiscretizer{N, Int}
end
function NamedCategorical(items::AbstractVector{N}, probs::Vector{Float64}) where {N<:MapableTypes}
    cat = Categorical(probs./sum(probs))
    map = CategoricalDiscretizer(items)
    NamedCategorical{N}(cat, map)
end

function Base.show(io::IO, d::NamedCategorical)
    println(io, "NamedCategorical with entries:")
    for val in keys(d.map.n2d)
        @printf(io, "\t%8.4f:  %s\n", pdf(d, val), string(val))
    end
end

Distributions.ncategories(d::NamedCategorical) = Distributions.ncategories(d.cat)
Distributions.probs(d::NamedCategorical) = Distributions.probs(d.cat)
Distributions.params(d::NamedCategorical) = Distributions.params(d.cat)

Distributions.pdf(d::NamedCategorical{N}, x::N) where {N<:MapableTypes} = Distributions.pdf(d.cat, encode(d.map, x))
Distributions.logpdf(d::NamedCategorical{N}, x::N) where {N<:MapableTypes} = Distributions.logpdf(d.cat, encode(d.map, x))
Base.rand(d::NamedCategorical) = rand(sampler(d))

struct MappedAliasTable <: Sampleable{Univariate,Discrete}
   alias::Distributions.AliasTable
   map::CategoricalDiscretizer
end
Distributions.ncategories(s::MappedAliasTable) = Distributions.ncategories(s.alias)
Random.rand(s::MappedAliasTable) = decode(s.map, rand(s.alias))
Base.show(io::IO, s::MappedAliasTable) = @printf(io, "MappedAliasTable with %d entries", ncategories(s))

Distributions.sampler(d::NamedCategorical) = MappedAliasTable(Distributions.sampler(d.cat), d.map)

ext: jl


text:
"""
Case 12:
This case study a three bus system with 2 machines (Classic Model - Single Shaft: 2 State model) without loads.
The machine at bus 1 is used as a reference machine, while machine at bus 2 has a simplified droop governor (TGTypeII).
The perturbation trips one of the four (out of 5) circuits of line between buses 1 and 2, multiplying by 4 its impedance.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test12.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#Time span
tspan = (0.0, 5.0)

#Define Fault: Change of YBus
Ybus_change = ThreePhaseFault(
    1.0, #change at t = 1.0
    Ybus_fault,
) #New YBus

path = (joinpath(pwd(), "test-12"))
!isdir(path) && mkdir(path)
try
    sim = Simulation(
        path,
        threebus_sys, #system
        tspan, #time span
        Ybus_change, #Type of Fault
    )

    small_sig = small_signal_analysis(sim)

    #Run simulation
    run_simulation!(
        sim, #simulation structure
        IDA(),#Sundials DAE Solver
        dtmax = 0.02, #keywords arguments
    )

    series = get_state_series(sim, ("generator-102-1", :ω))

    diff = [0.0]
    res = get_init_values_for_comparison(sim)
    for (k, v) in test12_x0_init
        diff[1] += LinearAlgebra.norm(res[k] - v)
    end
    @test (diff[1] < 1e-3)
    @test sim.solution.retcode == :Success
    @test small_sig.stable
finally
    @info("removing test files")
    rm(path, force = true, recursive = true)
end


text:
#Define abstract planet type that functions as a wrapper of the CelestialBody struct of AstroDynBase
abstract type abstractCelestialBody end
abstract type abstractBarycenter <: abstractCelestialBody end
abstract type abstractPlanet <: abstractCelestialBody end
abstract type abstractNaturalSatellite <: abstractCelestialBody end
abstract type abstractStar <: abstractCelestialBody end

struct Planet <: abstractPlanet
    Name::Type{<:abstractPlanet}
    Environment::T where {T <: AbstractPlanetEnvironment}
end

struct NaturalSatellite <: abstractNaturalSatellite
    Name::Type{<:abstractNaturalSatellite}
    Environment::T where {T <: AbstractPlanetEnvironment}
end

struct Star <: abstractStar
    Name::Type{<:abstractStar}
end

#Unless otherwise defined, define initial values for each of the variables
j2(::Type{C}) where {C <: abstractCelestialBody} = 0.
j3(::Type{C}) where {C <: abstractCelestialBody} = 0.
j4(::Type{C}) where {C <: abstractCelestialBody} = 0.
frame(::Type{C}) where {C <: abstractCelestialBody} = "J2000"

#Wrap the functions
naif_id(Planet::T) where {T <: abstractCelestialBody} = naif_id(Planet.Name)
mu(Planet::T) where {T <: abstractCelestialBody} = mu(Planet.Name)
mean_radius(Planet::T) where {T <: abstractCelestialBody} = mean_radius(Planet.Name)
equatorial_radius(Planet::T) where {T <: abstractCelestialBody} = equatorial_radius(Planet.Name)
polar_radius(Planet::T) where {T <: abstractCelestialBody} = polar_radius(Planet.Name)
parent(Planet::T) where {T <: abstractCelestialBody} = parent(Planet.Name)
right_ascension(Planet::T) where {T <: abstractCelestialBody} = right_ascension(Planet.Name)
right_ascension_rate(Planet::T) where {T <: abstractCelestialBody} = right_ascension_rate(Planet.Name)
declination(Planet::T) where {T <: abstractCelestialBody} = declination(Planet.Name)
declination_rate(Planet::T) where {T <: abstractCelestialBody} = declination_rate(Planet.Name)
rotation_angle(Planet::T) where {T <: abstractCelestialBody} = rotation_angle(Planet.Name)
rotation_rate(Planet::T) where {T <: abstractCelestialBody} = rotation_rate(Planet.Name)
euler_angles(Planet::T) where {T <: abstractCelestialBody} = euler_angles(Planet.Name)
euler_derivatives(Planet::T) where {T <: abstractCelestialBody} = euler_derivatives(Planet.Name)
j2(Planet::T) where {T <: abstractCelestialBody} = j2(Planet.Name)
j3(Planet::T) where {T <: abstractCelestialBody} = j3(Planet.Name)
j4(Planet::T) where {T <: abstractCelestialBody} = j4(Planet.Name)
frame(Planet::T) where {T <: abstractCelestialBody} = frame(Planet.Name)

#Compute the flattening of a CelestialBody
flattening(Body::T) where {T <: abstractCelestialBody} = 1 - (polar_radius(Body) / equatorial_radius(Body))

#Helper state-transformation matrix function
function BodyFixed2InertialState(Body::T, t::S) where {T <: abstractCelestialBody, S <: Real}
    sxform(frame(Body),"J2000",Float64(t))
end

function Inertial2BodyFixedState(Body::T, t::S) where {T <: abstractCelestialBody, S <: Real}
    sxform("J2000",frame(Body),Float64(t))
end

function BodyFixed2Inertial3(Body::T, t::S) where {T <: abstractCelestialBody, S <: Real}
    pxform(frame(Body),"J2000",Float64(t))
end

function Inertial2BodyFixed3(Body::T, t::S) where {T <: abstractCelestialBody, S <: Real}
    pxform("J2000",frame(Body),Float64(t))
end

function CelestialBodyState(FromBody::T, ToBody::P, t::S) where {T,P <: abstractCelestialBody, S <: Real}
    spkezr(naif_id(FromBody), Float64(t), "J2000", naif_id(ToBody))[1]
end

function CelestialBodyPosition(FromBody::T, ToBody::P, t::S) where {T,P <: abstractCelestialBody, S <: Real}
    spkpos(naif_id(FromBody), Float64(t), "J2000", naif_id(ToBody))[1]
end


text:
@doc raw"""
    MultinomialMatrices{n,m} <: AbstractPowerManifold{ℝ}

The multinomial manifold consists of `m` column vectors, where each column is of length
`n` and unit norm, i.e.

````math
\mathcal{MN}(n,m) \coloneqq \bigl\{ p ∈ ℝ^{n×m}\ \big|\ p_{i,j} > 0 \text{ for all } i=1,…,n, j=1,…,m \text{ and } p^{\mathrm{T}}\mathbb{1}_m = \mathbb{1}_n\bigr\},
````
where $\mathbb{1}_k$ is the vector of length $k$ containing ones.

This yields exactly the same metric as
considering the product metric of the probablity vectors, i.e. [`PowerManifold`](@ref) of the
$(n-1)$-dimensional [`ProbabilitySimplex`](@ref).

The [`ProbabilitySimplex`](@ref) is stored internally within `M.manifold`, such that all functions of
[`AbstractPowerManifold`](@ref) can be used directly.

# Constructor

    MultinomialMatrices(n, m)

Generate the manifold of matrices $\mathbb R^{n×m}$ such that the $m$ columns are
discrete probability distributions, i.e. sum up to one.
"""
struct MultinomialMatrices{N,M,S} <:
       AbstractPowerManifold{ℝ,ProbabilitySimplex{S},ArrayPowerRepresentation} where {N,M}
    manifold::ProbabilitySimplex{S}
end

function MultinomialMatrices(n::Int, m::Int)
    return MultinomialMatrices{n,m,n - 1}(ProbabilitySimplex(n - 1))
end

function Base.:^(M::ProbabilitySimplex{N}, m::Int) where {N}
    return MultinomialMatrices{manifold_dimension(M) + 1,m,N}(M)
end

@doc raw"""
    check_point(M::MultinomialMatrices, p)

Checks whether `p` is a valid point on the [`MultinomialMatrices`](@ref)`(m,n)` `M`, i.e. is a matrix
of `m` discrete probability distributions as columns from $\mathbb R^{n}$, i.e. each column is a point from
[`ProbabilitySimplex`](@ref)`(n-1)`.
"""
check_point(::MultinomialMatrices, ::Any)
function check_point(M::MultinomialMatrices{n,m}, p; kwargs...) where {n,m}
    if size(p) != (n, m)
        return DomainError(
            length(p),
            "The matrix in `p` ($(size(p))) does not match the dimensions of $(M).",
        )
    end
    return check_point(PowerManifold(M.manifold, m), p; kwargs...)
end
@doc raw"""
    check_vector(M::MultinomialMatrices p, X; kwargs...)

Checks whether `X` is a valid tangent vector to `p` on the [`MultinomialMatrices`](@ref) `M`.
This means, that `p` is valid, that `X` is of correct dimension and columnswise
a tangent vector to the columns of `p` on the [`ProbabilitySimplex`](@ref).
"""
function check_vector(M::MultinomialMatrices{n,m}, p, X; kwargs...) where {n,m}
    if size(X) != (n, m)
        return DomainError(
            length(X),
            "The matrix `X` ($(size(X))) does not match the required dimension ($(representation_size(M))) for $(M).",
        )
    end
    return check_vector(PowerManifold(M.manifold, m), p, X; kwargs...)
end

get_iterator(::MultinomialMatrices{n,m}) where {n,m} = Base.OneTo(m)

@generated manifold_dimension(::MultinomialMatrices{n,m}) where {n,m} = (n - 1) * m

@generated representation_size(::MultinomialMatrices{n,m}) where {n,m} = (n, m)

function Base.show(io::IO, ::MultinomialMatrices{n,m}) where {n,m}
    return print(io, "MultinomialMatrices($(n),$(m))")
end


text:
#=
The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number 600851475143 ?
=#

using Primes
@time println(maximum(keys(factor(600851475143))))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
"""
    bipartite_map(g)

For a bipartite graph `g`, return a vector `c` of size ``|V|`` containing
the assignment of each vertex to one of the two sets (``c_i == 1`` or c_i == 2``).
If `g` is not bipartite, return an empty vector.

### Implementation Notes
Note that an empty vector does not necessarily indicate non-bipartiteness.
An empty graph will return an empty vector but is bipartite.
"""
function bipartite_map(g::AbstractGraph{T}) where T
    nvg = nv(g)
    if !is_directed(g)
        ccs = filter(x -> length(x) >= 2, connected_components(g))
    else
        ccs = filter(x -> length(x) >= 2, weakly_connected_components(g))
    end
    seen = zeros(Bool, nvg)
    colors = zeros(Bool, nvg)
    for cc in ccs
        Q = Vector{T}()
        s = cc[1]
        push!(Q, s)
        bipartitemap = zeros(UInt8, nvg)
        while !isempty(Q)
            u = popfirst!(Q)
            for v in outneighbors(g, u)
                if !seen[v]
                    colors[v] = !colors[u]
                    push!(Q, v)
                    seen[v] = true
                elseif colors[v] == colors[u]
                    return Vector{UInt8}()
                end
            end
        end
    end
    return UInt8.(colors).+(one(UInt8))
end

"""
    is_bipartite(g)

Return `true` if graph `g` is [bipartite](https://en.wikipedia.org/wiki/Bipartite_graph).
"""
is_bipartite(g::AbstractGraph) = length(bipartite_map(g)) == nv(g)


text:
include("../src/GISSMOReader.jl")
import .GISSMOReader

import Unicode
import JSON
import JSON3

output_folder = "/home/roy/Documents/repo/NMRData/input/compound_mapping"
mkpath(output_folder)
save_path = joinpath(output_folder, "GISSMO_names.json")

JSON_folder = "/home/roy/Documents/repo/GISSMOReader.jl/GISSMO_info"
JSON_filename = "GISSMO_entries_metadata_Jan_2022.json"
JSON_load_path = joinpath(JSON_folder, JSON_filename)

function createGISSMOcompoundnamesdict(JSON_load_path)
    db_dict = JSON.parse(open(JSON_load_path))

    dict_out = Dict()
    for (key, entry_dict) in db_dict

        compound_name = entry_dict["molecule_name"]

        U = split(key, "/")
        entry_name = U[end-1]
        simulation_name = U[end]

        key_out = "$(compound_name) $(entry_name) $(simulation_name)"
        file_name = "$(entry_name)_$(simulation_name).json"
        notes = key

        dict_out[key_out] = Dict( "file name" => file_name, "notes" => notes)
    end

    return dict_out
end

dict_out = createGISSMOcompoundnamesdict(JSON_load_path)

stringdata = JSON.json(dict_out)

open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end

#### common compounds.

key_names = ["DSS";
"L-Tyrosine";
"L-Leucine";
"Glycine";
"L-Isoleucine";
"L-Valine";
"L-Tryptophan";
"L-Threonine";
"L-Serine";
"L-Phenylalanine";
"L-Methionine";
"L-Glutamine";
"L-(+) Lactic acid";
"D-(+)-Glucose";
"Caffeine";
"L-Histidine";
"L-Proline";
"Choline";
"L-Asparagine";
"L-Aspartic acid";
"L-Lysine";
"L-Alanine";
"L-Arginine";
"Folate";
"L-Cysteine";
"Putrescine";
"Succinic acid";
"Purine";
"L-Glutamic acid";
"Ethanol";
"Epinephrine";
"Norepinephrine";
"3,4-Dihydroxy-L-phenylalanine";
"Dopamine";
"5-Hydroxy-L-tryptophan";
"Serotonin";
"L-Kynurenine";
"Tryptamine";
"Creatine";
"Creatinine";
"Myo-Inositol";
"Glycerol";
"Gamma-Aminobutyric acid";
"3-Hydroxybutyrate";
"HEPES";
"Acetylcholine";
"DL-Alanine";
"Methanol";
"Triethanolamine";
"DL-Serine";
"Diethanolamine";
"N-Alpha-Acetyl-L-Lysine";
"Phenylacetylglycine";
"2-Phenylethanol"]

search_names = ["DSS";
"L-Tyrosine";
"L-Leucine";
"Glycine";
"L-Isoleucine";
"L-Valine";
"L-Tryptophan";
"L-Threonine";
"L-Serine";
"L-Phenylalanine";
"L-Methionine";
"L-Glutamine";
"Lactic acid";
"Glucose";
"Caffeine";
"L-Histidine";
"L-Proline";
"Choline";
"L-Asparagine";
"L-Aspartic acid";
"L-Lysine";
"L-Alanine";
"L-Arginine";
"Folate";
"L-Cysteine";
"Putrescine";
"Succinic acid";
"Purine";
"L-Glutamic acid";
"Ethanol";
"Epinephrine";
"Norepinephrine";
"Dihydroxy-L-phenylalanine";
"Dopamine";
"5-Hydroxy-L-tryptophan";
"Serotonin";
"L-Kynurenine";
"Tryptamine";
"Creatine";
"Creatinine";
"Myo-Inositol";
"Glycerol";
"Gamma-Aminobutyric acid";
"3-Hydroxybutyrate";
"HEPES";
"Acetylcholine";
"DL-Alanine";
"Methanol";
"Triethanolamine";
"DL-Serine";
"Diethanolamine";
"Acetyl-L-lysine";
"Phenylacetylglycine";
"2-Phenylethanol"]

function createGISSMOcompoundnamesdict(JSON_load_path, search_names::Vector{String})
    db_dict = JSON.parse(open(JSON_load_path))

    dict_out = Dict()
    for (key, entry_dict) in db_dict

        compound_name = entry_dict["molecule_name"]

        U = split(key, "/")
        entry_name = U[end-1]
        simulation_name = U[end]

        key_out = "$(compound_name) $(entry_name) $(simulation_name)"
        file_name = "$(entry_name)_$(simulation_name).json"
        notes = key

        compound_name0 = Unicode.normalize(compound_name, casefold = true)

        if any(occursin(Unicode.normalize(search_names[i], casefold = true), compound_name0) for i = 1:length(search_names))

            dict_out[key_out] = Dict( "file name" => file_name, "notes" => notes)
        end
    end

    return dict_out
end

dict_out_subset = createGISSMOcompoundnamesdict(JSON_load_path, search_names)

stringdata = JSON.json(dict_out_subset)

save_path = joinpath(output_folder, "subset_names.json")
open(save_path, "w") do f
    JSON3.pretty(f, stringdata)
    println(f)
end


text:
using Compat

if VERSION < v"0.3.3"
  function blasfunc(name::Symbol)
    string(name)
  end
elseif VERSION < v"0.5.0-dev+1915"
  function blasfunc(name::Symbol)
    Base.blasfunc(name)
  end
else
  function blasfunc(name::Symbol)
    str_name = string(name)
    fnc_symb = eval(:(Base.BLAS.@blasfunc $str_name))
    fnc_name = string(fnc_symb)
    return fnc_name
  end
end


text:
using HPAT
using HPAT.API.LinearRegression
#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DataTablePass.set_debug_level(3)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

#using ParallelAccelerator
#ParallelAccelerator.DomainIR.set_debug_level(3)

@acc hpat function q05(category, education, gender, file_name)
    web_clickstreams = DataSource(DataTable{:wcs_item_sk=Int64, :wcs_user_sk=Int64}, HDF5, file_name)
    item = DataSource(DataTable{:i_item_sk=Int64,:i_category_id=Int64,:i_category=Int64}, HDF5, file_name)
    customer = DataSource(DataTable{:c_customer_sk=Int64,:c_current_cdemo_sk=Int64}, HDF5, file_name)
    customer_demographics = DataSource(DataTable{:cd_demo_sk=Int64,:cd_gender=Int64,:cd_education_status=Int64}, HDF5, file_name)
    # TODO: following two filters are added by spark as an optimization;
    # I am adding manually. Need to add optimization in datatable pass
    web_clickstreams = web_clickstreams[:wcs_user_sk<typemax(Int32)]
    customer = customer[:c_current_cdemo_sk<typemax(Int32)]
    # typemax(Int32) used for NULL; typemin(Int32) was not working
    web_clickstreams = web_clickstreams[:wcs_item_sk<typemax(Int32)]

    user_items = join(web_clickstreams, item, :wcs_item_sk==:i_item_sk, :user_items_sk)

    user_clicks_in_cat = aggregate(user_items, :wcs_user_sk, :clicks_in_category = sum(:i_category==category),
                                                         :clicks_in_1 = sum(:i_category_id==1),
                                                         :clicks_in_2 = sum(:i_category_id==2),
                                                         :clicks_in_3 = sum(:i_category_id==3),
                                                         :clicks_in_4 = sum(:i_category_id==4),
                                                         :clicks_in_5 = sum(:i_category_id==5),
                                                         :clicks_in_6 = sum(:i_category_id==6),
                                                         :clicks_in_7 = sum(:i_category_id==7))
    customer_clicks = join(user_clicks_in_cat, customer, :wcs_user_sk==:c_customer_sk, :customer_clicks_sk)
    customer_demo_clicks = join(customer_clicks, customer_demographics, :c_current_cdemo_sk==:cd_demo_sk , :customer_demo_clicks_sk)
    # To make int array from bool
    college_education = (customer_demo_clicks[:cd_education_status].==education)
    male = (customer_demo_clicks[:cd_gender].==gender)
    responses = 1.0 .*(customer_demo_clicks[:clicks_in_category])
    data = transpose(typed_hcat(Float64, college_education,
            male,
            customer_demo_clicks[:clicks_in_1],
            customer_demo_clicks[:clicks_in_2],
            customer_demo_clicks[:clicks_in_3],
            customer_demo_clicks[:clicks_in_4],
            customer_demo_clicks[:clicks_in_5],
            customer_demo_clicks[:clicks_in_6],
            customer_demo_clicks[:clicks_in_7]))
    #pointsF = convert(Matrix{Float64},data)
    model = LinearRegression(data, responses)
    #return customer_demo_clicks[:customer_demo_clicks_sk], customer_demo_clicks[:cd_gender]
    return model
end

println(q05(3, 8303423, 1, "test_q05.hdf5"))


text:
module VSCode

import ..JET:
    @jetconfigurable,
    with_logger,
    gen_postprocess,
    tofullpath,
    print_signature,
    AbstractAnalyzer,
    JETToplevelResult,
    ToplevelErrorReport,
    JETCallResult,
    InferenceErrorReport,
    get_reports,
    print_report

# common
# ======

isuntitled(path::AbstractString)   = startswith(path, "Untitled")
tovscodepath(path::Symbol)         = tovscodepath(string(path))
tovscodepath(path::AbstractString) = isuntitled(path) ? path : tofullpath(path)

"""
    vscode_source(analyzer::AbstractAnalyzer, source::AbstractString) -> String

Returns "source" of [diagnostics](https://code.visualstudio.com/api/references/vscode-api#Diagnostic),
which represents the identity of diagnostics.
"""
function vscode_source(analyzer::Analyzer, source::AbstractString) where {Analyzer<:AbstractAnalyzer}
    return string(nameof(Analyzer), ": ", source)
end

"""
    vscode_diagnostics_order(analyzer::AbstractAnalyzer) -> Bool

If `true` (default) a diagnostic will be reported at entry site.
Otherwise it's reported at error point.
"""
vscode_diagnostics_order(analyzer::AbstractAnalyzer) = true

# configuration
# =============

"""
Configurations for the VSCode integration.
These configurations are active only when used in [the integrated Julia REPL](https://www.julia-vscode.org/docs/dev/userguide/runningcode/).

---
- `vscode_console_output::Union{Nothing,IO} = stdout` \\
  JET will show analysis result in VSCode's "PROBLEMS" pane and inline annotations.
  If `vscode_console_output::IO` is specified, JET will also print the result into the
  specified output stream in addition to showing the result in the integrated views.
  When `nothing`, the result will be only shown in the integrated views.
---
"""
struct VSCodeConfig end

forward_to_console_output(res::Union{JETToplevelResult,JETCallResult};
                          vscode_console_output::Union{Nothing,IO} = nothing,
                          __jetconfigs...) =
    isa(vscode_console_output, IO) && show(vscode_console_output, res)

# top-level
# =========

Base.showable(::MIME"application/vnd.julia-vscode.diagnostics", ::JETToplevelResult) = true
function Base.show(io::IO, ::MIME"application/vnd.julia-vscode.diagnostics",
                   res::JETToplevelResult)
    forward_to_console_output(res; res.jetconfigs...)
    return vscode_diagnostics(res.analyzer,
                              get_reports(res),
                              res.source,
                              gen_postprocess(res.res.actual2virtual),
                              )
end
function vscode_diagnostics(analyzer::Analyzer,
                            reports::Vector{ToplevelErrorReport},
                            source::AbstractString,
                            postprocess = identity) where {Analyzer<:AbstractAnalyzer}
    return (; source = vscode_source(analyzer, source),
              items = map(reports) do report
                  return (; msg = postprocess(sprint(print_report, report)),
                            path = tovscodepath(report.file),
                            line = report.line,
                            severity = 0, # 0: Error, 1: Warning, 2: Information, 3: Hint
                            )
              end,
              )
end

# inference
# =========

Base.showable(::MIME"application/vnd.julia-vscode.diagnostics", ::JETCallResult) = true
function Base.show(io::IO, ::MIME"application/vnd.julia-vscode.diagnostics",
                   res::JETCallResult)
    forward_to_console_output(res; res.jetconfigs...)
    return vscode_diagnostics(res.analyzer,
                              get_reports(res),
                              res.source,
                              )
end
function vscode_diagnostics(analyzer::Analyzer,
                            reports::Vector{InferenceErrorReport},
                            source::AbstractString,
                            postprocess = identity) where {Analyzer<:AbstractAnalyzer}
    order = vscode_diagnostics_order(analyzer)
    return (; source = vscode_source(analyzer, source),
              items = map(reports) do report
                  showpoint = (order ? first : last)(report.vst)
                  return (; msg = postprocess(report.msg),
                            path = tovscodepath(showpoint.file),
                            line = showpoint.line,
                            severity = 1, # 0: Error, 1: Warning, 2: Information, 3: Hint
                            relatedInformation = map((order ? identity : reverse)(report.vst)) do frame
                                sig = sprint(print_signature, frame.sig, (; annotate_types = true))
                                return (; msg = postprocess(sig),
                                          path = tovscodepath(frame.file),
                                          line = frame.line,
                                          )
                            end,
                            )
              end,
              )
end

end # module VSCode


text:
using Documenter

is_ci_build = get(ENV, "CI", nothing) == "true"

makedocs(
    sitename="metal-stack",
    format = Documenter.HTML(
        prettyurls = is_ci_build,
        assets = ["assets/favicon.ico", "assets/youtube.css"],
        canonical = "https://docs.metal-stack.io/",
        highlights = ["yaml"],
    ),
    authors = "metal-stack authors and contributors.",
    linkcheck = is_ci_build,
    linkcheck_ignore = [
        r"^(?!http)",
    ],
    #strict = true, # TODO: Should be enabled soon, links from repos to CONTRIBUTING have to be updated though...
    clean = true,
    pages = [
        "Introduction" => "index.md",
        "Overview" => Any[
            "Architecture" => "overview/architecture.md",
            "Networking" => "overview/networking.md",
            "Hardware Support" => "overview/hardware.md",
            "Operating Systems" => "overview/os.md",
            "Kubernetes Integration" => "overview/kubernetes.md",
            "Storage" => "overview/storage.md",
            "Comparison" => "overview/comparison.md",
        ],
        "Quickstart" => "quickstart.md",
        "Installation & Administration" => Any[
            "Installation" => "installation/deployment.md",
            "Monitoring" => "installation/monitoring.md",
            "Troubleshoot" => "installation/troubleshoot.md",
        ],
        "User Guides" => Any[
            "mini-lab" => "external/mini-lab/README.md",
            "metalctl" => "external/metalctl/README.md",
            "csi-driver-lvm" => "external/csi-driver-lvm/README.md",
            "firewall-controller" => "external/firewall-controller/README.md",
        ],
        "API Documentation" => "apidocs/apidocs.md",
        "Development" => Any[
            "development/client_libraries.md",
            "development/roadmap.md",
            "Enhancement Proposals" => "development/proposals/index.md",
            "development/contributing.md",
        ],
    ]
)

if is_ci_build
    deploydocs(
        repo = "github.com/metal-stack/docs.git",
        push_preview = true,
    )
end

max_stars_repo_name: metal-stack/docs

max_issues_repo_name: metal-stack/docs

max_forks_repo_name: metal-stack/docs


text:
macro mlput(ml,ex)
    if typeof(ex) == Expr
        if ex.head == :call
            Expr(:block,
                 :(mlputfunction($(esc(ml)),$(string(ex.args[1])),$(esc(length(ex.args)-1)))),
                 [:(@mlput $(esc(ml)) $(esc(u))) for u = ex.args[2:end]]...)
        elseif ex.head == :vcat || ex.head == :cell1d
            Expr(:block,
                 :(mlputfunction($(esc(ml)),"List",$(esc(length(ex.args))))),
                 [:(@mlput $(esc(ml)) $(esc(u))) for u = ex.args[1:end]]...)
        elseif ex.head == :(=)
            Expr(:block,
                 :(mlputfunction($(esc(ml)),"Set",$(esc(length(ex.args))))),
                 [:(@mlput $(esc(ml)) $(esc(u))) for u = ex.args[1:end]]...)
        elseif ex.head == :(=>)
            Expr(:block,
                 :(mlputfunction($(esc(ml)),"Rule",$(esc(length(ex.args))))),
                 [:(@mlput $(esc(ml)) $(esc(u))) for u = ex.args[1:end]]...)
        elseif ex.head == :quote
            :(mlput($(esc(ml)),$(esc(ex))))
        end
    else
        return :(mlput($(esc(ml)),$(esc(ex))))
    end
end


text:
@testset "emcee.jl" begin
    @testset "example" begin
        @testset "untransformed space" begin
            # define model
            function logprob(θ)
                s, m = θ
                s > 0 || return -Inf

                mdist = Normal(0, sqrt(s))
                obsdist = Normal(m, sqrt(s))

                return logpdf(InverseGamma(2, 3), s) + logpdf(mdist, m) +
                    logpdf(obsdist, 1.5) + logpdf(obsdist, 2.0)
            end
            model = DensityModel(logprob)

            # perform stretch move and sample from prior in initial step
            Random.seed!(100)
            sampler = Ensemble(1_000, StretchProposal([InverseGamma(2, 3), Normal(0, 1)]))
            chain = sample(model, sampler, 1_000;
                           param_names = ["s", "m"], chain_type = Chains)

            @test mean(chain["s"]) ≈ 49/24 atol=0.1
            @test mean(chain["m"]) ≈ 7/6 atol=0.1
        end

        @testset "transformed space" begin
            # define model
            function logprob(θ)
                logs, m = θ
                s = exp(logs)
                sqrts = sqrt(s)

                mdist = Normal(0, sqrts)
                obsdist = Normal(m, sqrts)

                return logpdf(InverseGamma(2, 3), s) + logpdf(mdist, m) +
                    logpdf(obsdist, 1.5) + logpdf(obsdist, 2.0) + logs
            end
            model = DensityModel(logprob)

            # perform stretch move and sample from normal distribution in initial step
            Random.seed!(100)
            sampler = Ensemble(1_000, StretchProposal(MvNormal(2, 1)))
            chain = sample(model, sampler, 1_000;
                           param_names = ["logs", "m"], chain_type = Chains)

            @test mean(exp, chain["logs"]) ≈ 49/24 atol=0.1
            @test mean(chain["m"]) ≈ 7/6 atol=0.1
        end
    end
end


text:
export bgreedy, bgreedy_contraction!

import StatsBase.sample
import StatsBase.weights

# *************************************************************************** #
#                         Boltzmann greedy functions
# *************************************************************************** #

"""
    function bgreedy(network::TensorNetworkCircuit, α::Real, τ::Real)

This function returns a contraction plan for a network by sequentially choosing
random pairs of nodes to contract according to a Boltzmann distribution where
expensive contractions are less likely to be chosen.
"""
function bgreedy(network::TensorNetworkCircuit, α::Real, τ::Real)

    tensors = copy(network.nodes)
    contraction_plan = Array{Array{Symbol, 1}, 1}()
    number_unique_tensors = network.counters["node"]
    time_cost = 0; space_cost = 0

    while length(tensors) > 1
        # Randomly choose two tensors A and B to contract. The resulting
        # tensor is denoted C.
        C, A_label, B_label, costs = get_random_pair(tensors, α, τ)
        number_unique_tensors += 1
        C_label = Symbol("node_$number_unique_tensors")

        # Remove A and B from tensors and add C
        # Append the contraction (A, B) to the contraction plan.
        delete!(tensors, A_label)
        delete!(tensors, B_label)
        tensors[C_label] = C
        append!(contraction_plan, [[A_label, B_label]])

        time_cost += costs[1]
        if space_cost < costs[2]
            space_cost = costs[2]
        end
    end

    contraction_plan, time_cost, space_cost
end

"""
    function bgreedy(network::TensorNetworkCircuit, α::Real, τ::Real, N::Int)

This function samples N contraction plans by calling the bgreedy function
N times and returns the best contraction plan found.
"""
function bgreedy(network::TensorNetworkCircuit, α::Real, τ::Real, N::Int)
    best_plan = []
    best_time_cost = Inf

    for n = 1:N
        plan, time_cost, space_cost = bgreedy(network, α, τ)
        if time_cost < best_time_cost
            best_time_cost = time_cost
            best_plan = plan
        end
    end

    best_plan, best_time_cost
end

"""
    function contraction_cost(A::Node, B::Node, α::Real)

This function returns the heuristic cost of contracting tensors A and B given
by Gray and Kourtis in 'arXiv:2002.01935'. The labels and dimensions for the
uncontracted indices and both the time cost (flops) and space cost (memory) of
the contraction are also returned.
"""
function contraction_cost(A::Node, B::Node, α::Real)
    contracted_indices, C_indices = sort_indices(A, B)

    A_open_dims = [A.dims[i] for i = 1:length(A.dims)
                   if A.indices[i] in C_indices]
    B_open_dims = [B.dims[i] for i = 1:length(B.dims)
                   if B.indices[i] in C_indices]
    C_dims = [A_open_dims; B_open_dims]

    size_C = prod(C_dims)
    size_A = prod(A.dims)
    size_B = prod(B.dims)

    # Note: sqrt(|A||B||C|) = |open dims| x |contracted dims| ∝ flops
    costs = [sqrt(size_A*size_B*size_C), size_C]

    # return the Gray and Kourtis heuristic cost for the given value of α.
    return size_C - α*(size_A + size_B), C_indices, C_dims, costs
end

"""
    function get_random_pair(tensors::Dict{Symbol, Node}, α::Real, τ::Real)

This function will return a pair of tensors A and B sampled from the boltzmann
distribution over the given dictionary of tensors. The cost of contracting A and
B plays the role of energy. (Cheaper contractions are more likely to by chosen.)
A Node instance for the tensor C resulting from contracting A and B is also
returned.
"""
function get_random_pair(tensors::OrderedDict{Symbol, Node}, α::Real, τ::Real)
    # Allocate memory
    N = length(tensors); N_pairs = N*(N-1)÷2
    tensors = collect(tensors)
    tensor_pairs = Array{Array{Symbol, 1}, 1}(undef, N_pairs)
    energy = Array{Real, 1}(undef, N_pairs)
    C_indices = Array{Array{Symbol, 1}, 1}(undef, N_pairs)
    C_dims = Array{Array{Int, 1}, 1}(undef, N_pairs)
    costs = Array{Array{Real, 1}, 1}(undef, N_pairs)

    # Loop over all possible pairings of tensors and compute heuristic cost
    k = 1
    for i = 1:N, j = i+1:N
        (A_label, A) = tensors[i]; (B_label, B) = tensors[j]
        tensor_pairs[k] = [A_label, B_label]
        energy[k], C_indices[k], C_dims[k], costs[k] = contraction_cost(A, B, α)
        k += 1
    end

    # Sample a pair according to the Boltzmann distribution and return
    # the result
    weights_vector = exp.(-(energy)./τ)
    k = sample(1:N_pairs, weights(weights_vector))
    A_label, B_label = tensor_pairs[k]
    C = Node(C_indices[k], C_dims[k], :intermediate_tensor)
    return C, A_label, B_label, costs[k]
end

"""
    function bgreedy_contraction!(network::TensorNetworkCircuit,
                                  α::Real, τ::Real, N::Int)

Function to contract a TensorNetworkCircuit with a contraction plan found by
the bgreedy method with the given parameters.
"""
function bgreedy_contraction!(network::TensorNetworkCircuit,
                              α::Real=1, τ::Real=1, N::Int=10,
                              output_shape::Union{String, Array{Int, 1}}="")
    contraction_plan, cost = bgreedy(network, α, τ, N)
    contract_network!(network, contraction_plan, output_shape)
end


text:
# Implementation of the LOBPCG algorithm, seeking optimal performance
# and stability (it should be very, very hard to break, even at tight
# tolerance criteria). The code is not 100% optimized yet, but the
# most time-consuming parts (BLAS3 and matvecs) should be OK. The
# implementation follows the scheme of Hetmaniuk & Lehoucq (see also
# the refinements in Duersch et. al.), with the following
# modifications:

# - Cholesky factorizations are used instead of the eigenvalue
# decomposition of Stathopolous & Wu for the orthogonalization.
# Cholesky orthogonalization is fast but has an unwarranted bad
# reputation: when applied twice to a matrix X with κ(X) <~
# sqrt(1/ε), where ε is the machine epsilon, it will produce a set
# of very orthogonal vectors, just as the eigendecomposition-based
# method. It can fail when κ >~ sqrt(1/ε), but that can be fixed by
# shifting the overlap matrix. This is very reliable while being
# much faster than eigendecompositions.

# - Termination criteria for the orthogonalization are based on cheap
# estimates instead of costly explicit checks.

# - The default tolerances are very tight, yielding a very stable
# algorithm. This can be tweaked (see keyword ortho_tol) to
# compromise on stability for less Cholesky factorizations.

# - Implicit product updates (reuse of matvecs) are performed whenever
# it is safe to do so, ie only with orthogonal transformations. An
# exception is the B matrix reuse, which seems to be OK even with very
# badly conditioned B matrices. The code is easy to modify if this is
# not the case.

# - The locking is performed carefully to ensure minimal impact on the
# other eigenvectors (which is not the case in many - all ? - other
# implementations)

## TODO micro-optimization of buffer reuse
## TODO write a version that doesn't assume that B is well-conditionned, and doesn't reuse B applications at all
## TODO it seems there is a lack of orthogonalization immediately after locking, maybe investigate this to save on some choleskys
## TODO debug orthogonalizations when A=I

# vprintln(args...) = println(args...)  # Uncomment for output
vprintln(args...) = nothing

using LinearAlgebra
using BlockArrays # used for the `mortar` command which makes block matrices

# when X or Y are BlockArrays, this makes the return value be a proper array (not a BlockArray)
function array_mul(X::AbstractArray{T}, Y) where T
    Z = Array{T}(undef, size(X, 1), size(Y, 2))
    mul!(Z, X, Y)
end

# Perform a Rayleigh-Ritz for the N first eigenvectors.
@timing function rayleigh_ritz(X, AX, N)
    F = eigen(Hermitian(array_mul(X', AX)))
    F.vectors[:,1:N], F.values[1:N]
end

# B-orthogonalize X (in place) using only one B apply.
# This uses an unstable method which is only OK if X is already
# orthogonal (not B-orthogonal) and B is relatively well-conditioned
# (which implies that X'BX is relatively well-conditioned, and
# therefore that it is safe to cholesky it and reuse the B apply)
function B_ortho!(X, BX)
    O = Hermitian(X'*BX)
    U = cholesky(O).U
    rdiv!(X, U)
    rdiv!(BX, U)
end

normest(M) = maximum(abs.(diag(M))) + norm(M - Diagonal(diag(M)))
# Orthogonalizes X to tol
# Returns the new X, the number of Cholesky factorizations algorithm, and the
# growth factor by which small perturbations of X can have been
# magnified
@timing function ortho!(X::AbstractArray{T}; tol=2eps(real(T))) where T
    local R

    # # Uncomment for "gold standard"
    # U,S,V = svd(X)
    # return U*V', 1, 1

    growth_factor = 1

    success = false
    nchol = 0
    while true
        O = Hermitian(X'X)
        try
            R = cholesky(O).U
            nchol += 1
            success = true
        catch err
            @assert isa(err, PosDefException)
            vprintln("fail")
            # see https://arxiv.org/pdf/1809.11085.pdf for a nice analysis
            # We are not being very clever here; but this should very rarely happen
            # so it should be OK
            α = 100
            nbad = 0
            while true
                O += α*eps(real(T))*norm(X)^2*I
                α *= 10
                try
                    R = cholesky(O).U
                    nchol += 1
                    break
                catch err
                    @assert isa(err, PosDefException)
                end
                nbad += 1
                if nbad > 10
                    error("Cholesky shifting is failing badly, this should never happen")
                end
            end
            success = false
        end
        invR = inv(R)
        rmul!(X, invR) # we do not use X/R because we use invR next

        # We would like growth_factor *= opnorm(inv(R)) but it's too
        # expensive, so we use an upper bound which is sharp enough to
        # be close to 1 when R is close to I, which we need for the
        # outer loop to terminate
        # ||invR|| = ||D + E|| ≤ ||D|| + ||E|| ≤ ||D|| + ||E||_F,
        # where ||.|| is the 2-norm and ||.||_F the Frobenius

        norminvR = normest(invR)
        growth_factor *= norminvR

        # condR = 1/LAPACK.trcon!('I', 'U', 'N', Array(R))
        condR = normest(R)*norminvR # in practice this seems to be an OK estimate

        vprintln("Ortho(X) success? $success ", eps(real(T))*condR^2, " < $tol")

        # a good a posteriori error is that X'X - I is eps()*κ(R)^2;
        # in practice this seems to be sometimes very overconservative
        success && eps(real(T))*condR^2 < tol && break

        nchol > 10 && error("Ortho(X) is failing badly, this should never happen")
    end

    # @assert norm(X'X - I) < tol

    X, nchol, growth_factor
end

# Randomize the columns of X if the norm is below tol
function drop!(X::AbstractArray{T}, tol=2eps(real(T))) where T
    dropped = Int[]
    for i=1:size(X,2)
        n = norm(@views X[:,i])
        if n <= tol
            X[:,i] = randn(T, size(X,1))
            push!(dropped, i)
        end
    end
    dropped
end

# Find X that is orthogonal, and B-orthogonal to Y, up to a tolerance tol.
@timing "ortho! X vs Y" function ortho!(X::AbstractArray{T}, Y, BY; tol=2eps(real(T))) where T
    # normalize to try to cheaply improve conditioning
    Threads.@threads for i=1:size(X,2)
        n = norm(@views X[:,i])
        @views X[:,i] ./= n
    end

    niter = 1
    ninners = zeros(Int,0)
    while true
        BYX = BY'X
        # XXX the one(T) instead of plain old 1 is because of https://github.com/JuliaArrays/BlockArrays.jl/issues/176
        mul!(X, Y, BYX, -one(T), one(T)) # X -= Y*BY'X
        # If the orthogonalization has produced results below 2eps, we drop them
        # This is to be able to orthogonalize eg [1;0] against [e^iθ;0],
        # as can happen in extreme cases in the ortho!(cP, cX)
        dropped = drop!(X)
        if dropped != []
            @views mul!(X[:, dropped], Y, BY' * (X[:, dropped]), -one(T), one(T)) # X -= Y*BY'X
        end
        if norm(BYX) < tol && niter > 1
            push!(ninners, 0)
            break
        end
        X, ninner, growth_factor = ortho!(X, tol=tol)
        push!(ninners, ninner)

        # norm(BY'X) < tol && break should be the proper check, but
        # usually unnecessarily costly. Instead, we predict the error
        # according to the growth factor of the orthogonalization.
        # Assuming BY'Y = 0 to machine precision, after the
        # Y-orthogonalization, BY'X is O(eps), so BY'X will be bounded
        # by O(eps * growth_factor).

        # If we're at a fixed point, growth_factor is 1 and if tol >
        # eps(), the loop will terminate, even if BY'Y != 0
        growth_factor*eps(real(T)) < tol && break

        niter > 10 && error("Ortho(X,Y) is failing badly, this should never happen")
        niter += 1
    end
    vprintln("ortho choleskys: ", ninners) # get how many Choleskys are performed

    # @assert (norm(BY'X)) < tol
    # @assert (norm(X'X-I)) < tol

    X
end

function final_retval(X, AX, resid_history, niter, n_matvec)
    λ = real(diag(X' * AX))
    residuals = AX .- X*Diagonal(λ)
    (λ=λ, X=X,
     residual_norms=[norm(residuals[:, i]) for i in 1:size(residuals, 2)],
     residual_history=resid_history[:, 1:niter+1],
     n_matvec=n_matvec)
end

### The algorithm is Xn+1 = rayleigh_ritz(hcat(Xn, A*Xn, Xn-Xn-1))
### We follow the strategy of Hetmaniuk and Lehoucq, and maintain a B-orthonormal basis Y = (X,R,P)
### After each rayleigh_ritz step, the B-orthonormal X and P are deduced by an orthogonal rotation from Y
### R is then recomputed, and orthonormalized explicitly wrt BX and BP
### We reuse applications of A/B when it is safe to do so, ie only with orthogonal transformations

@timing function LOBPCG(A, X, B=I, precon=I, tol=1e-10, maxiter=100;
                        miniter=1, ortho_tol=2eps(real(eltype(X))),
                        n_conv_check=nothing, display_progress=false)
    N, M = size(X)

    # If N is too small, we will likely get in trouble
    N >= 3M || @warn "Your problem is too small, and LOBPCG might
        fail; use a full diagonalization instead"

    n_conv_check === nothing && (n_conv_check = M)
    resid_history = zeros(real(eltype(X)), M, maxiter+1)

    # B-orthogonalize X
    X = ortho!(copy(X), tol=ortho_tol)[1]
    if B != I
        BX = similar(X)
        BX = mul!(BX, B, X)
        B_ortho!(X, BX)
    end

    n_matvec = M   # Count number of matrix-vector products
    AX = similar(X)
    AX = mul!(AX, A, X)
    # full_X/AX/BX will always store the full (including locked) X.
    # X/AX/BX only point to the active part
    P = zero(X)
    AP = zero(X)
    R = zero(X)
    AR = zero(X)
    if B != I
        BR = zero(X)
        # BX was already computed
        BP = zero(X)
    else
        # The B arrays are just pointers to the same data
        BR = R
        BX = X
        BP = P
    end
    nlocked = 0
    niter = 0  # the first iteration is fake
    λs = @views [(X[:,n]'*AX[:,n]) / (X[:,n]'BX[:,n]) for n=1:M]
    new_X = X
    new_AX = AX
    new_BX = BX
    full_X = X
    full_AX = AX
    full_BX = BX

    while true
        if niter > 0 # first iteration is just to compute the residuals (no X update)
            ###  Perform the Rayleigh-Ritz
            mul!(AR, A, R)
            n_matvec += size(R, 2)

            # Form Rayleigh-Ritz subspace
            if niter > 1
                Y = mortar((X, R, P))
                AY = mortar((AX, AR, AP))
                BY = mortar((BX, BR, BP))  # data shared with (X, R, P) in non-general case
            else
                Y  = mortar((X, R))
                AY = mortar((AX, AR))
                BY = mortar((BX, BR))  # data shared with (X, R) in non-general case
            end
            cX, λs = rayleigh_ritz(Y, AY, M-nlocked)

            # Update X. By contrast to some other implementations, we
            # wait on updating P because we have to know which vectors
            # to lock (and therefore the residuals) before computing P
            # only for the unlocked vectors. This results in better convergence.
            new_X  = array_mul(Y, cX)
            new_AX = array_mul(AY, cX)  # no accuracy loss, since cX orthogonal
            new_BX = (B == I) ? new_X : array_mul(BY, cX)
        end

        ### Compute new residuals
        new_R = new_AX .- new_BX .* λs'
        # it is actually a good question of knowing when to
        # precondition. Here seems sensible, but it could plausibly be
        # done before or after
        @views for i=1:size(X,2)
            resid_history[i + nlocked, niter+1] = norm(new_R[:, i])
        end
        vprintln(niter, "   ", resid_history[:, niter+1])
        if precon !== I
            @timing "preconditioning" begin
                precondprep!(precon, X) # update preconditioner if needed; defaults to noop
                ldiv!(precon, new_R)
            end
        end

        ### Compute number of locked vectors
        prev_nlocked = nlocked
        if niter ≥ miniter  # No locking if below miniter
            for i=nlocked+1:M
                if resid_history[i, niter+1] < tol
                    nlocked += 1
                    vprintln("locked $nlocked")
                else
                    # We lock in order, assuming that the lowest
                    # eigenvectors converge first; might be tricky otherwise
                    break
                end
            end
        end

        if display_progress
            println("Iter $niter, converged $(nlocked)/$(n_conv_check), resid ",
                    norm(resid_history[1:n_conv_check, niter+1]))
        end

        if nlocked >= n_conv_check  # Converged!
            X .= new_X  # Update the part of X which is still active
            AX .= new_AX
            return final_retval(full_X, full_AX, resid_history, niter, n_matvec)
        end
        newly_locked = nlocked - prev_nlocked
        active = newly_locked+1:size(X,2) # newly active vectors

        if niter > 0
            ### compute P = Y*cP only for the newly active vectors
            Xn_indices = newly_locked+1:M-prev_nlocked
            # TODO understand this for a potential save of an
            # orthogonalization, see Hetmaniuk & Lehoucq, and Duersch et. al.
            # cP = copy(cX)
            # cP[Xn_indices,:] .= 0
            e = zeros(eltype(X), size(cX, 1), M - prev_nlocked)
            for i in 1:length(Xn_indices)
                e[Xn_indices[i], i] = 1
            end
            cP = cX .- e
            cP = cP[:, Xn_indices]
            # orthogonalize against all Xn (including newly locked)
            ortho!(cP, cX, cX, tol=ortho_tol)

            # Get new P
            new_P  = array_mul( Y, cP)
            new_AP = array_mul(AY, cP)
            if B != I
                new_BP = array_mul(BY, cP)
            else
                new_BP = new_P
            end
        end

        # Update all X, even newly locked
        X  .= new_X
        AX .= new_AX
        if B != I
            BX .= new_BX
        end

        # Quick sanity check
        for i = 1:size(X, 2)
            @views if abs(BX[:, i]'X[:, i] - 1) >= sqrt(eps(real(eltype(X))))
                error("LOBPCG is badly failing to keep the vectors normalized; this should never happen")
            end
        end

        # Restrict all views to active
        @views begin
            X = X[:, active]
            AX = AX[:, active]
            BX = BX[:, active]
            R = new_R[:, active]
            AR = AR[:, active]
            BR = BR[:, active]
            P = P[:, active]
            AP = AP[:, active]
            BP = BP[:, active]
            λs = λs[active]
        end

        # Update newly active P
        if niter > 0
            P .= new_P
            AP .= new_AP
            if B != I
                BP .= new_BP
            end
        end

        # Orthogonalize R wrt all X, newly active P
        if niter > 0
            Z  = mortar((full_X, P))
            BZ = mortar((full_BX, BP))  # data shared with (full_X, P) in non-general case
        else
            Z  = full_X
            BZ = full_BX
        end
        ortho!(R, Z, BZ; tol=ortho_tol)
        if B != I
            mul!(BR, B, R)
            B_ortho!(R, BR)
        end

        niter < maxiter || break
        niter = niter + 1
    end

    final_retval(full_X, full_AX, resid_history, maxiter, n_matvec)
end


text:
##
##!------------------------------------------------------------------------------------------------!
##  generator.jl
##
##  Functions for generating graphs.
##
##  Creation date: 28/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##


text:
function addcols(model::Model, mstart::Vector{Cint}, mrwind::Vector{Cint}, dmatval::Vector{Float64}, objx::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64})
    fixinfinity!(lb)
    fixinfinity!(ub)

    newcols = Cint(length(lb))
    _chklen(ub, newcols)
    _chklen(objx, newcols)
    _chklen(mstart, newcols)
    newnz = Cint(length(dmatval))
    _chklen(mrwind, newnz)

    ret = @xprs_ccall(addcols, Cint, (
        Ptr{Nothing},    # model
        Cint,         # number of new cols
        Cint,         # number of nonzeros in added cols
        Ptr{Float64}, # coefs in objective
        Ptr{Cint},     # mstart ????
        Ptr{Cint},     # mrwind: row indexes of each col
        Ptr{Float64},      # elemnts values
        Ptr{Float64},      # lb
        Ptr{Float64}    # ub
        ),
        model.ptr_model, newcols, newnz, objx, mstart .- Cint(1), mrwind .- Cint(1), dmatval, lb, ub)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end
function addcols(model::Model, objx::Vector{Float64}, lb::Vector{Float64}, ub::Vector{Float64})
    fixinfinity!(lb)
    fixinfinity!(ub)

    newcols = Cint(length(lb))
    _chklen(ub, newcols)
    _chklen(objx, newcols)

    ret = @xprs_ccall(addcols, Cint, (
        Ptr{Nothing},    # model
        Cint,         # number of new cols
        Cint,         # number of nonzeros in added cols
        Ptr{Float64}, # coefs in objective
        Ptr{Cint},     # mstart ????
        Ptr{Cint},     # mrwind: row indexes of each col
        Ptr{Float64},      # elemnts values
        Ptr{Float64},      # lb
        Ptr{Float64}    # ub
        ),
        model.ptr_model, newcols, Cint(0), objx, inds32(newcols) .- Cint(1), C_NULL, C_NULL, lb, ub)
    if ret != 0
        throw(XpressError(model))
    end
    return nothing
end

"""
    add_var!(model::Model, mrwind::Vector, dmatval::Vector, objx::Real, lb::Real, ub::Real)
    add_var!(model::Model, vtype::Cchar, c::Real, lb::Real, ub::Real)
    add_var!(model::Model, objx::Real, lb::Real, ub::Real)
    add_var!(model::Model, vtype::Cchar, c::Real)

Add single variable with values in constraints (for column generation) value in objective `objx` with bounds `lb` and `ub`
"""
function add_var!(model::Model, nnz::Integer, mrwind::Vector, dmatval::Vector, objx::Real, lb::Real, ub::Real)

    add_var!(model::Model, mrwind::Vector, dmatval::Vector, objx::Real, lb::Real, ub::Real)
    return nothing
end
function add_var!(model::Model, mrwind::Vector, dmatval::Vector, objx::Real, lb::Real, ub::Real)
    addcols(model, Cint[1], ivec(mrwind), fvec(dmatval), fvec(objx), fvec(lb), fvec(ub))
    return nothing
end
function add_var!(model::Model, objx::Real, lb::Real, ub::Real)
    addcols(model, fvec(Float64(objx)), fvec(lb), fvec(ub))
    return nothing
end
function add_var!(model::Model, vtype::Cchar, c::Real, lb::Real, ub::Real)
    add_var!(model, c, lb, ub)
    if vtype == XPRS_INTEGER
        chgcoltype!(model, num_vars(model), XPRS_INTEGER)
    elseif vtype == XPRS_BINARY
        chgcoltype!(model, num_vars(model), XPRS_BINARY)
    end
    nothing
end
add_var!(model::Model, vtype::GChars, c::Real, lb::Real, ub::Real) = add_var!(model, cchar(vtype), Float64(c), Float64(lb), Float64(ub))
add_var!(model::Model, vtype::GChars, c::Real) = add_var!(model, vtype, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

"""
    add_cvar!(model::Model, c::Real, lb::Real, ub::Real)
    add_cvar!(model::Model, c::Real)

Add continuous variable.
"""
add_cvar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, XPRS_CONTINUOUS, c, lb, ub)
add_cvar!(model::Model, c::Real) = add_cvar!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

"""
    add_bvar!(model::Model, c::Real)

Add binary variable
"""
add_bvar!(model::Model, c::Real) = add_var!(model, XPRS_BINARY, c, 0., 1.)

"""
    add_ivar!(model::Model, c::Real, lb::Real, ub::Real)
    add_ivar!(model::Model, c::Real)

Add integer variable
"""
add_ivar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, XPRS_INTEGER, c, lb, ub)
add_ivar!(model::Model, c::Real) = add_ivar!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

"""
    add_vars!(model::Model, vtypes::Vector{Cchar}, c::Vector, lb::Vector, ub::Vector)
    add_vars!(model::Model, vtypes::GCharOrVec, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}})

Add multiple variables at once
"""
function add_vars!(model::Model, vtypes::Vector{Cchar}, c::Vector, lb::Vector, ub::Vector)
    n = length(vtypes)
    _chklen(c,n)

    prev_cols = num_vars(model)
    addcols(model, c, lb, ub)

    for i in 1:n
        if vtypes[i] == XPRS_INTEGER
            chgcoltype!(model,prev_cols+i,XPRS_INTEGER)
        elseif vtypes[i] ==XPRS_BINARY
            chgcoltype!(model,prev_cols+i,XPRS_BINARY)
        end
    end
    nothing
end
function add_vars!(model::Model, vtypes::GChars, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}}) where T<:Real
    n = length(c)
    add_vars!(model, cvecx(vtypes, n), fvec(c), fvecx(lb, n), fvecx(ub, n))
end

"""
    add_cvars!(model::Model, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}})
    add_cvars!(model::Model, c::Vector)

Add multiple continuous variables at once.
"""
add_cvars!(model::Model, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}}) where {T<:Real} = add_vars!(model, XPRS_CONTINUOUS, c, lb, ub)
add_cvars!(model::Model, c::Vector) = add_cvars!(model, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY)

"""
    add_bvars!(model::Model, c::Vector)

Add multiple binary variables at once.
"""
add_bvars!(model::Model, c::Vector) = add_vars!(model, XPRS_BINARY, c, 0, 1)

"""
    add_ivars!(model::Model, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}})
    add_ivars!(model::Model, c::Vector)

Add multiple integer variables at once.
"""
add_ivars!(model::Model, c::Vector, lb::Union{T,Vector{T}}, ub::Union{T,Vector{T}}) where {T<:Real} = add_vars!(model, XPRS_INTEGER, c, lb, ub)
add_ivars!(model::Model, c::Vector) = add_ivars!(model, XPRS_INTEGER, c, XPRS_MINUSINFINITY, XPRS_PLUSINFINITY) #

"""
    del_vars!{T<:Real}(model::Model, idx::T)
    del_vars!{T<:Real}(model::Model, idx::Vector{T})

Delete from models the varaible(s) `idx`
"""
del_vars!(model::Model, idx::T) where {T<:Real} = del_vars!(model, Cint[idx])
del_vars!(model::Model, idx::Vector{T}) where {T<:Real} = del_vars!(model, convert(Vector{Cint},idx))
function del_vars!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @xprs_ccall(delcols, Cint, (
                     Ptr{Nothing},
                     Cint,
                     Ptr{Cint}),
                     model, convert(Cint,numdel), idx .- Cint(1))
    if ret != 0
        throw(XpressError(model))
    end
end

"""
    chgcoltype!(model::Model, colnum::Integer, vtype::Cchar)

Modify type of variable `colnum` to `vtype`.
Possible values for vtype are `XPRS_CONTINUOUS`, `XPRS_INTEGER` and `XPRS_BINARY`
"""
chgcoltype!(model::Model, colnum::Integer, vtype::Cchar) = chgcoltypes!(model, Cint[colnum], Cchar[vtype])

"""
    chgcoltypes!(model::Model, colnums::Vector{Integer}, vtypes::Vector{Cchar})
    chgcoltypes!(model::Model, colnums::Vector{Integer}, vtypes::Cchar)

Modify type of variables `colnums` to `vtypes`.
Possible values for vtypes are `XPRS_CONTINUOUS`, `XPRS_INTEGER` and `XPRS_BINARY`
"""
chgcoltypes!(model::Model, colnums::Vector{Integer}, vtypes::Cchar) = chgcoltypes!(model, colnums, cvecx(vtypes))
function chgcoltypes!(model::Model, colnums::Vector{I}, vtypes::Vector{Cchar}) where I<:Integer
    n = length(colnums)
    ret = @xprs_ccall(chgcoltype,Cint,(
        Ptr{Nothing},
        Cint,
        Ptr{Cint},
        Ptr{Cchar}),
        model.ptr_model, n, ivec(colnums .- Cint(1)), vtypes )
    if 0 != ret
        throw(XpressError(model))
    end
    return nothing
end

"""
    chgsemilb!{I<:Integer, R<:Real}(model::Model, colnums::Vector{I}, lb::Vector{R})

Change values of semi continuous/integer variables `colnums` to `lb`
"""
function chgsemilb!(model::Model, colnums::Vector{I}, lb::Vector{R}) where {I<:Integer, R<:Real}
#int XPRS_CC XPRSchgglblimit(XPRSprob prob, int ncols, const int mindex[], const double dlimit[]);

    n=length(colnums)
    ret = @xprs_ccall(chgglblimit,Cint,(
        Ptr{Nothing},
        Cint,
        Ptr{Cint},
        Ptr{Float64}),
        model.ptr_model, n, ivec(colnums .- 1), fvec(lb))
    if 0 != ret
        throw(XpressError(model))
    end

end


text:
include("TwoPartCost.jl")
include("ThreePartCost.jl")
include("TechHydro.jl")
include("TechRenewable.jl")
include("TechThermal.jl")
include("Bus.jl")
include("Arc.jl")
include("Line.jl")
include("MonitoredLine.jl")
include("PhaseShiftingTransformer.jl")
include("TapTransformer.jl")
include("Transformer2W.jl")
include("HVDCLine.jl")
include("VSCDCLine.jl")
include("InterruptibleLoad.jl")
include("FixedAdmittance.jl")
include("PowerLoad.jl")
include("HydroDispatch.jl")
include("HydroFix.jl")
include("HydroStorage.jl")
include("RenewableDispatch.jl")
include("RenewableFix.jl")
include("ThermalStandard.jl")
include("LoadZones.jl")
include("GenericBattery.jl")
include("ProportionalReserve.jl")
include("StaticReserve.jl")
include("Transfer.jl")

export get_Y
export get_activepower
export get_activepower_flow
export get_activepowerlimits
export get_activepowerlimits_from
export get_activepowerlimits_to
export get_angle
export get_anglelimits
export get_arc
export get_available
export get_b
export get_basevoltage
export get_bus
export get_buses
export get_bustype
export get_capacity
export get_contributingdevices
export get_efficiency
export get_energy
export get_fixed
export get_flowlimits
export get_from
export get_fuel
export get_initial_storage
export get_inputactivepowerlimits
export get_internal
export get_inverter_firing_angle
export get_inverter_taplimits
export get_inverter_xrc
export get_loss
export get_maxactivepower
export get_maxreactivepower
export get_model
export get_name
export get_number
export get_op_cost
export get_outputactivepowerlimits
export get_powerfactor
export get_primaryshunt
export get_primemover
export get_r
export get_ramplimits
export get_rate
export get_rating
export get_reactivepower
export get_reactivepower_flow
export get_reactivepowerlimits
export get_reactivepowerlimits_from
export get_reactivepowerlimits_to
export get_rectifier_firing_angle
export get_rectifier_taplimits
export get_rectifier_xrc
export get_requirement
export get_shutdn
export get_startup
export get_storagecapacity
export get_tap
export get_tech
export get_timeframe
export get_timelimits
export get_to
export get_variable
export get_voltage
export get_voltagelimits
export get_x
export get_α


text:
# Decomente a linha abaixo para instalar os pacotes
#using Pkg; Pkg.add(["Metis","LightGraphs","GraphPlot", "Colors"])

using Metis, LightGraphs, GraphPlot, Colors

# Cria um grafo
g = LightGraphs.smallgraph(:tutte)
nodelabel = 1:nv(g)

# Faz plot inicial
p1 = gplot(g, nodelabel=nodelabel)
display(p1)

# Executa o METIS
k = 2
partition = Metis.partition(g, k, alg = :KWAY)

# Faz o plot do grafo
nodecolor = distinguishable_colors(k, colorant"blue")
nodefillc = nodecolor[partition]
p2 = gplot(g, nodefillc=nodefillc, nodelabel=nodelabel, nodelabeldist=1.8, nodelabelangleoffset=π/4)
display(p2)

max_stars_repo_name: LuizHNLorena/PFCLPMetis

max_issues_repo_name: LuizHNLorena/PFCLPMetis

max_forks_repo_name: LuizHNLorena/PFCLPMetis


text:
# ### Formatter

#md # ![FILE_NAME.png](images/FILE_NAME.png)

plot(exp, 0, 10, yformatter=:scientific)


text:
struct Model
    fw
    fc
    fs
end

function Model{T}(wordembeds::Matrix{T}, charembeds::Matrix{T}, ntags::Int)
    fw = @graph x begin
        Lookup(wordembeds)(x)
    end

    fc = @graph x begin
        d = size(charembeds, 1)
        x = Lookup(charembeds)(x)
        x = Conv1D(T,5,d,5d,2,1)(x)
        max(x, 2)
    end

    d = size(wordembeds,1) + size(charembeds,1)*5
    fs = @graph x begin
        dh = 300

        #c = Linear(T,d,dh)(x)
        #c = relu(c)
        #c = Standardize(T,(dh,10))(c)
        #f = Linear(T,d,dh)(x)
        #f = sigmoid(f)
        #x = Linear(T,d,dh)(x)
        #x = f .* c + (1-f) .* x

        #x1 = dropout(x, 0.3)
        #x1 = Conv1D(T,5,dh,dh,2,1)(x1)
        #x1 = sigmoid(x1)
        # x1 = Standardize(T,(dh,10))(x1)
        #x = x1 .* x
        x += x1

        Linear(T,dh,ntags)(x)
    end
    Model(fw, fc, fs)
end

function (m::Model)(data::Tuple)
    w, c, t = data
    w = m.fw(w)
    c = m.fc(c)
    c.batchdims = w.batchdims
    x = cat(1, w, c)
    for b = 1:length(w.batchdims)
        n = w.batchdims[b]
        for i = 1:n

        end
    end

    y = m.fs(x)
    if Merlin.config.train
        softmax_crossentropy(t, y)
    else
        vec(argmax(y.data,1))
    end
end

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# test/runtests.jl
# Test the ECOS.jl solver wrapper
#############################################################################

tests = ["direct.jl",
         "mpb_linear.jl",
         "mpb_conic.jl"]

println("Running tests:")

for curtest in tests
    println(" Test: $(curtest)")
    include(curtest)
end


text:
export AbstractSites, sitetype, data

"""
    AbstractSites{Label, T, N} <: AbstractArray{T, N}

Sites are arrays with certain labels
"""
abstract type AbstractSites{Label, T, N} <: AbstractArray{T, N} end

"""
get sitetype from subtypes of `AbstractSites`
"""
function sitetype(x::AbstractSites) end

"""
    data(sites)

get data of this `sites`
"""
function data(x::AbstractSites) end

import Compat.Random: rand, rand!, GLOBAL_RNG
export rand!, rand, flip!, randflip!

# TODO: add documents
function rand! end
function rand end
function flip! end
function randflip! end


text:
import TrajectoryOptimization.TrajectoryVariable
n,m,N = 3,2,501
X = TrajectoryVariable(N,n)
X[1] == zeros(n)
size(X) == (n,N)
to_array(X.x)

K = TrajectoryVariable(N,n,m)
K[1] == zeros(n,m)
size(K) == (n,m,N)

X = TrajectoryVariable(ones(n,N))
X *2

K = TVar1(N,(n,m),size_N=(1,1))
TVar1([Diagonal(I,4) for k = 1:N])

X = rand(n,N)
U = ones(m,N)
Xv = to_dvecs(X)
Uv = to_dvecs(U)
Xt = TVar1(Xv)
Ut = TVar1(Uv)
TX = TVar(X)
TU = TVar(U)

model,obj = Dynamics.dubinscar
solver = Solver(model,obj,N=N)
U = ones(m,N)

solver.opts.restype = Matrix
b1 = @benchmark res,stat = solve(solver,U)
norm(res.X[N] - obj.xf)

solver.opts.restype = TVar1
b_tv1 = @benchmark res_tv1,stat_tv1 = solve(solver,U)
norm(res_tv1.X[N] - obj.xf)

solver.opts.restype = TVar
b_tv = @benchmark res_tv,stat_tv = solve(solver,U)
norm(res_tv.X[N] - obj.xf)

judge(median(b_tv1),median(b1))
judge(median(b_tv),median(b1))

model,obj = Dynamics.dubinscar
n,m,N = model.n,model.m,51
x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]
obj_con = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max)
solver = Solver(model, obj_con, N=N)
U = ones(m,N-1)

solver.opts.restype = Matrix
res,stat = solve(solver,U)
solver.opts.restype = TrajectoryVariable
res_t,stat_t = solve(solver,U)
res.X[N] == res_t.X[N]

solver.opts.restype = Matrix
b1 = @benchmark solve(solver,U)
solver.opts.restype = TrajectoryVariable
b_tv1 = @benchmark solve(solver,U)
judge(median(b_tv1),median(b1))

cost(solver,X,U)
cost(solver,Xv,Uv)
cost(solver,Xt,Ut)
cost(solver,TX,TU)

b  = @benchmark(cost($solver,$X,$U))
bv = @benchmark(cost($solver,$Xv,$Ut))
bt = @benchmark(cost($solver,$Xt,$Uv))
btv = @benchmark(cost($solver,$Xt,$Uv))

judge(median(bv),median(b))
judge(median(bt),median(bv))
judge(median(btv),median(bv))

X2 = [view(X,1:n,k) for k = 1:N]
X2[1] .= zeros(n)

X[:,1]
X[:,1] = ones(n)
X2[1]

X = rand(n,N)
@btime copyto!($TX,$X)
@btime copy2!($TX,$X)
@btime copy3!($TX,$X)
@btime copyto!($Xt,$X)

@btime copyto!($TX,$TX)
@btime copyto!($Xt,$Xt)

b = rand(n)
@btime $TX[5] = $b
@btime $Xt[5] = $b

@btime $TX[5]
@btime $Xt[5]

n,m,N = 3,2,501
UnconstrainedVectorResults(n,m,N)
UnconstrainedVectorResults(n,m,N,TrajectoryVariable) isa UnconstrainedVectorResults{TV,TM} where {TV <: TrajectoryVariable, TM <: TrajectoryVariable}
ConstrainedVectorResults(n,m,N,10,4,Matrix)
ConstrainedVectorResults(n,m,N,10,4)

T = TVar1
T(N,n)

TVar(N,n)


text:
export viewImage2D

"""
	function jInv.Vis.viewImage2D

	visualizes 2D image on mesh

	Input:

	I - image data
	M - 2 dimensional AbstractTensorMesh

	kwargs get piped to pcolormesh
"""
function viewImage2D(I,M::AbstractTensorMesh;kwargs...)
	if M.dim!=2
		error("viewImage2D supports only 2-dimensional images")
	end
	I = reshape(I,tuple(M.n...))
	x1,x2 = getCellCenteredAxes(M)
	return heatmap(x1,x2,I';kwargs...)
end


text:
module DiffEqJump

using Reexport
@reexport using DiffEqBase

using Compat, RandomNumbers, TreeViews, LinearAlgebra
using DataStructures, PoissonRandom, Random, ArrayInterface
using FunctionWrappers, UnPack
using LightGraphs

import DiffEqBase: DiscreteCallback, init, solve, solve!, plot_indices, initialize!
import Base: size, getindex, setindex!, length, similar, show, merge!, merge
import DataStructures: update!
import LightGraphs: neighbors, outdegree

import RecursiveArrayTools: recursivecopy!
using StaticArrays, Base.Threads

abstract type AbstractJump end
abstract type AbstractAggregatorAlgorithm end
abstract type AbstractJumpAggregator end

import Base.Threads
@static if VERSION < v"1.3"
  seed_multiplier() = Threads.threadid()
else
  seed_multiplier() = 1
end

include("jumps.jl")
include("massaction_rates.jl")
include("aggregators/aggregators.jl")
include("aggregators/ssajump.jl")
include("aggregators/direct.jl")
include("aggregators/frm.jl")
include("aggregators/sortingdirect.jl")
include("aggregators/nrm.jl")
include("aggregators/bracketing.jl")
include("aggregators/rssa.jl")
include("aggregators/prioritytable.jl")
include("aggregators/directcr.jl")
include("aggregators/rssacr.jl")
include("aggregators/rdirect.jl")

# spatial:
include("spatial/topology.jl")
include("spatial/hop_rates.jl")
include("spatial/reaction_rates.jl")
include("spatial/flatten.jl")
include("spatial/utils.jl")

include("spatial/nsm.jl")
include("spatial/directcrdirect.jl")

include("aggregators/aggregated_api.jl")

include("extended_jump_array.jl")
include("problem.jl")
include("solve.jl")
include("coupled_array.jl")
include("coupling.jl")
include("SSA_stepper.jl")
include("simple_regular_solve.jl")

export ConstantRateJump, VariableRateJump, RegularJump, MassActionJump,
       JumpSet

export JumpProblem

export SplitCoupledJumpProblem

export Direct, DirectFW, SortingDirect, DirectCR
export BracketData, RSSA
export FRM, FRMFW, NRM
export RSSACR, RDirect

export get_num_majumps, needs_depgraph, needs_vartojumps_map

export init, solve, solve!

export reset_aggregated_jumps!

export ExtendedJumpArray

# spatial structs and functions
export CartesianGrid, CartesianGridRej, CartesianGridIter
export outdegree, num_sites, neighbors
export NSM, DirectCRDirect

end # module


text:
using Dash, DashBootstrapComponents, DashHtmlComponents

function make_item(i)
    # we use this function to make the example items to avoid code duplication
    return dbc_card([
        dbc_cardheader(
            html_h2(
                dbc_button(
                    "Collapsible group #$i",
                    color="link",
                    id="group-$i-toggle",
                    n_clicks=0,
                ),
            ),
        ),
        dbc_collapse(
            dbc_cardbody("This is the content of group $i..."),
            id="collapse-$i",
            is_open=false,
        ),
    ])
end;

accordion = html_div([make_item(1), make_item(2), make_item(3)], className="accordion");

callback!(
    app,
    Output("collapse-1", "is_open"),
    Output("collapse-2", "is_open"),
    Output("collapse-3", "is_open"),
    Input("group-1-toggle", "n_clicks"),
    Input("group-2-toggle", "n_clicks"),
    Input("group-3-toggle", "n_clicks"),
    State("collapse-1", "is_open"),
    State("collapse-2", "is_open"),
    State("collapse-3", "is_open"),
) do n1, n2, n3, is_open1, is_open2, is_open3
    ctx = callback_context()

    if length(ctx.triggered) == 0
        return false, false, false
    else
        button_id = split(ctx.triggered[1][:prop_id], ".")[1]
    end

    if button_id == "group-1-toggle" && n1 > 0
        return is_open1 == 0, false, false
    elseif button_id == "group-2-toggle" && n2 > 0
        return false, is_open2 == 0, false
    elseif button_id == "group-3-toggle" && n3 > 0
        return false, false, is_open3 == 0
    else
        return false, false, false
    end
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using LSMH
using Test, Random, BioSequences

@testset "LSMH" begin

include("murmurhash3.jl")
include("xxhash.jl")
include("hashing.jl")
include("proteins.jl")

end


text:
"""
    DataFrame(mtg::Node)
    DataFrame(mtg::Node, key)

Convert an MTG into a DataFrame.

# Arguments

- `mtg::Node`: An mtg node (usually the root node).
- `key`: The attribute(s) name(s). Select a list of variables given either as a Symbol
(faster), a String, or an Array of (or a Tuple).

# Examples

```julia
# Importing an mtg from the package:
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

# Full DataFrame:
DataFrame(mtg)

# Select just :Length:
DataFrame(mtg, :Length)

# Select just :Length and :Width:
DataFrame(mtg, [:Length, :Width])
```
"""
function DataFrame(mtg::Node, key::T) where {T<:Union{AbstractArray,Tuple}}

    tree_vars = [:node_id, :node_symbol, :node_scale, :node_index, :parent_id, :node_link] ∪ key

    # Extract the MTG info:
    nodes_info =
        traverse(
            mtg,
            node -> (
                node_id = node.id,
                node_symbol = node.MTG.symbol,
                node_scale = node.MTG.scale,
                node_index = node.MTG.index,
                parent_id = MultiScaleTreeGraph.get_parent_id(node),
                node_link = node.MTG.link
            )
        )

    # Build the DataFrame:
    df = DataFrame(nodes_info)
    insertcols!(df, 1, :tree => get_printing(mtg))

    for var in key
        insertcols!(df, var => [descendants(mtg, var, self = true)...])
    end

    # Replace the nothing values by missing values as it is the standard in DataFrames:
    for i in names(df)
        df[!, i] = replace(df[!, i], nothing => missing)
    end

    rename!(
        df,
        Dict(
            :node_id => "id",
            :node_symbol => "symbol",
            :node_scale => "scale",
            :node_index => "index",
            :node_link => "link"
        )
    )

    return df
end

function DataFrame(mtg::Node, key::T) where {T<:Symbol}
    DataFrame(mtg, [key])
end

function DataFrame(mtg::Node, key::T) where {T<:AbstractString}
    DataFrame(mtg, Symbol(key))
end

function DataFrame(mtg::Node)
    DataFrame(mtg, get_attributes(mtg))
end

function get_parent_id(x)
    if !isroot(x)
        return x.parent.id
    end
end


text:
using EntropyHub
using Test

@testset "EntropyHub.jl" begin
    #@test EntropyHub.ApEn(rand(100))
end


text:
using Pkg; Pkg.activate(@__DIR__)
# add Lyceum registry: add registry https://github.com/Lyceum/LyceumRegistry.git
using MuJoCo
mj_activate("/home/taylor/.mujoco/bin/mjkey.txt") # set location to MuJoCo key path

using LyceumMuJoCo, LyceumMuJoCoViz

# ## load MuJoCo
path = joinpath(@__DIR__, "../../../environments/atlas/deps/atlas_v5_10Hz.xml")
# path = joinpath(@__DIR__, "../../../environments/atlas/deps/atlas_v5_50Hz.xml")
# path = joinpath(@__DIR__, "../../../environments/atlas/deps/atlas_v5_100Hz.xml")
# path = joinpath(@__DIR__, "../../../environments/atlas/deps/atlas_v5_250Hz.xml")
# path = joinpath(@__DIR__, "../../../environments/atlas/deps/atlas_v5_1000Hz.xml")

include("mujoco_model.jl")
atlas = MuJoCoModel(path)
sim = LyceumMuJoCo.MJSim(atlas.m, atlas.d)

T = 10 * 1
states = Array(undef, statespace(sim), T)
# sim.d.qpos = copy(sim.m.qpos0)
sim.d.qpos[3] += 0.5 # = zero(sim.d.qpos)
@show sim.d.qpos
sim.d.qvel .= 0.0
max_viol = Float64[]
for t = 1:T
    LyceumMuJoCo.step!(sim)
    states[:, t] .= getstate(sim)
    push!(max_viol, !isempty(sim.d.contact) ? minimum([c.dist for c in sim.d.contact]) : 0.0)
end

visualize(sim, trajectories=[states])

@show minimum(max_viol)


text:
# Exercise 6-1
function compare(x, y)
    return x > y ? 1 : x == y ? 0 : -1
end

# println(compare(5, 7))
# println(compare(6, 6))
# println(compare(8, 6))

# Exercise 6-2
function hypotenuse(a, b)
    return sqrt(a^2 + b^2)
end

# println(hypotenuse(5, 8))

# Exercise 6-3
function isbetween(x, y, z)
    return x <= y <= z
end

# println(isbetween(1, 5, 10))
# println(isbetween(1, 12, 10))

# Exercise 6-4
function b(z)
    prod = a(z, z)
    println(z, " ", prod)
    prod
end

function a(x, y)
    x = x + 1
    x * y
end

function c(x, y, z)
    total = x + y + z
    square = b(total)^2
    square
end

x = 1
y = x + 1
# println(c(x, y+3, x+y))

# Exercise 6-5
function ack(m, n)
    if m == 0
        return n + 1
    elseif m > 0 && n == 0
        return ack(m - 1, 1)
    else
        return ack(m - 1, ack(m, n - 1))
    end
end

# println(ack(3, 4))

# Exercise 6-6
function first(word)
    first = firstindex(word)
    word[first]
end

function last(word)
    last = lastindex(word)
    word[last]
end

function middle(word)
    first = firstindex(word)
    last = lastindex(word)
    word[nextind(word, first) : prevind(word, last)]
end

word = "ab"
# println(middle(word))

function palindrome(word)
    if length(word) <= 1
        return true
    elseif word[begin] == word[end]
        return palindrome(middle(word))
    else
        return false
    end
end

# println(palindrome("abcdef"))
# println(palindrome("abcdeffedcba"))

# More efficient
function palindrome(word)
    l = length(word)
    for i = 1:div(l, 2)
        if word[i] != word[l - i + 1]
            return false
        end
    end
    return true
end

# Exercise 6-7
function ispower(a, b)
    if a < b
        return false
    elseif a == b
        return true
    elseif a % b == 0
        return ispower(a/b, b)
    else
        return false
    end
end

# println(ispower(32, 2))
# println(ispower(1, 1))
# println(ispower(4, 5))
# println(ispower(6, 6))
# println(ispower(32, 3))

# Exercise 6-8
function gcd(a, b)
    if b == 0
        return a
    else
        return gcd(b, a % b)
    end
end

println(gcd(5, 10))
println(gcd(5, 7))
println(gcd(12, 16))


text:
module ReinforcementLearning

export learn;

using Scenarios;

#using HDF5, JLD
using Cairo;
using DataFrames;
using Gadfly;

const SEASONS_NO = 501;
const EPISODES_NO = 400;

const EVAL_EVERY       = 250;
const EVAL_SEASONS_NO  = 20;
const EVAL_EPISODES_NO = 1000;

const EXCHANGE = false;
const LEARNS_ONE = false;

type Results{State, Action}
    bestScore::Float64
    scores::Array{Float64,2}
    statesNo::Array{Int64,2}
    drops::Array{Int64,2}
    dirty::Bool
    Qs::Array{Dict{State,Dict{Action,Float64}},1}
end

log_decay(min_val, max_val, step, steps_no) =
    max_val - log(1 + (10-1) * step / steps_no) / log(10) * (max_val - min_val);

function bestScore{State,Action}(Qs::Dict{State,Dict{Action,Float64}},
                                 state::State)
    return maximum(values(get(Qs,state,
                              Dict{Action,Float64}([(zero(Action),0.0)]))));
end

function bestAction{State,Action}(Qs::Dict{State,Dict{Action,Float64}},
                                  validActions::Array{Action,1},
                                  state::State)
    bestA = rand(validActions);
    if haskey(Qs, state)
        AsQs = Qs[state];
        bestQ = -Inf;
        for (a, Q) in AsQs
            if Q > bestQ
                bestA = a;
                bestQ = Q;
            end
        end
    end
    return bestA;
end

function bestPair{State,Action}(Qs::Dict{State,Dict{Action,Float64}},
                                validActions::Array{Action,1},
                                state::State)
    bestA = rand(validActions);
    bestQ = 0.0;
    if haskey(Qs, state)
        AsQs = Qs[state];
        for (a, Q) in AsQs
            if Q > bestQ
                bestA = a;
                bestQ = Q;
            end
        end
    end
    return (bestA, bestQ);
end

function ϵGreedy{State, Action}(Qs::Dict{State,Dict{Action,Float64}},
                                validActions::Array{Action,1},
                                state::State,
                                ϵ::Float64)
    if rand() < ϵ
        return validActions[rand(1:end)];
    else
        return bestAction(Qs, validActions, state);
    end
end

function ϵGreedy2{State, Action}(Qs::Dict{State,Dict{Action,Float64}},
                                 validActions::Array{Action, 1},
                                 state::State,
                                 ϵ::Float64,
                                 neighbourA::Action,
                                 neighbourQ::Float64)
    if rand() < ϵ
        return validActions[rand(1:end)];
    else
        A, Q = bestPair(Qs, validActions, state);
        return Q > neighbourQ ? A : neighbourA;
    end
end

function updateQs{State,Action}(Qs::Dict{State,Dict{Action,Float64}},
                                s1::State, a::Action, r::Float64, s2::State,
                                γ::Float64, α::Float64)
    AsQs = get!(Qs, s1, Dict{Action, Float64}());
    oldQ = get!(AsQs, a, 0.0);
    AsQs[a] = (1-α) * oldQ + α * (r + γ * bestScore(Qs, s2));
    nothing
end

function learn{State,Action}(s::Scenario{State,Action})
    const AGENTS_NO = s.AGENTS_NO;
    const REWARDERS_NO = s.REWARDERS_NO;

    #= Q values =#
    Qs = Array(Dict{State, Dict{Action, Float64}}, AGENTS_NO);
    for ag in 1:AGENTS_NO
        Qs[ag] = Dict{State, Dict{Action, Float64}}();
    end

    #= Actions, States, Rewards =#
    actions = Array(Action, AGENTS_NO);
    prevStates = Array(State, AGENTS_NO);
    states = Array(State, AGENTS_NO);
    #= Evaluation =#
    const evals_no = div(SEASONS_NO - 1 ,EVAL_EVERY) + 1;
    results = Results(0.0,
                      zeros(Float64, EVAL_SEASONS_NO, evals_no),
                      zeros(Int64, AGENTS_NO, evals_no),
                      zeros(Int64, 8, evals_no),
                      true,
                      deepcopy(Qs));
    # ---
    for season in 1:SEASONS_NO
        #println("S$(season)");
        gs = s.init()
        for ag in 1:AGENTS_NO
            states[ag] = s.perceive(gs, ag);
        end

        max_ϵ = log_decay(0.05, 0.3, season, SEASONS_NO);
        #total_score = 0.0;
        for episode in 1:EPISODES_NO
            # Choose actions
            ϵ = log_decay(0.003, max_ϵ, episode, EPISODES_NO);
            for ag in 1:AGENTS_NO
                if LEARNS_ONE && (mod(div(season,10), AGENTS_NO) + 1) != ag
                    actions[ag] = bestAction(Qs[ag],s.validActions,states[ag]);
                else
                    if EXCHANGE
                        ns = s.neighbours(gs, ag);
                        neighbourA = rand(s.validActions);
                        neighbourQ = 0.0;
                        for n in ns
                            nA, nQ = bestPair(Qs[n],s.validActions,states[ag]);
                            if nQ > neighbourQ
                                neighbourA = nA;
                                neighbourQ = nQ;
                            end
                        end
                        actions[ag] = ϵGreedy2(Qs[ag], s.validActions,
                                               states[ag], ϵ,
                                               neighbourA, neighbourQ);
                    else
                        actions[ag] = ϵGreedy(Qs[ag], s.validActions,
                                              states[ag], ϵ);
                    end
                end
            end

            # Do actions
            rewards, rewarders = s.doActions!(gs, actions);

            states, prevStates = prevStates, states;
            for ag in 1:AGENTS_NO
                states[ag] = s.perceive(gs, ag, prevStates[ag]);
                updateQs(Qs[ag],
                         prevStates[ag], actions[ag], rewards[ag], states[ag],
                         0.9, 0.04);
            end
        end

        if mod(season-1, EVAL_EVERY) == 0 || isfile("do_exit")
            assess(s, SEASONS_NO, EPISODES_NO, season, Qs, results);
        end
        if isfile("do_exit")
            break
        end
    end
    nothing
end

function assess{State, Action}(s::Scenario{State, Action},
                               SEASONS_NO,
                               EPISODES_NO,
                               season,
                               Qs::Array{Dict{State, Dict{Action, Float64}}, 1},
                               res::Results)
    const AGENTS_NO = s.AGENTS_NO;
    const REWARDERS_NO = s.REWARDERS_NO;

    const idx = div(season - 1, EVAL_EVERY) + 1;

    #= Number of states =#
    res.statesNo[:, idx] = [length(keys(Qs[ag])) for ag in 1:AGENTS_NO];

    #= Actions, States, Rewards =#
    actions = Array(Action, AGENTS_NO);
    prevStates = Array(State, AGENTS_NO);
    states = Array(State, AGENTS_NO);

    # ---
    for evalSeason in 1:EVAL_SEASONS_NO
        gs = s.init()
        for ag in 1:AGENTS_NO
            states[ag] = s.perceive(gs, ag);
        end

        for episode in 1:EVAL_EPISODES_NO
            # Choose actions
            for ag in 1:AGENTS_NO
                actions[ag] = bestAction(Qs[ag], s.validActions, states[ag]);
            end
            # Do actions
            rewards, rewarders = s.doActions!(gs, actions);

            prevStates, states = states, prevStates;
            for ag in 1:AGENTS_NO
                states[ag] = s.perceive(gs, ag, prevStates[ag]);
            end

            res.scores[evalSeason,idx] += sum(rewards);

            for w in 1:REWARDERS_NO
                if rewarders[w] > 0
                    res.drops[rewarders[w],idx] += 1;
                end
            end
        end
    end
        #= Save plots =#
    # Plot number of states
    @time savePlots(season, idx, res, Qs);
    if sum(res.scores[:,idx]) > res.bestScore &&
        (season > 5000 || isfile("do_exit") || isfile("do_save"))
        res.bestScore = sum(res.scores[:,idx]);
        #@time JLD.save("results/qs.jld", "Qs", Qs);
        res.Qs = deepcopy(Qs);
        res.dirty = true;
    end
    if (isfile("do_exit") || isfile("do_save") || season == SEASONS_NO) &&
        res.dirty
        @time for ag in 1:length(res.Qs)
            println(length(res.Qs[ag]));
            open("results/policy_$(ag)", "w") do f
                for (s, d) in res.Qs[ag]
                    if isempty(d) == false
                        besta = convert(Uint8, 0);
                        bestq = -Inf;
                        for (a, q) in d
                            if q > bestq
                                bestq = q
                                besta = a
                            end
                        end
                        write(f, convert(Uint8,length(s)));
                        write(f, s);
                        write(f, besta);
                    end
                end
                flush(f);
            end
        end
        open("results/scores", "w") do f
            write(f, mean(res.scores[:,1:idx], 1));
        end
        res.dirty = false;
    end
    println(season);
    flush(STDOUT);
    nothing
end

function savePlots{State,Action}(season::Int64, idx::Int64,
                                 res::Results,
                                 Qs::Array{Dict{State,Dict{Action,Float64}},1})
    dfStates = DataFrame(Season=1:EVAL_EVERY:season,
                         Max=1.0.*maximum(res.statesNo[:,1:idx],1)[:],
                         Min=1.0.*minimum(res.statesNo[:,1:idx],1)[:],
                         Mean = mean(res.statesNo[:,1:idx],1)[:]);
    statesPlot = plot(dfStates, x=:Season, y=:Mean, ymax=:Max, ymin=:Min,
                      Geom.line, Geom.ribbon,
                      Guide.ylabel("No. of states"),
                      Guide.title("Evolution of dictionary size"),
                      Guide.xlabel("Learning Season"));
    draw(PDF("results/states.pdf", 6inch, 3inch), statesPlot);
    draw(SVG("results/states.svg", 6inch, 3inch), statesPlot);
    # Plot score
    dfScore = DataFrame(Season = 1:EVAL_EVERY:season,
                        Max = 1.0 .* maximum(res.scores[:,1:idx], 1)[:],
                        Min = 1.0 .* minimum(res.scores[:,1:idx], 1)[:],
                        Mean = mean(res.scores[:,1:idx], 1)[:]);
    scorePlot = plot(dfScore, x=:Season, y=:Mean, ymax=:Max, ymin=:Min,
                     Geom.line, Geom.ribbon,
                     Guide.ylabel("Score"),
                     Guide.title("Evolution of score"),
                     Guide.xlabel("Learning Season"));

    draw(PDF("results/scores.pdf", 6inch, 3inch), scorePlot);
    draw(SVG("results/scores.svg", 6inch, 3inch), scorePlot);

    # Plot drops
    dfDrops=DataFrame(Season=[1 + div(i,8) * EVAL_EVERY for i=0:(idx*8-1)],
                      Drops = [mod(i,8)+1 for i=0:(idx*8-1)],
                      Count = res.drops[1:(idx * 8)])
    #println(dfDrops);
    const m = maximum(sum(res.drops,1));
    const m10 = max(1, div(m, 10));
    dropsPlot = plot(dfDrops, x=:Season, y=:Count, color=:Drops,
                     Geom.bar(position=:stack),
                     Guide.yticks(ticks=[0:m10:m]));
    draw(PDF("results/drops.pdf", 6inch, 3inch), dropsPlot);
    draw(SVG("results/drops.svg", 6inch, 3inch), dropsPlot);
    nothing
end

end

max_stars_repo_name: tudor-berariu/phd

max_issues_repo_name: tudor-berariu/phd

max_forks_repo_name: tudor-berariu/phd